From d69ae7d6f37abcb1bd717cbe3adb147115adb9c4 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 9 Dec 2024 14:10:25 +0100 Subject: [PATCH 01/21] Mercury module added --- bfmv5/BFMtab.xml | 36 ++++++++ src/BIO/BIO_mem.f90 | 9 ++ src/BIO/SED_mem.f90 | 38 +++++---- src/BIO/trcbio.f90 | 19 +++-- src/BIO/trcsed.f90 | 10 +-- src/General/TimeManager.f90 | 14 ++++ src/General/Timers.f90 | 1 + src/General/ogstm.f90 | 3 + src/General/step.f90 | 2 + src/IO/bc_atm_Hg0.f90 | 163 ++++++++++++++++++++++++++++++++++++ 10 files changed, 265 insertions(+), 30 deletions(-) create mode 100644 src/IO/bc_atm_Hg0.f90 diff --git a/bfmv5/BFMtab.xml b/bfmv5/BFMtab.xml index 8a7053da..34a57d87 100644 --- a/bfmv5/BFMtab.xml +++ b/bfmv5/BFMtab.xml @@ -2,6 +2,10 @@ + + + + @@ -16,33 +20,42 @@ + + + + + + + + + @@ -70,6 +83,7 @@ + @@ -78,6 +92,7 @@ + @@ -164,6 +179,16 @@ + + + + + + + + + + @@ -197,6 +222,15 @@ + + + + + + + + + @@ -236,6 +270,8 @@ + + diff --git a/src/BIO/BIO_mem.f90 b/src/BIO/BIO_mem.f90 index 45342594..101dc368 100644 --- a/src/BIO/BIO_mem.f90 +++ b/src/BIO/BIO_mem.f90 @@ -21,6 +21,7 @@ MODULE BIO_mem double precision, allocatable :: ogstm_ph(:,:,:) ! GUESS for FOLLOWS algorithm double precision, allocatable :: NPPF2(:,:,:) double precision, allocatable :: ogstm_co2(:,:), co2_IO(:,:,:) + double precision, allocatable :: atm_Hg0(:,:), Hg0_IO(:,:,:) double precision:: ice @@ -49,6 +50,12 @@ subroutine myalloc_BIO() allocate(co2_IO(jpj,jpi,2)) co2_IO = huge(co2_IO(1,1,1)) + + allocate(atm_Hg0(jpj,jpi)) + atm_Hg0 = huge(atm_Hg0(1,1)) + allocate(Hg0_IO(jpj,jpi,2)) + Hg0_IO = huge(Hg0_IO(1,1,1)) + allocate(ogstm_sedipi(jpk,jpj,jpi,4)) ogstm_sedipi = huge(ogstm_sedipi(1,1,1,1)) allocate(ogstm_ph(jpk,jpj,jpi)) @@ -74,6 +81,8 @@ subroutine clean_memory_bio() deallocate(surf_mask) deallocate(ogstm_co2) deallocate(co2_IO) + deallocate(atm_Hg0) + deallocate(Hg0_IO) deallocate(ogstm_sedipi) deallocate(ogstm_ph) deallocate(NPPF2) diff --git a/src/BIO/SED_mem.f90 b/src/BIO/SED_mem.f90 index 6b0a87a2..d1a4afcc 100644 --- a/src/BIO/SED_mem.f90 +++ b/src/BIO/SED_mem.f90 @@ -16,7 +16,7 @@ MODULE SED_mem INTEGER :: dimen_jvsed - INTEGER :: nsed=26 + INTEGER :: nsed=30 INTEGER, allocatable :: sed_idx(:) INTEGER, allocatable :: jarr_sed(:,:),jarr_sed_flx(:,:) double precision, allocatable :: ztra(:,:) @@ -64,22 +64,26 @@ subroutine myalloc_SED() sed_idx(11) = ppP1p sed_idx(12) = ppP1s sed_idx(13) = ppP1l - - sed_idx(14) = ppP2c - sed_idx(15) = ppP2n - sed_idx(16) = ppP2p - sed_idx(17) = ppP2l - - sed_idx(18) = ppP3c - sed_idx(19) = ppP3n - sed_idx(20) = ppP3p - sed_idx(21) = ppP3l - - sed_idx(22) = ppP4c - sed_idx(23) = ppP4n - sed_idx(24) = ppP4p - sed_idx(25) = ppP4l - sed_idx(26) = ppO5c + sed_idx(14) = ppP1h + + sed_idx(15) = ppP2c + sed_idx(16) = ppP2n + sed_idx(17) = ppP2p + sed_idx(18) = ppP2l + sed_idx(19) = ppP2h + + sed_idx(20) = ppP3c + sed_idx(21) = ppP3n + sed_idx(22) = ppP3p + sed_idx(23) = ppP3l + sed_idx(24) = ppP3h + + sed_idx(25) = ppP4c + sed_idx(26) = ppP4n + sed_idx(27) = ppP4p + sed_idx(28) = ppP4l + sed_idx(29) = ppP4h + sed_idx(30) = ppO5c allocate(jarr_sed(2, jpi*jpj)) jarr_sed = huge(jarr_sed(1,1)) allocate(jarr_sed_flx(jpk,jpi*jpj)) diff --git a/src/BIO/trcbio.f90 b/src/BIO/trcbio.f90 index b2af56ba..f6c452df 100644 --- a/src/BIO/trcbio.f90 +++ b/src/BIO/trcbio.f90 @@ -36,7 +36,7 @@ SUBROUTINE trcbio USE myalloc USE BIO_mem - USE OPT_mem, ONLY: PAR, RMU + USE OPT_mem, ONLY: PAR, RMU,SWR_RT USE BC_mem USE mpi @@ -61,7 +61,7 @@ SUBROUTINE trcbio double precision,dimension(jpk,jptra) :: a double precision,dimension(4,jpk) :: c double precision,dimension(jptra_dia,jpk) :: d - double precision,dimension(jpk,16) :: er + double precision,dimension(jpk,18) :: er double precision,dimension(jptra_dia_2d) :: d2 @@ -144,11 +144,12 @@ SUBROUTINE trcbio er(1:bottom,8) = PAR(1:bottom,jj,ji,3) ! PAR for pico phytoplankton er(1:bottom,9) = PAR(1:bottom,jj,ji,4) ! PAR for dinoflagellates er(1:bottom,10) = PAR(1:bottom,jj,ji,5) ! total PAR for CDOM - er(1 ,11) = DAY_LENGTH(jj,ji) ! fotoperiod expressed in hours - er(1:bottom,12) = e3t(1:bottom,jj,ji) ! depth in meters of the given cell - er(1 ,13) = vatm(jj,ji) ! wind speed (m/s) - er(1:bottom,14) = ogstm_PH(1:bottom,jj,ji) ! 8.1 - er(1 ,15) = RMU(jj,ji) ! avg. cosine direct + er(1:bottom,11) = SWR_RT(1:bottom,jj,ji) + er(1 ,12) = DAY_LENGTH(jj,ji) ! fotoperiod expressed in hours + er(1:bottom,13) = e3t(1:bottom,jj,ji) ! depth in meters of the given cell + er(1 ,14) = vatm(jj,ji) ! wind speed (m/s) + er(1:bottom,15) = ogstm_PH(1:bottom,jj,ji) ! 8.1 + er(1 ,16) = RMU(jj,ji) ! avg. cosine direct #ifndef gdept1d do jk=1,bottom @@ -161,9 +162,11 @@ SUBROUTINE trcbio correct_fact= 0.0D0 endif - er(jk,16) = correct_fact * ( gdept(jpk,jj,ji)-gdept(jk,jj,ji) ) /gdept(jpk,jj,ji) + er(jk,17) = correct_fact * ( gdept(jpk,jj,ji)-gdept(jk,jj,ji) ) /gdept(jpk,jj,ji) enddo #endif + + er(1 ,18) = atm_Hg0(jj,ji) ! ATM Hg0 call BFM1D_Input_EcologyDynamics(bottom,a,jtrmax,er) call BFM1D_reset() diff --git a/src/BIO/trcsed.f90 b/src/BIO/trcsed.f90 index 2ada62c4..883fbc47 100644 --- a/src/BIO/trcsed.f90 +++ b/src/BIO/trcsed.f90 @@ -146,32 +146,32 @@ SUBROUTINE trcsed END DO END DO ! Diatoms - DO js =9,13 + DO js =9,14 DO jk = 2,jpkm1 zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,1) * trn(jk-1,jj,ji, sed_idx(js)) END DO END DO ! Flagellates - DO js =14,17 + DO js =15,19 DO jk = 2,jpkm1 zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,2) * trn(jk-1,jj,ji, sed_idx(js)) END DO END DO ! Picophytoplankton - DO js =18,21 + DO js =20,24 DO jk = 2,jpkm1 zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,3) * trn(jk-1,jj,ji, sed_idx(js)) END DO END DO ! Dinoflagellates - DO js =22,25 + DO js =25,29 DO jk = 2,jpkm1 zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,4) * trn(jk-1,jj,ji, sed_idx(js)) END DO END DO ! Calcite - DO js =26,26 + DO js =30,30 DO jk = 2,jpkm1 zwork(jk,js,1) = - vsedO5c * trn(jk-1,jj,ji, sed_idx(js)) END DO diff --git a/src/General/TimeManager.f90 b/src/General/TimeManager.f90 index 3031bf22..3ae7261b 100644 --- a/src/General/TimeManager.f90 +++ b/src/General/TimeManager.f90 @@ -40,6 +40,7 @@ MODULE TIME_MANAGER TYPE (TIME_CONTAINER) :: TC_OPTAERO TYPE (TIME_CONTAINER) :: TC_OPTCLIM TYPE (TIME_CONTAINER) :: TC_CO2 + TYPE (TIME_CONTAINER) :: TC_Hg0 TYPE (DUMP_CONTAINER) :: RESTARTS TYPE (DUMP_CONTAINER) :: AVE_FREQ1, AVE_FREQ2 @@ -211,6 +212,9 @@ LOGICAL FUNCTION CheckStartEnd() if (.not.CheckDatelist(DATESTART,TC_CO2)) then B = .false. endif + if (.not.CheckDatelist(DATESTART,TC_Hg0)) then + B = .false. + endif if (.not.CheckDatelist(DATE__END,TC_FOR)) then B = .false. @@ -233,6 +237,10 @@ LOGICAL FUNCTION CheckStartEnd() B = .false. endif + if (.not.CheckDatelist(DATE__END,TC_Hg0)) then + B = .false. + endif + CheckStartEnd = B END FUNCTION CheckStartEnd @@ -385,6 +393,9 @@ SUBROUTINE Load_Timestrings TC_CO2%Filename = 'carbonTimes' TC_CO2%Name = 'Co2 surface value' + TC_Hg0%Filename = 'Hg0Times' + TC_Hg0%Name = 'Hg0 surface value' + RESTARTS%Filename = 'restartTimes' RESTARTS%Name = 'BFM files' @@ -410,6 +421,7 @@ SUBROUTINE Load_Timestrings call Load_Time_container(TC_FOR) ! call Load_Time_container(TC_TIN) call Load_Time_container(TC_ATM) + call Load_Time_container(TC_Hg0) ! call Load_Time_container(TC_GIB) call Load_Time_container(TC_OPTATM) call Load_Time_container(TC_OPTAERO) @@ -436,6 +448,7 @@ subroutine unload_timestrings() call unload_time_container(TC_OPTAERO) call unload_time_container(TC_OPTCLIM) call unload_time_container(TC_CO2) + call unload_time_container(TC_Hg0) call unload_dump_container(RESTARTS) call unload_dump_container(AVE_FREQ1) @@ -577,6 +590,7 @@ SUBROUTINE YEARLY(datestring) call TimeExtension(datestring,TC_OPTCLIM) call TimeExtension(datestring,TC_OPTAERO) call TimeExtension(datestring,TC_CO2) + call TimeExtension(datestring,TC_Hg0) endif END SUBROUTINE YEARLY diff --git a/src/General/Timers.f90 b/src/General/Timers.f90 index 7234220b..9ef7f936 100644 --- a/src/General/Timers.f90 +++ b/src/General/Timers.f90 @@ -8,6 +8,7 @@ MODULE Timers double precision :: forcing_phys_partTime=0.0, forcing_phys_TotTime = 0.0 double precision :: forcing_kext_partTime=0.0, forcing_kext_TotTime = 0.0 double precision :: bc_Co2_partTime =0.0, bc_Co2_TotTime = 0.0 + double precision :: bc_atm_Hg0_partTime =0.0, bc_atm_Hg0_TotTime = 0.0 double precision :: bc_atm_partTime =0.0, bc_atm_TotTime = 0.0 double precision :: bc_tin_partTime =0.0, bc_tin_TotTime = 0.0 double precision :: bc_gib_partTime =0.0, bc_gib_TotTime = 0.0 diff --git a/src/General/ogstm.f90 b/src/General/ogstm.f90 index f93cfd74..b218ee6f 100644 --- a/src/General/ogstm.f90 +++ b/src/General/ogstm.f90 @@ -328,6 +328,7 @@ SUBROUTINE time_init call TimeExtension(DATESTART,TC_OPTAERO) call TimeExtension(DATESTART,TC_OPTCLIM) call TimeExtension(DATESTART,TC_CO2) + call TimeExtension(DATESTART,TC_Hg0) call TimeInterpolation(sec,TC_FOR, TC_FOR%Before, TC_FOR%After, t_interp) @@ -337,12 +338,14 @@ SUBROUTINE time_init call TimeInterpolation(sec,TC_OPTCLIM, TC_OPTCLIM%Before, TC_OPTCLIM%After, t_interp) call TimeInterpolation(sec,TC_CO2, TC_CO2%Before, TC_CO2%After, t_interp) + call TimeInterpolation(sec,TC_Hg0, TC_Hg0%Before, TC_Hg0%After, t_interp) if (lwp) then write(*,*) 'BeforeForcings', TC_FOR%Before, 'AfterForcing', TC_FOR%After write(*,*) 'BeforeAtm', TC_ATM%Before, 'AfterAtm', TC_ATM%After write(*,*) 'BeforeCo2', TC_CO2%Before, 'AfterCo2', TC_CO2%After write(*,*) 'BeforeOpt', TC_OPTATM%Before, 'AfterOpt', TC_OPTATM%After + write(*,*) 'BeforeHg0', TC_Hg0%Before, 'AfterHg0', TC_Hg0%After endif END SUBROUTINE time_init diff --git a/src/General/step.f90 b/src/General/step.f90 index 48efea80..b4133b39 100644 --- a/src/General/step.f90 +++ b/src/General/step.f90 @@ -8,6 +8,7 @@ MODULE module_step USE mpi USE mod_atmbc USE mod_cbc + USE mod_Hgbc ! ---------------------------------------------------------------------- @@ -136,6 +137,7 @@ SUBROUTINE step CALL bc_atm (DATEstring) ! CALL dtatrc(istp,2) CALL bc_co2 (DATEstring) CALL eos () ! Water density + CALL BC_atm_Hg0 (DATEstring) diff --git a/src/IO/bc_atm_Hg0.f90 b/src/IO/bc_atm_Hg0.f90 new file mode 100644 index 00000000..a1aa47df --- /dev/null +++ b/src/IO/bc_atm_Hg0.f90 @@ -0,0 +1,163 @@ +MODULE mod_Hgbc + +USE myalloc +USE BIO_mem +USE TIME_MANAGER +USE calendar +USE mpi + +implicit NONE + +contains + + SUBROUTINE BC_atm_Hg0(datestring) + + IMPLICIT NONE + + character(LEN=17), INTENT(IN) :: datestring + +! local declarations +! ================== + double precision sec,zweigh + integer Before, After + INTEGER iswap + + bc_atm_Hg0_partTime = MPI_WTIME() + + sec=datestring2sec(DATEstring) + call TimeInterpolation(sec,TC_Hg0, BEFORE, AFTER, zweigh) + + iswap = 0 + +! ----------------------- INITIALIZATION ------------- + IF (datestring.eq.DATESTART) then +! scrivi qualche log sul numout + + CALL LOAD_atm_Hg0(TC_Hg0%TimeStrings(TC_Hg0%Before)) ! CALL trcrea(iperm1,bc) + iswap = 1 + call swap_Hg0 + + + CALL LOAD_atm_Hg0(TC_Hg0%TimeStrings(TC_Hg0%After)) ! CALL trcrea(iper,bc) + + + ENDIF + +! -------------------------------------------------------- +! and now what we have to DO at every time step +! -------------------------------------------------------- + +! check the validity of the period in memory + IF (BEFORE.ne.TC_Hg0%Before) then + TC_Hg0%Before = BEFORE + TC_Hg0%After = AFTER + + call swap_Hg0 + iswap = 1 + + + CALL LOAD_atm_Hg0(TC_Hg0%TimeStrings(TC_Hg0%After)) + + IF(lwp) WRITE (numout,*) 'Hg0 atm DATA READ for Time = ', TC_Hg0%TimeStrings(TC_Hg0%After) +! ******* LOADED NEW FRAME ************* + END IF + +! compute the DATA at the given time step + + SELECT CASE (nsptint) + CASE (0) ! ------- no time interpolation +! we have to initialize DATA IF we have changed the period + IF (iswap.eq.1) THEN + zweigh = 1.0 + call ACTUALIZE_Hg0(zweigh)! initialize now fields with the NEW DATA READ + END IF + + CASE (1) ! ------------linear interpolation --------------- + call ACTUALIZE_Hg0(zweigh) + + END SELECT + + + bc_atm_Hg0_partTime = MPI_WTIME() - bc_atm_Hg0_partTime + bc_atm_Hg0_TotTime = bc_atm_Hg0_TotTime + bc_atm_Hg0_partTime + + END SUBROUTINE BC_atm_Hg0 + + + + +! ****************************************************** +! SUBROUTINE LOAD_atm_Hg0(datestring) +! loads BC/atm_Hg0_yyyy0107-00:00:00.nc in BC_mem.Hg0_dtatrc(:,2,:) +! ****************************************************** + SUBROUTINE LOAD_atm_Hg0(datestring) + + + + IMPLICIT NONE + + CHARACTER(LEN=17), INTENT(IN) :: datestring + +! local + character(LEN=31) nomefile + + nomefile='BC/atm_Hg0_yyyy0107-00:00:00.nc' + + + ! Starting I/O + ! ********************************************************** + nomefile(1:31)='BC/atm_Hg0_'//datestring//'.nc' + if(lwp) write(*,'(A,I4,A,A)') "LOAD_atm_Hg0 --> I am ", myrank, " starting reading forcing fields from ", nomefile(1:31) + + call readnc_slice_float_2d(nomefile,'hgconc',buf2,0) + Hg0_IO(:,:,2) = buf2*tmask(1,:,:) + + + + END SUBROUTINE LOAD_atm_Hg0 + + +! ****************************************************** +! SUBROUTINE ACTUALIZE_Hg0(zweigh) +! performs time interpolation +! x(1)*(1-zweigh) + x(2)*zweigh +! ****************************************************** + SUBROUTINE ACTUALIZE_Hg0(zweigh) + + IMPLICIT NONE + + double precision, INTENT(IN) :: zweigh + + ! local + INTEGER jj,ji + + DO ji=1,jpi + DO jj=1,jpj + atm_Hg0(jj,ji) = ( (1. - zweigh) * Hg0_IO(jj,ji,1)+ zweigh * Hg0_IO(jj,ji,2) ) + END DO + END DO + + + END SUBROUTINE ACTUALIZE_Hg0 + + +! **************************************************** + + + SUBROUTINE swap_Hg0 + + IMPLICIT NONE + + + INTEGER jj,ji + + DO ji=1,jpi + DO jj=1,jpj + Hg0_IO(jj,ji,1) = Hg0_IO(jj,ji,2) + END DO + END DO + + + END SUBROUTINE swap_Hg0 + +END MODULE mod_Hgbc From eb563fb8482dac33acda8fe3fb910a6b1c288a5b Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 9 Dec 2024 14:11:34 +0100 Subject: [PATCH 02/21] set boundary conditions with Mercury --- src/IO/BC_mem.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/IO/BC_mem.f90 b/src/IO/BC_mem.f90 index 0cae6af8..5294e5c3 100644 --- a/src/IO/BC_mem.f90 +++ b/src/IO/BC_mem.f90 @@ -112,9 +112,9 @@ SUBROUTINE alloc_DTATRC ! N3n jn = 3 nitrates (climatological) ! N5s jn = 4 Silicates (seasonal) - jn_gib = 6 + jn_gib = 7 ! jn_riv = 6 - jn_atm = 2 + jn_atm = 4 ! resto is kept just to provide compliance with bfmv2, but should be removed with bfmv5 allocate(resto(jpk,jpj,jpi,jn_gib)) @@ -175,6 +175,8 @@ SUBROUTINE alloc_DTATRC tra_matrix_atm(1) = ppN1p ! phosphates tra_matrix_atm(2) = ppN3n ! nitrates + tra_matrix_atm(3) = ppHg2 ! HgII + tra_matrix_atm(4) = ppMMHg ! Mono Methyl Mercurio ENDIF From 08773a3afc2230e530bf1ae43d8fb41bf50fd4c7 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 9 Dec 2024 14:12:24 +0100 Subject: [PATCH 03/21] Adjustments for filenames (TRIM) --- src/IO/bc_atm.f90 | 4 ++-- src/IO/domrea.f90 | 6 +++--- src/IO/trcrst.f90 | 50 +++++++++++++++++++++++------------------------ src/IO/trcwri.f90 | 6 +++--- 4 files changed, 33 insertions(+), 33 deletions(-) diff --git a/src/IO/bc_atm.f90 b/src/IO/bc_atm.f90 index 9291586a..90759e6e 100644 --- a/src/IO/bc_atm.f90 +++ b/src/IO/bc_atm.f90 @@ -93,7 +93,7 @@ SUBROUTINE LOAD_ATM(datestring) CHARACTER(LEN=17), INTENT(IN) :: datestring ! local - character(LEN=7) :: nomevar + character(LEN=100) :: nomevar character(LEN=27) :: nomefile INTEGER(4) jn,jv,j,i double precision,allocatable,dimension(:,:) :: M1 @@ -112,7 +112,7 @@ SUBROUTINE LOAD_ATM(datestring) M1 = 0 nomevar = 'atm_'//ctrcnm(tra_matrix_atm(jn)) !CALL ioogsnc_bc_1d2(nomefile,nomevar,Asizeglo,atm_aux) - call readnc_slice_float_2d(nomefile,nomevar,M1,0) + call readnc_slice_float_2d(nomefile,TRIM(nomevar),M1,0) !CALL readnc_double_1d(nomefile,nomevar, Asizeglo,atm_aux) DO i=1,jpi diff --git a/src/IO/domrea.f90 b/src/IO/domrea.f90 index 6e6eda32..e241dfd1 100644 --- a/src/IO/domrea.f90 +++ b/src/IO/domrea.f90 @@ -49,8 +49,8 @@ SUBROUTINE domrea character(len=10) bfmmask_file CHARACTER(LEN=50) filename - CHARACTER(LEN=3), DIMENSION(7) :: var_nc - CHARACTER(LEN=5) nomevar01 + CHARACTER(LEN=20), DIMENSION(7) :: var_nc + CHARACTER(LEN=22) nomevar01 LOGICAL B ! ------------------- @@ -209,7 +209,7 @@ SUBROUTINE domrea IF (NWATERPOINTS.GT.0) THEN do jn=1,jn_gib - nomevar01='re'//var_nc(jn) + nomevar01='re'//TRIM(var_nc(jn)) call readnc_slice_float('bounmask.nc',nomevar01,resto(:,:,:,jn),0) enddo diff --git a/src/IO/trcrst.f90 b/src/IO/trcrst.f90 index 4e3baa29..754c7b2a 100644 --- a/src/IO/trcrst.f90 +++ b/src/IO/trcrst.f90 @@ -22,7 +22,7 @@ SUBROUTINE trcrst ! local declarations ! ================== INTEGER jn, jn_high - CHARACTER(LEN=37) filename + CHARACTER(LEN=54) filename CHARACTER(LEN=100) bkpname logical existFile logical bkp1hasbeenread,bkp2hasbeenread @@ -49,7 +49,7 @@ SUBROUTINE trcrst filename = 'RESTARTS/RST.'//DateStart//'.'//trim(ctrcnm(jn))// & '.nc' - CALL readnc_slice_double(filename, 'TRN'//trim(ctrcnm(jn)), trn(:,:,:,jn) ) + CALL readnc_slice_double(TRIM(filename), 'TRN'//trim(ctrcnm(jn)), trn(:,:,:,jn) ) @@ -67,10 +67,10 @@ SUBROUTINE trcrst if (existFilebkp) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(bkpname,ctrcnm(jn), traIO(:,:,:,jn) ) + CALL readnc_slice_double(TRIM(bkpname),TRIM(ctrcnm(jn)), traIO(:,:,:,jn) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -82,14 +82,14 @@ SUBROUTINE trcrst IF (ctr_hf(jn).eq.1) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(ctrcnm(jn))//'.nc.bkp' - INQUIRE(FILE=bkpname, EXIST=existFilebkp) + INQUIRE(FILE=TRIM(bkpname), EXIST=existFilebkp) if (existFilebkp) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(bkpname,ctrcnm(jn), traIO_HIGH(:,:,:,jn_high) ) + CALL readnc_slice_double(TRIM(bkpname),TRIM(ctrcnm(jn)), traIO_HIGH(:,:,:,jn_high) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else @@ -106,14 +106,14 @@ SUBROUTINE trcrst DO jn=1, jptra_dia bkpname = 'AVE_FREQ_2/ave.'//DateStart//'.'//trim(dianm(jn))//'.nc.bkp' - INQUIRE(FILE=bkpname, EXIST=existFile) + INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(bkpname,dianm(jn), tra_DIA_IO(jn,:,:,:) ) + CALL readnc_slice_double(TRIM(bkpname),TRIM(dianm(jn)), tra_DIA_IO(jn,:,:,:) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -123,13 +123,13 @@ SUBROUTINE trcrst IF ((diahf(jn).eq.1).and.(diaWR(jn).eq.1)) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(dianm(jn))//'.nc.bkp' - INQUIRE(FILE=bkpname, EXIST=existFile) + INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(bkpname,trim(dianm(jn)), tra_DIA_IO_HIGH(jn_high,:,:,:) ) + CALL readnc_slice_double(TRIM(bkpname),trim(dianm(jn)), tra_DIA_IO_HIGH(jn_high,:,:,:) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else @@ -146,14 +146,14 @@ SUBROUTINE trcrst DO jn=1, jptra_dia_2d bkpname = 'AVE_FREQ_2/ave.'//DateStart//'.'//trim(dianm_2d(jn))//'.nc.bkp' - INQUIRE(FILE=bkpname, EXIST=existFile) + INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double_2d(bkpname,dianm_2d(jn), tra_DIA_2d_IO(jn,:,:) ) + CALL readnc_slice_double_2d(TRIM(bkpname),dianm_2d(jn), tra_DIA_2d_IO(jn,:,:) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -163,13 +163,13 @@ SUBROUTINE trcrst IF ((diahf_2d(jn).eq.1).and.(diaWR_2d(jn).eq.1)) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(dianm_2d(jn))//'.nc.bkp' - INQUIRE(FILE=bkpname, EXIST=existFile) + INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double_2d(bkpname,trim(dianm_2d(jn)), tra_DIA_2d_IO_HIGH(jn_high,:,:) ) + CALL readnc_slice_double_2d(TRIM(bkpname),trim(dianm_2d(jn)), tra_DIA_2d_IO_HIGH(jn_high,:,:) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) - call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) + call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else diff --git a/src/IO/trcwri.f90 b/src/IO/trcwri.f90 index c5ff89e5..055977d1 100644 --- a/src/IO/trcwri.f90 +++ b/src/IO/trcwri.f90 @@ -26,9 +26,9 @@ SUBROUTINE trcwri(datestring) double precision julian - CHARACTER(LEN=37) filename + CHARACTER(LEN=54) filename - CHARACTER(LEN=3) varname + CHARACTER(LEN=20) varname INTEGER idrank, ierr, istart, jstart, iPe, iPd, jPe, jPd, status(MPI_STATUS_SIZE) INTEGER irange, jrange @@ -121,7 +121,7 @@ SUBROUTINE trcwri(datestring) filename = 'RESTARTS/RST.'//datestring//'.'//trim(var_to_store)//'.nc' - CALL write_restart(filename,var_to_store,julian, deflate_rst, deflate_level_rst) + CALL write_restart(TRIM(filename),TRIM(var_to_store),julian, deflate_rst, deflate_level_rst) END IF END IF From d0c9b518b52d275cda58c07aed0eeb45033ed737 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 9 Dec 2024 14:13:25 +0100 Subject: [PATCH 04/21] Compute shortwave for Mercury Photochemistry --- src/BIO-OPTICS/OPT_mem.f90 | 53 ++++++++- src/BIO-OPTICS/RT_exact.F90 | 209 +++++++++++++++++++++++++++++---- src/BIO-OPTICS/aasack.F90 | 21 ++-- src/BIO-OPTICS/edeseu.F90 | 7 +- src/BIO-OPTICS/trc3streams.f90 | 12 +- 5 files changed, 263 insertions(+), 39 deletions(-) diff --git a/src/BIO-OPTICS/OPT_mem.f90 b/src/BIO-OPTICS/OPT_mem.f90 index 6ef3b6af..a5e11969 100644 --- a/src/BIO-OPTICS/OPT_mem.f90 +++ b/src/BIO-OPTICS/OPT_mem.f90 @@ -95,7 +95,7 @@ MODULE OPT_mem double precision,allocatable :: Edaux(:), Esaux(:) double precision,allocatable :: cd(:,:),Cs(:,:),Bu(:,:),Cu(:,:),Bs(:,:),Fd(:,:),Bd(:,:) double precision,allocatable :: au(:,:),as(:,:),bquad(:,:),cquad(:,:),sqarg(:,:) - double precision,allocatable :: inhoD(:,:),inhox(:,:),inhoy(:,:) + double precision,allocatable :: inhoD(:,:),inhox(:,:),inhoy(:,:),inhoD_r(:,:) double precision,allocatable :: D(:,:),a_m(:,:),a_p(:,:) double precision,allocatable :: r_m(:,:),r_p(:,:) double precision,allocatable :: e_m(:,:),e_p(:,:) @@ -113,6 +113,7 @@ MODULE OPT_mem ! Outputs of radiative transfer model double precision,allocatable :: Ed(:,:,:,:), Es(:,:,:,:), Eu(:,:,:,:) ! depth, lat, lon, wave length double precision,allocatable :: PAR(:,:,:,:) ! depth, lat, lon, phyto + double precision,allocatable :: SWR_RT(:,:,:) ! depth, lat, lon double precision,allocatable :: RMU(:,:) ! lat, lon double precision,allocatable :: Ed_IO(:,:,:,:),Es_IO(:,:,:,:),Eu_IO(:,:,:,:) logical, allocatable, dimension(:,:,:) :: opt_mask @@ -210,28 +211,75 @@ subroutine myalloc_OPT() call lidata() allocate(Edaux(nlt)) + Edaux=huge(Edaux(1)) allocate(Esaux(nlt)) + Esaux=huge(Esaux(1)) allocate(cd(jpk,nlt),Cs(jpk,nlt),Bu(jpk,nlt),Cu(jpk,nlt),Bs(jpk,nlt),Fd(jpk,nlt),Bd(jpk,nlt)) + cd=huge(cd(1,1)) + Cs=huge(Cs(1,1)) + Bu=huge(Bu(1,1)) + Cu=huge(Cu(1,1)) + Bs=huge(Bs(1,1)) + Fd=huge(Fd(1,1)) + Bd=huge(Bd(1,1)) allocate(au(jpk,nlt),as(jpk,nlt),bquad(jpk,nlt),cquad(jpk,nlt),sqarg(jpk,nlt)) - allocate(inhoD(jpk,nlt),inhox(jpk,nlt),inhoy(jpk,nlt)) + au=huge(au(1,1)) + as=huge(as(1,1)) + bquad=huge(bquad(1,1)) + cquad=huge(cquad(1,1)) + sqarg=huge(sqarg(1,1)) + allocate(inhoD(jpk,nlt),inhox(jpk,nlt),inhoy(jpk,nlt),inhoD_r(jpk,nlt)) + inhoD = huge(inhoD(1,1)) + inhoD_r = huge(inhoD_r(1,1)) + inhox = huge(inhox(1,1)) + inhoy = huge(inhoy(1,1)) allocate(D(jpk,nlt),a_m(jpk,nlt),a_p(jpk,nlt)) + D = huge(D(1,1)) + a_m = huge(a_m(1,1)) + a_p = huge(a_p(1,1)) allocate(r_m(jpk,nlt),r_p(jpk,nlt)) + r_m = huge(r_m(1,1)) + r_p = huge(r_p(1,1)) allocate(e_m(jpk,nlt),e_p(jpk,nlt)) + e_m = huge(e_m(1,1)) + e_p = huge(e_p(1,1)) allocate(zeta0(nlt),eta0(nlt)) + zeta0 = huge(zeta0(1)) + eta0 = huge(eta0(1)) allocate(alpha(jpk-1,nlt),beta(jpk-1,nlt),gamm(jpk-1,nlt),delta(jpk-1,nlt)) + alpha = huge(alpha(1,1)) + beta = huge(beta(1,1)) + gamm = huge(gamm(1,1)) + delta = huge(delta(1,1)) allocate(epsRT(jpk-1,nlt),zeta(jpk-1,nlt),eta(jpk-1,nlt),theta(jpk-1,nlt)) + epsRT = huge(epsRT(1,1)) + zeta = huge(zeta(1,1)) + eta = huge(eta(1,1)) + theta = huge(theta(1,1)) ! allocate(vD(2*jpk-1,nlt),vL(2*jpk-1,nlt),vU(2*jpk-1,nlt)) ! allocate(WW(2*jpk-1,nlt), WW1(2*jpk-1,nlt)) ! allocate(sol(2*jpk-1,nlt),sol_p(jpk,nlt),sol_m(jpk,nlt)) allocate(err_RT(nlt)) + err_RT = huge(err_RT(1)) ! Additional variables for approximate model allocate(a1(jpk,nlt),a2(jpk,nlt),S(jpk,nlt)) + a1 = huge(a1(1,1)) + a2 = huge(a2(1,1)) + S = huge(S(1,1)) allocate(SEdz(jpk,nlt),a2ma1(jpk,nlt),rM(jpk,nlt),rN(jpk,nlt),c2(jpk,nlt)) + SEdz = huge(SEdz(1,1)) + a2ma1 = huge(a2ma1(1,1)) + rM = huge(rM(1,1)) + rN = huge(rN(1,1)) + c2 = huge(c2(1,1)) allocate(Ta2z(jpk,nlt), Eutmp(jpk,nlt)) + Ta2z = huge(Ta2z(1,1)) + Eutmp = huge(Eutmp(1,1)) ! Allocate output variables allocate(Ed(jpk,jpj,jpi,nlt),Es(jpk,jpj,jpi,nlt),Eu(jpk,jpj,jpi,nlt)) allocate(PAR(jpk,jpj,jpi,nchl+1)) ! last index total par + allocate(SWR_RT(jpk,jpj,jpi)) allocate(RMU(jpj,jpi)) allocate(Ed_IO(jpk,jpj,jpi,nlt),Es_IO(jpk,jpj,jpi,nlt),Eu_IO(jpk,jpj,jpi,nlt)) allocate(opt_mask(jpk,jpj,jpi)) ; opt_mask=0 @@ -242,6 +290,7 @@ subroutine myalloc_OPT() Es(:,:,:,:) = 1.e-4 Eu(:,:,:,:) = 1.e-8 PAR(:,:,:,:) = 0.0d0 + SWR_RT(:,:,:) = 0.0d0 RMU(:,:) = 0.0d0 Ed_IO(:,:,:,:) = 1.e-4 Es_IO(:,:,:,:) = 1.e-4 diff --git a/src/BIO-OPTICS/RT_exact.F90 b/src/BIO-OPTICS/RT_exact.F90 index 9c0cacd6..527005ba 100644 --- a/src/BIO-OPTICS/RT_exact.F90 +++ b/src/BIO-OPTICS/RT_exact.F90 @@ -1,4 +1,4 @@ - subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0,PARz) + subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0,PARz,SWRz) USE OPT_mem USE Tridiagonal IMPLICIT NONE @@ -20,19 +20,41 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, CHARACTER(*), INTENT(IN) :: V_POSITION double precision, INTENT(IN) :: Edtop(nlt),Estop(nlt) double precision, INTENT(OUT) :: Edz(jpk,nlt),Esz(jpk,nlt),Euz(jpk,nlt),Eu0(nlt) - double precision, INTENT(OUT) :: PARz(jpk,nchl+1) + double precision, INTENT(OUT) :: PARz(jpk,nchl+1),SWRz(jpk) ! local variables double precision, parameter :: rd=1.5d0 !these are taken from Ackleson, et al. 1994 (JGR) double precision, parameter :: ru=3.0d0 integer :: k,ii,b,p,nl + logical :: is_Singular(jpk,nlt) + CHARACTER(20) :: V_POS double precision :: rmus, rmuu + double precision :: RT_depth(jpk) double precision :: Edzm(jpk,nlt),Edza(jpk,nlt) double precision :: vD(2*bottom-1,nlt),vL(2*bottom-1,nlt),vU(2*bottom-1,nlt) double precision :: WW(2*bottom-1,nlt),WW1(2*bottom-1,nlt),sol(2*bottom-1,nlt) double precision :: sol_p(2*bottom-1,nlt), sol_m(2*bottom-1,nlt) double precision :: aux1(nlt), aux2(nlt) + + Edzm=huge(Edzm(1,1)) + Edza=huge(Edza(1,1)) + vD=huge(vD(1,1)) + vL=huge(vL(1,1)) + vU=huge(vU(1,1)) + WW=huge(WW(1,1)) + WW1=huge(WW1(1,1)) + sol=huge(sol(1,1)) + sol_p=huge(sol_p(1,1)) + sol_m=huge(sol_m(1,1)) + aux1=huge(aux1(1)) + aux2=huge(aux2(1)) + Edz=huge(Edz(1,1)) + Esz=huge(Esz(1,1)) + Euz=huge(Euz(1,1)) + Eu0=huge(Eu0(1)) + V_POS = V_POSITION + ! Constants rmus = 1.0d0/0.83d0 !avg cosine diffuse down rmuu = 1.0d0/0.4d0 !avg cosine diffuse up @@ -51,13 +73,38 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, ! ! Solution of the inhomogenous system (inhox and inhoy) - inhoD(1:b,:) = 1.0d0/((cd(1:b,:)-Cs(1:b,:))*(cd(1:b,:)+Cu(1:b,:))+Bs(1:b,:)*Bu(1:b,:)) - inhox(1:b,:) = inhoD(1:b,:)*(- ( cd(1:b,:)+Cu(1:b,:) )*Fd(1:b,:) -Bu(1:b,:)*Bd(1:b,:) ) - inhoy(1:b,:) = inhoD(1:b,:)*(- Fd(1:b,:)*Bs(1:b,:) + ( cd(1:b,:)-Cs(1:b,:) )*Bd(1:b,:) ) +! inhoD(1:b,:) = 1.0d0/((cd(1:b,:)-Cs(1:b,:))*(cd(1:b,:)+Cu(1:b,:))+Bs(1:b,:)*Bu(1:b,:)) +! inhox(1:b,:) = inhoD(1:b,:)*(- ( cd(1:b,:)+Cu(1:b,:) )*Fd(1:b,:) -Bu(1:b,:)*Bd(1:b,:) ) +! inhoy(1:b,:) = inhoD(1:b,:)*(- Fd(1:b,:)*Bs(1:b,:) + ( cd(1:b,:)-Cs(1:b,:) )*Bd(1:b,:) ) ! write(*,*) 'inhoD',inhoD(1,1) ! write(*,*) 'inhox',inhox(1,1) ! write(*,*) 'inhoy',inhoy(1,1) + do nl=1,nlt + do k =1,bottom + inhoD(k,nl) = ((cd(k,nl)-Cs(k,nl))*(cd(k,nl)+Cu(k,nl))+Bs(k,nl)*Bu(k,nl)) + if (abs(inhoD(k,nl)) .LT. 1.0D-4) then + is_Singular(k,nl) = .True. + V_POS="MID" + else + is_Singular(k,nl) = .False. + endif + enddo + enddo + do nl=1,nlt + do k =1,bottom + if (is_Singular(k,nl) ) then + inhox(k,nl) = Fd(k,nl) + inhoy(k,nl) = -Bd(k,nl) + else + inhoD_r(k,nl) = 1.0d0/inhoD(k,nl) + inhox(k,nl) = inhoD_r(k,nl)*(- ( cd(k,nl)+Cu(k,nl) )*Fd(k,nl) -Bu(k,nl)*Bd(k,nl) ) + inhoy(k,nl) = inhoD_r(k,nl)*(- Fd(k,nl)*Bs(k,nl) + ( cd(k,nl)-Cs(k,nl) )*Bd(k,nl) ) + endif + enddo + enddo + + ! eigenvalues of the homogeneous system D(1:b,:) = 0.5d0*( Cs(1:b,:)+Cu(1:b,:) & + sqrt((Cs(1:b,:)+Cu(1:b,:))*(Cs(1:b,:)+Cu(1:b,:))-4.0d0*Bs(1:b,:)*Bu(1:b,:) ) ) @@ -87,8 +134,16 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, if (k > 1) then Edaux(:) = Edz(k-1,:) endif - Edz(k,:) = Edaux(:)*exp(-cd(k,:)*zd(k)) - Edzm(k,:) = Edaux(:)*exp(-0.5D0*cd(k,:)*zd(k)) ! mid cell +! Edz(k,:) = Edaux(:)*exp(-cd(k,:)*zd(k)) +! Edzm(k,:) = Edaux(:)*exp(-0.5D0*cd(k,:)*zd(k)) ! mid cell + + Edz(k,:) = Edaux(:)*exp(-cd(k,:)*zd(k)) + Edzm(k,:) = Edaux(:)*exp(-0.5D0*cd(k,:)*zd(k)) ! mid cell + if (k .EQ. 1) then + RT_depth(k) = zd(k) + else + RT_depth(k) = RT_depth(k-1)+zd(k) + endif enddo ! Imposing boundary conditions @@ -101,17 +156,60 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, alpha(k,:) = e_p(k,:)*(1.0d0-r_p(k,:)*r_m(k+1,:) ) beta(k,:) = r_m(k,:) - r_m(k+1,:) gamm(k,:) = - (1.0d0 -r_p(k+1,:)*r_m(k+1,:) ) - delta(k,:) = ( inhox(k+1,:) - inhox(k,:) - (inhoy(k+1,:) - inhoy(k,:) ) & - * r_m(k+1,:) ) * Edz(k,:) !! Check if Edz(k) is correct +! delta(k,:) = ( inhox(k+1,:) - inhox(k,:) - (inhoy(k+1,:) - inhoy(k,:) ) & +! * r_m(k+1,:) ) * Edz(k,:) !! Check if Edz(k) is correct epsRT(k,:) = 1.0d0-r_m(k,:)*r_p(k,:) zeta(k,:) = -(r_p(k+1,:)-r_p(k,:)) eta(k,:) = -e_m(k+1,:)*(1.0d0-r_p(k,:)*r_m(k+1,:) ) - theta(k,:) = ( inhoy(k+1,:) - inhoy(k,:) - (inhox(k+1,:) - inhox(k,:) ) & - * r_p(k,:) ) * Edz(k,:) !! Check if Edz(k) is correct +! theta(k,:) = ( inhoy(k+1,:) - inhoy(k,:) - (inhox(k+1,:) - inhox(k,:) ) & +! * r_p(k,:) ) * Edz(k,:) !! Check if Edz(k) is correct enddo +! Theta and Delta depends on the inhomogenous solution and their determination will be different according +! to presence of singularities + do nl=1,nlt + do k =1,bottom-1 +! k Singular +! - +! k+1 Singular + if (is_Singular(k,nl) .AND. is_Singular(k+1,nl)) then + delta(k,nl) = ( inhox(k+1,nl) - inhox(k,nl) - (inhoy(k+1,nl) - inhoy(k,nl) ) & + * r_m(k+1,nl) ) * RT_depth(k) * Edz(k,nl) !! + theta(k,nl) = ( inhoy(k+1,nl) - inhoy(k,nl) - (inhox(k+1,nl) - inhox(k,nl) ) & + * r_p(k,nl) ) * RT_depth(k) * Edz(k,nl) !! + endif +! k Non Singular +! - +! k+1 Singular + if( (.NOT. is_Singular(k,nl)) .AND. (is_Singular(k+1,nl))) then + delta(k,nl) = ( RT_depth(k) * inhox(k+1,nl) - inhox(k,nl) - (RT_depth(k) *inhoy(k+1,nl) - inhoy(k,nl) ) & + * r_m(k+1,nl) ) * Edz(k,nl) !! + theta(k,nl) = ( RT_depth(k) * inhoy(k+1,nl) - inhoy(k,nl) - (RT_depth(k) *inhox(k+1,nl) - inhox(k,nl) ) & + * r_p(k,nl) ) * Edz(k,nl) !! + endif +! k Singular +! - +! k+1 Non Singular + if( (is_Singular(k,nl)) .AND. ( .NOT. is_Singular(k+1,nl))) then + delta(k,nl) = ( inhox(k+1,nl) - RT_depth(k) *inhox(k,nl) - (inhoy(k+1,nl) - RT_depth(k) *inhoy(k,nl) ) & + * r_m(k+1,nl) ) * Edz(k,nl) !! + theta(k,nl) = ( inhoy(k+1,nl) - RT_depth(k) *inhoy(k,nl) - (inhox(k+1,nl) - RT_depth(k) *inhox(k,nl) ) & + * r_p(k,nl) ) * Edz(k,nl) !! + endif +! k Non Singular +! - +! k+1 Non Singular + if (( .NOT. is_Singular(k,nl)) .AND. ( .NOT. is_Singular(k+1,nl))) then + delta(k,nl) = ( inhox(k+1,nl) - inhox(k,nl) - (inhoy(k+1,nl) - inhoy(k,nl) ) & + * r_m(k+1,nl) ) * Edz(k,nl) + theta(k,nl) = ( inhoy(k+1,nl) - inhoy(k,nl) - (inhox(k+1,nl) - inhox(k,nl) ) & + * r_p(k,nl) ) * Edz(k,nl) + endif + enddo + enddo + ! contruction of the vectors of the tri-diagonal matrix ! diagonal vector vD @@ -148,7 +246,14 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, ! constant vector - WW(1,:)= Estop(:) - inhox(1,:) * Edtop(:) +! WW(1,:)= Estop(:) - inhox(1,:) * Edtop(:) + do nl=1,nlt + if (is_Singular(1,nl) ) then + WW(1,nl)= Estop(nl) + else + WW(1,nl)= Estop(nl) - inhox(1,nl) * Edtop(nl) + endif + enddo do k =1,bottom-1 ii = 2*k @@ -234,16 +339,27 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, enddo CASE ("MID") - + do nl=1,nlt do k =1,bottom - aux1(:) = exp(-a_p(k,:)*zd(k)*0.5D0) - aux2(:) = exp(-a_m(k,:)*zd(k)*0.5D0) - Edz(k,:) = Edzm(k,:) - Esz(k,:) = sol_p(k,:) * aux1(:) + sol_m(k,:)*r_m(k,:)*aux2(:) + inhox(k,:)*Edzm(k,:) - Euz(k,:) = sol_p(k,:)*r_p(k,:) * aux1(:) + sol_m(k,:)*aux2(:) + inhoy(k,:)*Edzm(k,:) +! aux1(:) = exp(-a_p(k,:)*zd(k)*0.5D0) +! aux2(:) = exp(-a_m(k,:)*zd(k)*0.5D0) +! Edz(k,:) = Edzm(k,:) +! Esz(k,:) = sol_p(k,:) * aux1(:) + sol_m(k,:)*r_m(k,:)*aux2(:) + inhox(k,:)*Edzm(k,:) +! Euz(k,:) = sol_p(k,:)*r_p(k,:) * aux1(:) + sol_m(k,:)*aux2(:) + inhoy(k,:)*Edzm(k,:) + aux1(nl) = exp(-a_p(k,nl)*zd(k)*0.5D0) + aux2(nl) = exp(-a_m(k,nl)*zd(k)*0.5D0) + Edz(k,nl) = Edzm(k,nl) + if (is_Singular(k,nl) ) then + Esz(k,nl) = sol_p(k,nl) * aux1(nl) + sol_m(k,nl)*r_m(k,nl)*aux2(nl) + inhox(k,nl)*Edzm(k,nl) * (RT_depth(k)-zd(k)*0.5D0) + Euz(k,nl) = sol_p(k,nl)*r_p(k,nl) * aux1(nl) + sol_m(k,nl) *aux2(nl) + inhoy(k,nl)*Edzm(k,nl) * (RT_depth(k)-zd(k)*0.5D0) + else + Esz(k,nl) = sol_p(k,nl) * aux1(nl) + sol_m(k,nl)*r_m(k,nl)*aux2(nl) + inhox(k,nl)*Edzm(k,nl) + Euz(k,nl) = sol_p(k,nl)*r_p(k,nl) * aux1(nl) + sol_m(k,nl) *aux2(nl) + inhoy(k,nl)*Edzm(k,nl) + endif enddo + enddo CASE ("AVERAGE") @@ -273,11 +389,55 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, ! surface upward diffuse radiance at depth = 0 - Eu0(:) = sol_p(1,:)*r_p(1,:) + sol_m(1,:)*exp(-a_m(1,:)*zd(1)) + inhoy(1,:) * Edtop(:) +! Eu0(:) = sol_p(1,:)*r_p(1,:) + sol_m(1,:)*exp(-a_m(1,:)*zd(1)) + inhoy(1,:) * Edtop(:) ! write(*,*) 'Eu', Eu - + do nl=1,nlt + if (is_Singular(1,nl) ) then + Eu0(:) = sol_p(1,:)*r_p(1,:) + sol_m(1,:)*exp(-a_m(1,:)*zd(1)) + else + Eu0(:) = sol_p(1,:)*r_p(1,:) + sol_m(1,:)*exp(-a_m(1,:)*zd(1)) + inhoy(1,:) * Edtop(:) + endif + enddo + do nl=1,nlt + if (abs(Eu0(nl)) .GT. 1000.) then + write(1100,*) 'size PARz', size(PARz) + write(1101,*) 'V_POSITION', V_POSITION + write(1102,*) 'bottom', bottom + write(1103,*) 'zd', zd + write(1104,*) 'Edtop', Edtop + write(1105,*) 'Estop', Estop + write(1106,*) 'rmud', rmud + write(1107,*) 'a', a(:,nl) + write(1108,*) 'bt', bt(:,nl) + write(1109,*) 'bb', bb(:,nl) + write(1110,*) 'Edz', Edz(:,nl) + write(1111,*) 'Esz', Esz(:,nl) + write(1112,*) 'Euz', Euz(:,nl) + write(1113,*) 'Eu0', Eu0(:) + write(1114,*) 'PARz', PARz + write(1115,*) 'SWRz', SWRz + +! call SLAE3diag(2*bottom-1, nlt , vL, vD, vU, WW, sol) + write(1116,*) 'vL', vL(:,nl) + write(1117,*) 'vD', vL(:,nl) + write(1118,*) 'vU', vL(:,nl) + write(1119,*) 'WW', WW(:,nl) + write(1120,*) 'sol', sol(:,nl) + write(1121,*) 'cd', cd(:,nl) + write(1122,*) 'Cs', Cs(:,nl) + write(1123,*) 'Cu', Cu(:,nl) + write(1124,*) 'Bs', Bs(:,nl) + write(1125,*) 'Bu', Bu(:,nl) + write(1126,*) 'Fd', Fd(:,nl) + write(1127,*) 'Bd', Bd(:,nl) + write(1128,*) 'inhoD', inhoD(:,nl) + write(1129,*) 'inhox', inhox(:,nl) + STOP + endif + enddo ! compute PAR for each PFT using scalar irradiance PARz(:,:) = 0.0D0 + SWRz(:) = 0.0D0 do p =1,nchl do nl =1,nlt ! 400 to 700 nm @@ -296,6 +456,13 @@ subroutine radmod_ex(V_POSITION, bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz, enddo enddo - + + do nl =1,nlt ! from 0 to 4 um + do k =1,bottom + SWRz(k) = SWRz(k) +( Edz(k,nl) * rmud + Esz(k,nl) * rmus + Euz(k,nl) * rmuu ) + enddo + enddo + + return end diff --git a/src/BIO-OPTICS/aasack.F90 b/src/BIO-OPTICS/aasack.F90 index aa52d60b..5940566b 100644 --- a/src/BIO-OPTICS/aasack.F90 +++ b/src/BIO-OPTICS/aasack.F90 @@ -1,4 +1,4 @@ - subroutine radmod(V_POSITION,bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0,PARz) + subroutine radmod(V_POSITION,bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0,PARz,SWRz) USE OPT_mem IMPLICIT NONE @@ -27,6 +27,7 @@ subroutine radmod(V_POSITION,bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0, character(*), INTENT(IN) :: V_POSITION double precision, INTENT(OUT) :: Edz(jpk,nlt),Esz(jpk,nlt),Euz(jpk,nlt),Eu0(nlt) double precision, INTENT(OUT) :: PARz(jpk,nchl+1) + double precision, INTENT(OUT) :: SWRz(jpk) ! local variables double precision, parameter :: rd=1.5d0 !these are taken from Ackleson, et al. 1994 (JGR) @@ -197,25 +198,25 @@ subroutine radmod(V_POSITION,bottom,zd,Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eu0, do p =1,nchl do nl =5,19 ! 400 to 700 nm - - PARz(1,p) = PARz(1,p) + WtoQ(nl) * ac(p,nl) * ( Edz(1,nl) * rmud + Esz(1,nl) * rmus) - - do jk =2,bottom + do jk =1,bottom PARz(jk,p) = PARz(jk,p) + WtoQ(nl) * ac(p,nl) * ( Edz(jk,nl) * rmud + Esz(jk,nl) * rmus + Euz(jk,nl) * rmuu ) enddo enddo enddo - - do nl =5,19 ! 400 to 700 nm - - PARz(1,nchl+1) = PARz(1,nchl+1) + WtoQ(nl) * ( Edz(1,nl) * rmud + Esz(1,nl) * rmus ) - do jk =2,bottom + do nl =5,19 ! 400 to 700 nm + do jk =1,bottom PARz(jk,nchl+1) = PARz(jk,nchl+1) + WtoQ(nl) * ( Edz(jk,nl) * rmud + Esz(jk,nl) * rmus + Euz(jk,nl) * rmuu ) enddo enddo + do nl =1,nlt ! 0 to 4 um + do jk =1,bottom + SWRz(jk) = SWRz(jk) + Edz(jk,nl) * rmud + Esz(jk,nl) * rmus + Euz(jk,nl) * rmuu + enddo + enddo + return end diff --git a/src/BIO-OPTICS/edeseu.F90 b/src/BIO-OPTICS/edeseu.F90 index 6c83e48a..72de8fc6 100644 --- a/src/BIO-OPTICS/edeseu.F90 +++ b/src/BIO-OPTICS/edeseu.F90 @@ -1,4 +1,4 @@ - subroutine edeseu(MODE,V_POSITION,bottom,dzRT,Edtop,Estop,CHLz,PCz,CDOMz,POCz,rmud,Edz,Esz,Euz,Eutop,PARz) + subroutine edeseu(MODE,V_POSITION,bottom,dzRT,Edtop,Estop,CHLz,PCz,CDOMz,POCz,rmud,Edz,Esz,Euz,Eutop,PARz,SWRz) USE myalloc USE mpi USE OPT_mem @@ -24,6 +24,7 @@ subroutine edeseu(MODE,V_POSITION,bottom,dzRT,Edtop,Estop,CHLz,PCz,CDOMz,POCz,rm double precision, INTENT(OUT) :: Edz(jpk,nlt),Esz(jpk,nlt) double precision, INTENT(OUT) :: Euz(jpk,nlt) double precision, INTENT(OUT) :: PARz(jpk,nchl+1) + double precision, INTENT(OUT) :: SWRz(jpk) double precision, intent(OUT) :: Eutop(nlt) ! Local variables double precision :: Etop @@ -63,11 +64,11 @@ subroutine edeseu(MODE,V_POSITION,bottom,dzRT,Edtop,Estop,CHLz,PCz,CDOMz,POCz,rm enddo if (MODE == 1) then ! condition on level < 200 mt - call radmod(V_POSITION,bottom,dzRT(:),Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eutop,PARz) + call radmod(V_POSITION,bottom,dzRT(:),Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eutop,PARz,SWRz) endif ! condition on level > 200 mt if ( MODE ==0) then - call radmod_ex(V_POSITION,bottom,dzRT(:),Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eutop,PARz) + call radmod_ex(V_POSITION,bottom,dzRT(:),Edtop,Estop,rmud,a,bt,bb,Edz,Esz,Euz,Eutop,PARz,SWRz) endif diff --git a/src/BIO-OPTICS/trc3streams.f90 b/src/BIO-OPTICS/trc3streams.f90 index 01a40421..c123437d 100644 --- a/src/BIO-OPTICS/trc3streams.f90 +++ b/src/BIO-OPTICS/trc3streams.f90 @@ -26,10 +26,10 @@ SUBROUTINE trc3streams(datestring) INTEGER :: it_actual CHARACTER(LEN=20) :: V_POSITION double precision :: solz(jpj,jpi), rmud(jpj,jpi) - double precision :: Edz(jpk,nlt),Esz(jpk,nlt),Euz(jpk,nlt) + double precision :: Edz(jpk,nlt),Esz(jpk,nlt),Euz(jpk,nlt),SWRz(jpk) double precision, allocatable:: E(:,:,:) !(3,jpk+1,nlt) - double precision, allocatable :: PARz(:,:) !(jpk,nchl+1) + double precision, allocatable :: PARz(:,:) !, SWRz(:) !(jpk,nchl+1) double precision :: CHLz(jpk,nchl),PCz(jpk,nchl),CDOMz(jpk,3),POCz(jpk) double precision :: Eu_0m(nlt) double precision :: zgrid(jpk+1) @@ -96,6 +96,7 @@ SUBROUTINE trc3streams(datestring) Es(1:bottom,jj,ji,:) = 0.0001d0 Eu(2:bottom,jj,ji,:) = 1.0E-08 PAR(1:bottom,jj,ji,:) = 0.0001d0 + SWR_RT(1:bottom,jj,ji) = 0.0001d0 else @@ -116,7 +117,7 @@ SUBROUTINE trc3streams(datestring) POCz(1:bottom) = trn(1:bottom,jj,ji,ppR6c) IF ( (MODE .EQ. 0) .OR. (MODE .EQ. 1)) then - call edeseu(MODE,V_POSITION,bottom,e3t(:,jj,ji),Ed_0m(:,jj,ji),Es_0m(:,jj,ji),CHLz,PCz,CDOMz,POCz,rmud(jj,ji),Edz,Esz,Euz,Eu_0m,PARz) + call edeseu(MODE,V_POSITION,bottom,e3t(:,jj,ji),Ed_0m(:,jj,ji),Es_0m(:,jj,ji),CHLz,PCz,CDOMz,POCz,rmud(jj,ji),Edz,Esz,Euz,Eu_0m,PARz,SWRz) Ed(1,jj,ji,:) = Ed_0m(:,jj,ji) Es(1,jj,ji,:) = Es_0m(:,jj,ji) @@ -142,6 +143,11 @@ SUBROUTINE trc3streams(datestring) enddo enddo + do jk =1, bottom + SWR_RT(jk,jj,ji) = SWRz(jk) + enddo + + ENDIF IF (MODE .EQ. 2) then From 9130a524947e6ec64ee722a7697aff72f458ea35 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 9 Dec 2024 14:14:12 +0100 Subject: [PATCH 05/21] Adjustment for model depths --- GeneralCmake.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/GeneralCmake.cmake b/GeneralCmake.cmake index deefd27d..d0a04967 100644 --- a/GeneralCmake.cmake +++ b/GeneralCmake.cmake @@ -31,6 +31,7 @@ set(CMAKE_Fortran_COMPILER ${MPI_Fortran_COMPILER}) add_definitions(-Dkey_trahdfcoef1d -Dkey_trahdfbilap -Dkey_trc_smolar) add_definitions(-Dkey_trc_hdfbilap -Dkey_trc_dmp -Dkey_kef -Dkey_trc_sed ) add_definitions(-Dkey_mpp -Dkey_mpp_mpi) +add_definitions(-Dgdept1d) IF (BFMv2) add_definitions(-DBFMv2) ENDIF() From 20b1c0b674b46d2e77be4e87711c1f095009d25a Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Tue, 10 Dec 2024 16:38:26 +0100 Subject: [PATCH 06/21] Bug fix: adding #ifdef gdept1d in forcing_phys, about Eddy viscosity calculation --- src/IO/forcing_phys.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/IO/forcing_phys.f90 b/src/IO/forcing_phys.f90 index fc20ae87..622e7e2c 100644 --- a/src/IO/forcing_phys.f90 +++ b/src/IO/forcing_phys.f90 @@ -198,7 +198,11 @@ SUBROUTINE LOAD_PHYS(datestring) do jj=1,jpj do jk=1,jpk if (tmask(jk,jj,ji) .ne. 0.) then +#ifdef gdept1d + avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk)/(sigma*buf2(jj,ji)) ) **2) + DvBackground +#else avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk,jj,ji)/(sigma*buf2(jj,ji)) ) **2) + DvBackground +#endif endif enddo enddo From fac2b7b0477c9b0a58a94d2c9b45b26cbf9bc508 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 11 Dec 2024 18:20:51 +0100 Subject: [PATCH 07/21] Adjusted boundary conditions for mercury --- src/BC/gib.nml | 16 ++++++++++++++++ src/BC/nudging.f90 | 8 ++++++-- src/BC/rivers.f90 | 7 +++++-- src/BC/sponge.f90 | 9 ++++++--- src/IO/BC_mem.f90 | 2 +- 5 files changed, 34 insertions(+), 8 deletions(-) diff --git a/src/BC/gib.nml b/src/BC/gib.nml index 5fbf7563..cf73b44f 100644 --- a/src/BC/gib.nml +++ b/src/BC/gib.nml @@ -14,6 +14,14 @@ vars(6) = "O3h" vars(7) = "N6r" + var_names_idx(1) = 1 + var_names_idx(2) = 2 + var_names_idx(3) = 3 + var_names_idx(4) = 6 + var_names_idx(5) = 49 + var_names_idx(6) = 50 + var_names_idx(7) = 7 + alpha = 4.0d0 reduction_value_t = 1.0d-6 length = -7.5d0 @@ -38,6 +46,14 @@ vars(6) = "O3h" vars(7) = "N6r" + var_names_idx(1) = 1 + var_names_idx(2) = 2 + var_names_idx(3) = 3 + var_names_idx(4) = 6 + var_names_idx(5) = 49 + var_names_idx(6) = 50 + var_names_idx(7) = 7 + rst_corr(1) = 1.0d0 rst_corr(2) = 1.0d0 rst_corr(3) = 1.0d0 diff --git a/src/BC/nudging.f90 b/src/BC/nudging.f90 index 8fddad75..1dc4a802 100644 --- a/src/BC/nudging.f90 +++ b/src/BC/nudging.f90 @@ -80,6 +80,7 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) integer :: n_vars character(len=11) :: data_file ! 11 chars in order to handle names like 'bounmask.nc' character(len=20), allocatable, dimension(:) :: vars + integer(4), allocatable, dimension(:) :: var_names_idx double precision, allocatable, dimension(:) :: rst_corr integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files @@ -90,7 +91,7 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) integer :: i namelist /nudging_vars_dimension/ n_vars - namelist /nudging_core/ data_file, vars, rst_corr + namelist /nudging_core/ data_file, vars, var_names_idx, rst_corr ! pointer to bc_no_nudging association self%m_bc_no_nudging => bc_no_nudging @@ -104,6 +105,7 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) ! allocate local arrays allocate(vars(self%m_n_nudging_vars)) + allocate(var_names_idx(self%m_n_nudging_vars)) allocate(rst_corr(self%m_n_nudging_vars)) ! allocate class members @@ -126,13 +128,15 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) self%m_nudging_vars(i) = vars(i) self%m_nudging_vars_rst(i) = 're'//trim(self%m_nudging_vars(i)) call readnc_slice_float(self%m_data_file, trim(self%m_nudging_vars_rst(i)), self%m_rst(:, :, :, i), 0) - self%m_nudging_vars_idx(i) = find_index_var(self%m_nudging_vars(i)) + self%m_nudging_vars_idx(i) = var_names_idx(i) + ! find_index_var(self%m_nudging_vars(i)) self%m_rst_corr(i) = rst_corr(i) self%m_rst_tracers(:, :, :, self%m_nudging_vars_idx(i)) = self%m_rst_corr(i) * self%m_rst(:, :, :, i) enddo ! deallocation deallocate(vars) + deallocate(var_names_idx) deallocate(rst_corr) ! close file diff --git a/src/BC/rivers.f90 b/src/BC/rivers.f90 index 53808e9b..4052a928 100644 --- a/src/BC/rivers.f90 +++ b/src/BC/rivers.f90 @@ -140,10 +140,11 @@ subroutine init_members(self, bc_name, namelist_file) integer :: n_vars character(len=20), allocatable, dimension(:) :: vars + integer(4), allocatable, dimension(:) :: var_names_idx integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files integer :: i namelist /vars_dimension/ n_vars - namelist /core/ vars + namelist /core/ vars, var_names_idx self%m_name = bc_name @@ -156,6 +157,7 @@ subroutine init_members(self, bc_name, namelist_file) ! allocate local arrays allocate(vars(self%m_n_vars)) + allocate(var_names_idx(self%m_n_vars)) ! allocate class members allocate(self%m_var_names(self%m_n_vars)) @@ -169,7 +171,7 @@ subroutine init_members(self, bc_name, namelist_file) do i = 1, self%m_n_vars self%m_var_names(i) = vars(i) self%m_var_names_data(i) = "riv"//'_'//trim(self%m_var_names(i)) - self%m_var_names_idx(i) = find_index_var(self%m_var_names(i)) + self%m_var_names_idx(i) = var_names_idx(i) !find_index_var(self%m_var_names(i)) if (lwp) write(*,*) 'RIV ', self%m_var_names(i), self%m_var_names_idx(i) enddo @@ -184,6 +186,7 @@ subroutine init_members(self, bc_name, namelist_file) ! deallocation deallocate(vars) + deallocate(var_names_idx) ! close file close(unit=file_unit) diff --git a/src/BC/sponge.f90 b/src/BC/sponge.f90 index a912bf5c..71d6c1bf 100644 --- a/src/BC/sponge.f90 +++ b/src/BC/sponge.f90 @@ -150,13 +150,14 @@ subroutine init_members(self, bc_name, namelist_file) integer :: n_vars character(len=20), allocatable, dimension(:) :: vars + integer(4), allocatable, dimension(:) :: var_names_idx double precision :: alpha double precision :: reduction_value_t double precision :: length integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files integer :: i namelist /vars_dimension/ n_vars - namelist /core/ vars, alpha, reduction_value_t, length + namelist /core/ vars, var_names_idx, alpha, reduction_value_t, length self%m_name = bc_name @@ -169,6 +170,7 @@ subroutine init_members(self, bc_name, namelist_file) ! allocate local arrays allocate(vars(self%m_n_vars)) + allocate(var_names_idx(self%m_n_vars)) ! allocate class members allocate(self%m_var_names(self%m_n_vars)) @@ -181,8 +183,8 @@ subroutine init_members(self, bc_name, namelist_file) do i = 1, self%m_n_vars self%m_var_names(i) = vars(i) - self%m_var_names_data(i) = trim(self%m_var_names(i)) - self%m_var_names_idx(i) = find_index_var(self%m_var_names(i)) + self%m_var_names_data(i) = self%m_name//'_'//trim(self%m_var_names(i)) + self%m_var_names_idx(i) = var_names_idx(i) enddo ! call delegated constructor - related procedures @@ -200,6 +202,7 @@ subroutine init_members(self, bc_name, namelist_file) ! deallocation deallocate(vars) + deallocate(var_names_idx) ! close file close(unit=file_unit) diff --git a/src/IO/BC_mem.f90 b/src/IO/BC_mem.f90 index 5294e5c3..4ad3e796 100644 --- a/src/IO/BC_mem.f90 +++ b/src/IO/BC_mem.f90 @@ -112,7 +112,7 @@ SUBROUTINE alloc_DTATRC ! N3n jn = 3 nitrates (climatological) ! N5s jn = 4 Silicates (seasonal) - jn_gib = 7 + jn_gib = 6 ! jn_riv = 6 jn_atm = 4 From f9f6f98ac5183eadc7a87b1bce92687458033dfc Mon Sep 17 00:00:00 2001 From: ginevrar Date: Thu, 12 Dec 2024 11:33:39 +0100 Subject: [PATCH 08/21] Restored original indexing in trcbio.f90 --- src/BIO/trcbio.f90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/BIO/trcbio.f90 b/src/BIO/trcbio.f90 index f6c452df..e1b48ef6 100644 --- a/src/BIO/trcbio.f90 +++ b/src/BIO/trcbio.f90 @@ -115,7 +115,7 @@ SUBROUTINE trcbio correct_fact= 0.0D0 endif - er(jk,11) = correct_fact * ( gdept(jpk)-gdept(jk) ) /gdept(jpk) + er(jk,16) = correct_fact * ( gdept(jpk)-gdept(jk) ) /gdept(jpk) enddo #endif @@ -144,12 +144,11 @@ SUBROUTINE trcbio er(1:bottom,8) = PAR(1:bottom,jj,ji,3) ! PAR for pico phytoplankton er(1:bottom,9) = PAR(1:bottom,jj,ji,4) ! PAR for dinoflagellates er(1:bottom,10) = PAR(1:bottom,jj,ji,5) ! total PAR for CDOM - er(1:bottom,11) = SWR_RT(1:bottom,jj,ji) - er(1 ,12) = DAY_LENGTH(jj,ji) ! fotoperiod expressed in hours - er(1:bottom,13) = e3t(1:bottom,jj,ji) ! depth in meters of the given cell - er(1 ,14) = vatm(jj,ji) ! wind speed (m/s) - er(1:bottom,15) = ogstm_PH(1:bottom,jj,ji) ! 8.1 - er(1 ,16) = RMU(jj,ji) ! avg. cosine direct + er(1 ,11) = DAY_LENGTH(jj,ji) ! fotoperiod expressed in hours + er(1:bottom,12) = e3t(1:bottom,jj,ji) ! depth in meters of the given cell + er(1 ,13) = vatm(jj,ji) ! wind speed (m/s) + er(1:bottom,14) = ogstm_PH(1:bottom,jj,ji) ! 8.1 + er(1 ,15) = RMU(jj,ji) ! avg. cosine direct #ifndef gdept1d do jk=1,bottom @@ -162,11 +161,12 @@ SUBROUTINE trcbio correct_fact= 0.0D0 endif - er(jk,17) = correct_fact * ( gdept(jpk,jj,ji)-gdept(jk,jj,ji) ) /gdept(jpk,jj,ji) + er(jk,16) = correct_fact * ( gdept(jpk,jj,ji)-gdept(jk,jj,ji) ) /gdept(jpk,jj,ji) enddo #endif - er(1 ,18) = atm_Hg0(jj,ji) ! ATM Hg0 + er(1:bottom,17) = SWR_RT(1:bottom,jj,ji) + er(1 ,18) = atm_Hg0(jj,ji) call BFM1D_Input_EcologyDynamics(bottom,a,jtrmax,er) call BFM1D_reset() From 25123e83b5e9bf3df7df96a3d31c92e2f0a4fdc5 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:11:26 +0100 Subject: [PATCH 09/21] Removed var_names_idx for rivers --- src/BC/rivers.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/BC/rivers.f90 b/src/BC/rivers.f90 index 4052a928..48b38048 100644 --- a/src/BC/rivers.f90 +++ b/src/BC/rivers.f90 @@ -140,11 +140,11 @@ subroutine init_members(self, bc_name, namelist_file) integer :: n_vars character(len=20), allocatable, dimension(:) :: vars - integer(4), allocatable, dimension(:) :: var_names_idx + ! integer(4), allocatable, dimension(:) :: var_names_idx integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files integer :: i namelist /vars_dimension/ n_vars - namelist /core/ vars, var_names_idx + namelist /core/ vars self%m_name = bc_name @@ -157,7 +157,7 @@ subroutine init_members(self, bc_name, namelist_file) ! allocate local arrays allocate(vars(self%m_n_vars)) - allocate(var_names_idx(self%m_n_vars)) + ! allocate(var_names_idx(self%m_n_vars)) ! allocate class members allocate(self%m_var_names(self%m_n_vars)) @@ -171,7 +171,8 @@ subroutine init_members(self, bc_name, namelist_file) do i = 1, self%m_n_vars self%m_var_names(i) = vars(i) self%m_var_names_data(i) = "riv"//'_'//trim(self%m_var_names(i)) - self%m_var_names_idx(i) = var_names_idx(i) !find_index_var(self%m_var_names(i)) + self%m_var_names_idx(i) = find_index_var(self%m_var_names(i)) + ! self%m_var_names_idx(i) = var_names_idx(i) !find_index_var(self%m_var_names(i)) if (lwp) write(*,*) 'RIV ', self%m_var_names(i), self%m_var_names_idx(i) enddo @@ -186,7 +187,7 @@ subroutine init_members(self, bc_name, namelist_file) ! deallocation deallocate(vars) - deallocate(var_names_idx) + ! deallocate(var_names_idx) ! close file close(unit=file_unit) From c95a90399b3b7195a5920f8706c6ee34eb51c9d1 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:16:09 +0100 Subject: [PATCH 10/21] Restored original CHAR lenghts for var names and Removed TRIMs --- src/IO/domrea.f90 | 6 +++--- src/IO/trcrst.f90 | 50 +++++++++++++++++++++++------------------------ src/IO/trcwri.f90 | 6 +++--- 3 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/IO/domrea.f90 b/src/IO/domrea.f90 index e241dfd1..6e6eda32 100644 --- a/src/IO/domrea.f90 +++ b/src/IO/domrea.f90 @@ -49,8 +49,8 @@ SUBROUTINE domrea character(len=10) bfmmask_file CHARACTER(LEN=50) filename - CHARACTER(LEN=20), DIMENSION(7) :: var_nc - CHARACTER(LEN=22) nomevar01 + CHARACTER(LEN=3), DIMENSION(7) :: var_nc + CHARACTER(LEN=5) nomevar01 LOGICAL B ! ------------------- @@ -209,7 +209,7 @@ SUBROUTINE domrea IF (NWATERPOINTS.GT.0) THEN do jn=1,jn_gib - nomevar01='re'//TRIM(var_nc(jn)) + nomevar01='re'//var_nc(jn) call readnc_slice_float('bounmask.nc',nomevar01,resto(:,:,:,jn),0) enddo diff --git a/src/IO/trcrst.f90 b/src/IO/trcrst.f90 index 754c7b2a..4e3baa29 100644 --- a/src/IO/trcrst.f90 +++ b/src/IO/trcrst.f90 @@ -22,7 +22,7 @@ SUBROUTINE trcrst ! local declarations ! ================== INTEGER jn, jn_high - CHARACTER(LEN=54) filename + CHARACTER(LEN=37) filename CHARACTER(LEN=100) bkpname logical existFile logical bkp1hasbeenread,bkp2hasbeenread @@ -49,7 +49,7 @@ SUBROUTINE trcrst filename = 'RESTARTS/RST.'//DateStart//'.'//trim(ctrcnm(jn))// & '.nc' - CALL readnc_slice_double(TRIM(filename), 'TRN'//trim(ctrcnm(jn)), trn(:,:,:,jn) ) + CALL readnc_slice_double(filename, 'TRN'//trim(ctrcnm(jn)), trn(:,:,:,jn) ) @@ -67,10 +67,10 @@ SUBROUTINE trcrst if (existFilebkp) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(TRIM(bkpname),TRIM(ctrcnm(jn)), traIO(:,:,:,jn) ) + CALL readnc_slice_double(bkpname,ctrcnm(jn), traIO(:,:,:,jn) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -82,14 +82,14 @@ SUBROUTINE trcrst IF (ctr_hf(jn).eq.1) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(ctrcnm(jn))//'.nc.bkp' - INQUIRE(FILE=TRIM(bkpname), EXIST=existFilebkp) + INQUIRE(FILE=bkpname, EXIST=existFilebkp) if (existFilebkp) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(TRIM(bkpname),TRIM(ctrcnm(jn)), traIO_HIGH(:,:,:,jn_high) ) + CALL readnc_slice_double(bkpname,ctrcnm(jn), traIO_HIGH(:,:,:,jn_high) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else @@ -106,14 +106,14 @@ SUBROUTINE trcrst DO jn=1, jptra_dia bkpname = 'AVE_FREQ_2/ave.'//DateStart//'.'//trim(dianm(jn))//'.nc.bkp' - INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) + INQUIRE(FILE=bkpname, EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(TRIM(bkpname),TRIM(dianm(jn)), tra_DIA_IO(jn,:,:,:) ) + CALL readnc_slice_double(bkpname,dianm(jn), tra_DIA_IO(jn,:,:,:) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -123,13 +123,13 @@ SUBROUTINE trcrst IF ((diahf(jn).eq.1).and.(diaWR(jn).eq.1)) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(dianm(jn))//'.nc.bkp' - INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) + INQUIRE(FILE=bkpname, EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double(TRIM(bkpname),trim(dianm(jn)), tra_DIA_IO_HIGH(jn_high,:,:,:) ) + CALL readnc_slice_double(bkpname,trim(dianm(jn)), tra_DIA_IO_HIGH(jn_high,:,:,:) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else @@ -146,14 +146,14 @@ SUBROUTINE trcrst DO jn=1, jptra_dia_2d bkpname = 'AVE_FREQ_2/ave.'//DateStart//'.'//trim(dianm_2d(jn))//'.nc.bkp' - INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) + INQUIRE(FILE=bkpname, EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double_2d(TRIM(bkpname),dianm_2d(jn), tra_DIA_2d_IO(jn,:,:) ) + CALL readnc_slice_double_2d(bkpname,dianm_2d(jn), tra_DIA_2d_IO(jn,:,:) ) if (.not.bkp2hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_2) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_2) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_2) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_2) bkp2hasbeenread=.true. endif else @@ -163,13 +163,13 @@ SUBROUTINE trcrst IF ((diahf_2d(jn).eq.1).and.(diaWR_2d(jn).eq.1)) THEN jn_high = jn_high + 1 bkpname= 'AVE_FREQ_1/ave.'//DateStart//'.'//trim(dianm_2d(jn))//'.nc.bkp' - INQUIRE(FILE=TRIM(bkpname), EXIST=existFile) + INQUIRE(FILE=bkpname, EXIST=existFile) if (existFile) then if (lwp) write(*,*) 'reading ', bkpname - CALL readnc_slice_double_2d(TRIM(bkpname),trim(dianm_2d(jn)), tra_DIA_2d_IO_HIGH(jn_high,:,:) ) + CALL readnc_slice_double_2d(bkpname,trim(dianm_2d(jn)), tra_DIA_2d_IO_HIGH(jn_high,:,:) ) if (.not.bkp1hasbeenread) then - call readnc_scalar_double(TRIM(bkpname),'elapsed_time',elapsed_time_1) - call get_att_char(TRIM(bkpname),'DateStart' , BKPdatefrom_1) + call readnc_scalar_double(bkpname,'elapsed_time',elapsed_time_1) + call get_att_char(bkpname,'DateStart' , BKPdatefrom_1) bkp1hasbeenread=.true. endif else diff --git a/src/IO/trcwri.f90 b/src/IO/trcwri.f90 index 055977d1..c5ff89e5 100644 --- a/src/IO/trcwri.f90 +++ b/src/IO/trcwri.f90 @@ -26,9 +26,9 @@ SUBROUTINE trcwri(datestring) double precision julian - CHARACTER(LEN=54) filename + CHARACTER(LEN=37) filename - CHARACTER(LEN=20) varname + CHARACTER(LEN=3) varname INTEGER idrank, ierr, istart, jstart, iPe, iPd, jPe, jPd, status(MPI_STATUS_SIZE) INTEGER irange, jrange @@ -121,7 +121,7 @@ SUBROUTINE trcwri(datestring) filename = 'RESTARTS/RST.'//datestring//'.'//trim(var_to_store)//'.nc' - CALL write_restart(TRIM(filename),TRIM(var_to_store),julian, deflate_rst, deflate_level_rst) + CALL write_restart(filename,var_to_store,julian, deflate_rst, deflate_level_rst) END IF END IF From 572251b576bb3da2fbdd1148c65fd053fd43b731 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:20:45 +0100 Subject: [PATCH 11/21] Fixed var_names_idx in sponge and nudging.f90 --- src/BC/nudging.f90 | 8 ++------ src/BC/sponge.f90 | 9 +++------ 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/src/BC/nudging.f90 b/src/BC/nudging.f90 index 1dc4a802..8fddad75 100644 --- a/src/BC/nudging.f90 +++ b/src/BC/nudging.f90 @@ -80,7 +80,6 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) integer :: n_vars character(len=11) :: data_file ! 11 chars in order to handle names like 'bounmask.nc' character(len=20), allocatable, dimension(:) :: vars - integer(4), allocatable, dimension(:) :: var_names_idx double precision, allocatable, dimension(:) :: rst_corr integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files @@ -91,7 +90,7 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) integer :: i namelist /nudging_vars_dimension/ n_vars - namelist /nudging_core/ data_file, vars, var_names_idx, rst_corr + namelist /nudging_core/ data_file, vars, rst_corr ! pointer to bc_no_nudging association self%m_bc_no_nudging => bc_no_nudging @@ -105,7 +104,6 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) ! allocate local arrays allocate(vars(self%m_n_nudging_vars)) - allocate(var_names_idx(self%m_n_nudging_vars)) allocate(rst_corr(self%m_n_nudging_vars)) ! allocate class members @@ -128,15 +126,13 @@ subroutine init_members(self, bc_no_nudging, namelist_file, n_tracers) self%m_nudging_vars(i) = vars(i) self%m_nudging_vars_rst(i) = 're'//trim(self%m_nudging_vars(i)) call readnc_slice_float(self%m_data_file, trim(self%m_nudging_vars_rst(i)), self%m_rst(:, :, :, i), 0) - self%m_nudging_vars_idx(i) = var_names_idx(i) - ! find_index_var(self%m_nudging_vars(i)) + self%m_nudging_vars_idx(i) = find_index_var(self%m_nudging_vars(i)) self%m_rst_corr(i) = rst_corr(i) self%m_rst_tracers(:, :, :, self%m_nudging_vars_idx(i)) = self%m_rst_corr(i) * self%m_rst(:, :, :, i) enddo ! deallocation deallocate(vars) - deallocate(var_names_idx) deallocate(rst_corr) ! close file diff --git a/src/BC/sponge.f90 b/src/BC/sponge.f90 index 71d6c1bf..a912bf5c 100644 --- a/src/BC/sponge.f90 +++ b/src/BC/sponge.f90 @@ -150,14 +150,13 @@ subroutine init_members(self, bc_name, namelist_file) integer :: n_vars character(len=20), allocatable, dimension(:) :: vars - integer(4), allocatable, dimension(:) :: var_names_idx double precision :: alpha double precision :: reduction_value_t double precision :: length integer, parameter :: file_unit = 101 ! 100 for data files, 101 for boundary namelist files integer :: i namelist /vars_dimension/ n_vars - namelist /core/ vars, var_names_idx, alpha, reduction_value_t, length + namelist /core/ vars, alpha, reduction_value_t, length self%m_name = bc_name @@ -170,7 +169,6 @@ subroutine init_members(self, bc_name, namelist_file) ! allocate local arrays allocate(vars(self%m_n_vars)) - allocate(var_names_idx(self%m_n_vars)) ! allocate class members allocate(self%m_var_names(self%m_n_vars)) @@ -183,8 +181,8 @@ subroutine init_members(self, bc_name, namelist_file) do i = 1, self%m_n_vars self%m_var_names(i) = vars(i) - self%m_var_names_data(i) = self%m_name//'_'//trim(self%m_var_names(i)) - self%m_var_names_idx(i) = var_names_idx(i) + self%m_var_names_data(i) = trim(self%m_var_names(i)) + self%m_var_names_idx(i) = find_index_var(self%m_var_names(i)) enddo ! call delegated constructor - related procedures @@ -202,7 +200,6 @@ subroutine init_members(self, bc_name, namelist_file) ! deallocation deallocate(vars) - deallocate(var_names_idx) ! close file close(unit=file_unit) From f7d4aec80304b852de2d9e05874548b7bb902efa Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:32:50 +0100 Subject: [PATCH 12/21] Bug fix, added reset fro bc_atm_Hg0_TotTime --- src/General/Timers.f90 | 1 + src/IO/forcing_phys.f90 | 4 ---- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/src/General/Timers.f90 b/src/General/Timers.f90 index 9ef7f936..af522d2f 100644 --- a/src/General/Timers.f90 +++ b/src/General/Timers.f90 @@ -44,6 +44,7 @@ SUBROUTINE reset_Timers() forcing_phys_TotTime = 0.0 forcing_kext_TotTime = 0.0 bc_Co2_TotTime = 0.0 + bc_atm_Hg0_TotTime = 0.0 bc_atm_TotTime = 0.0 bc_tin_TotTime = 0.0 bc_gib_TotTime = 0.0 diff --git a/src/IO/forcing_phys.f90 b/src/IO/forcing_phys.f90 index 622e7e2c..fc20ae87 100644 --- a/src/IO/forcing_phys.f90 +++ b/src/IO/forcing_phys.f90 @@ -198,11 +198,7 @@ SUBROUTINE LOAD_PHYS(datestring) do jj=1,jpj do jk=1,jpk if (tmask(jk,jj,ji) .ne. 0.) then -#ifdef gdept1d - avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk)/(sigma*buf2(jj,ji)) ) **2) + DvBackground -#else avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk,jj,ji)/(sigma*buf2(jj,ji)) ) **2) + DvBackground -#endif endif enddo enddo From 848fc29a9497b0e3b7a42054b7e9118da5fe4b60 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:33:42 +0100 Subject: [PATCH 13/21] Bug fix for er(jk,11) --- src/BIO/trcbio.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/BIO/trcbio.f90 b/src/BIO/trcbio.f90 index e1b48ef6..26ae263e 100644 --- a/src/BIO/trcbio.f90 +++ b/src/BIO/trcbio.f90 @@ -115,7 +115,7 @@ SUBROUTINE trcbio correct_fact= 0.0D0 endif - er(jk,16) = correct_fact * ( gdept(jpk)-gdept(jk) ) /gdept(jpk) + er(jk,11) = correct_fact * ( gdept(jpk)-gdept(jk) ) /gdept(jpk) enddo #endif @@ -167,7 +167,8 @@ SUBROUTINE trcbio er(1:bottom,17) = SWR_RT(1:bottom,jj,ji) er(1 ,18) = atm_Hg0(jj,ji) - call BFM1D_Input_EcologyDynamics(bottom,a,jtrmax,er) + + call BFM1D_Input_EcologyDynamics(bottom,a,jtrmax,er) call BFM1D_reset() From 35dc3e9b27370f3b086241c3eb717c174dc3a9b9 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 12:56:12 +0100 Subject: [PATCH 14/21] Removed var_names_idx in gib.nml --- src/BC/gib.nml | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/BC/gib.nml b/src/BC/gib.nml index cf73b44f..5fbf7563 100644 --- a/src/BC/gib.nml +++ b/src/BC/gib.nml @@ -14,14 +14,6 @@ vars(6) = "O3h" vars(7) = "N6r" - var_names_idx(1) = 1 - var_names_idx(2) = 2 - var_names_idx(3) = 3 - var_names_idx(4) = 6 - var_names_idx(5) = 49 - var_names_idx(6) = 50 - var_names_idx(7) = 7 - alpha = 4.0d0 reduction_value_t = 1.0d-6 length = -7.5d0 @@ -46,14 +38,6 @@ vars(6) = "O3h" vars(7) = "N6r" - var_names_idx(1) = 1 - var_names_idx(2) = 2 - var_names_idx(3) = 3 - var_names_idx(4) = 6 - var_names_idx(5) = 49 - var_names_idx(6) = 50 - var_names_idx(7) = 7 - rst_corr(1) = 1.0d0 rst_corr(2) = 1.0d0 rst_corr(3) = 1.0d0 From f91b3e8483dc37e7f95d07b825b10d3dbfb4a3c9 Mon Sep 17 00:00:00 2001 From: ginevrar Date: Wed, 22 Jan 2025 14:24:03 +0100 Subject: [PATCH 15/21] Bug fix following commit 20b1c0b on neccton branch --- src/IO/forcing_phys.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/IO/forcing_phys.f90 b/src/IO/forcing_phys.f90 index fc20ae87..18f2f01c 100644 --- a/src/IO/forcing_phys.f90 +++ b/src/IO/forcing_phys.f90 @@ -198,7 +198,11 @@ SUBROUTINE LOAD_PHYS(datestring) do jj=1,jpj do jk=1,jpk if (tmask(jk,jj,ji) .ne. 0.) then - avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk,jj,ji)/(sigma*buf2(jj,ji)) ) **2) + DvBackground +#ifdef gdept1d + avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk)/(sigma*buf2(jj,ji)) ) **2) + DvBackground +#else + avtdta(jk,jj,ji,2) = DvMld * exp(-0.5* ( gdept(jk,jj,ji)/(sigma*buf2(jj,ji)) ) **2) + DvBackgroun +#endif endif enddo enddo From 32466244b7b7ec199841d7d18f796906c67f9c9d Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Thu, 30 Jan 2025 17:54:15 +0100 Subject: [PATCH 16/21] Increasing RST filename lengths --- src/IO/trcrst.f90 | 2 +- src/IO/trcwri.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/IO/trcrst.f90 b/src/IO/trcrst.f90 index 4e3baa29..6a2eced9 100644 --- a/src/IO/trcrst.f90 +++ b/src/IO/trcrst.f90 @@ -22,7 +22,7 @@ SUBROUTINE trcrst ! local declarations ! ================== INTEGER jn, jn_high - CHARACTER(LEN=37) filename + CHARACTER(LEN=54) filename CHARACTER(LEN=100) bkpname logical existFile logical bkp1hasbeenread,bkp2hasbeenread diff --git a/src/IO/trcwri.f90 b/src/IO/trcwri.f90 index c5ff89e5..22a4fe08 100644 --- a/src/IO/trcwri.f90 +++ b/src/IO/trcwri.f90 @@ -26,7 +26,7 @@ SUBROUTINE trcwri(datestring) double precision julian - CHARACTER(LEN=37) filename + CHARACTER(LEN=54) filename CHARACTER(LEN=3) varname @@ -38,7 +38,7 @@ SUBROUTINE trcwri(datestring) CHARACTER(LEN=20) var_to_store INTEGER :: COUNTER_VAR_TRCWRI, n_dumping_cycles, jv, ivar, writing_rank, ind_col - filename = 'RST.20111231-15:30:00.N1p.nc' + filename = 'RESTARTS/RST.20111231-15:30:00.N1p.nc' julian=datestring2sec(datestring) trcwriparttime = MPI_WTIME() ! cronometer-start From 64220e81007c9f5efe985c17d3487bb60e415d1e Mon Sep 17 00:00:00 2001 From: ginevrar Date: Mon, 3 Feb 2025 12:26:20 +0100 Subject: [PATCH 17/21] Added two remPOC diagnostic to BFMtab.xml - NECCTON dependencies --- bfmv5/BFMtab.xml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bfmv5/BFMtab.xml b/bfmv5/BFMtab.xml index 34a57d87..f8cb4b59 100644 --- a/bfmv5/BFMtab.xml +++ b/bfmv5/BFMtab.xml @@ -179,7 +179,8 @@ - + + From 16af5c619ce11ba6a474e6b5ec907d6c18d50bbb Mon Sep 17 00:00:00 2001 From: ginevrar Date: Thu, 6 Feb 2025 18:35:30 +0100 Subject: [PATCH 18/21] Modified to read 4 string Vars in RIV boundaries --- src/BC/rivers.f90 | 4 ++-- src/IO/BC_mem.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/BC/rivers.f90 b/src/BC/rivers.f90 index 48b38048..1a497cad 100644 --- a/src/BC/rivers.f90 +++ b/src/BC/rivers.f90 @@ -19,8 +19,8 @@ module rivers_mod ! TO DO: review names character(len=3) :: m_name ! ex: 'riv' integer :: m_n_vars ! BC_mem.f90:95 - character(len=20), allocatable, dimension(:) :: m_var_names - character(len=20), allocatable, dimension(:) :: m_var_names_data ! bc_tin.f90:116 + character(len=26), allocatable, dimension(:) :: m_var_names + character(len=26), allocatable, dimension(:) :: m_var_names_data ! bc_tin.f90:116 integer(4), allocatable, dimension(:) :: m_var_names_idx ! tra_matrix_riv double precision, allocatable, dimension(:, :) :: m_buffer ! replaces m_aux, now it is a 2D matrix integer(4) :: m_size diff --git a/src/IO/BC_mem.f90 b/src/IO/BC_mem.f90 index 4ad3e796..982ff62d 100644 --- a/src/IO/BC_mem.f90 +++ b/src/IO/BC_mem.f90 @@ -67,7 +67,7 @@ SUBROUTINE alloc_DTATRC INTEGER :: err double precision :: aux_mem - CHARACTER(LEN=50) atmfile, rivfile, gibfile + CHARACTER(LEN=54) atmfile, rivfile, gibfile #ifdef Mem_Monitor aux_mem = get_mem(err) From 50f1b4b9e93064263474b141582d09898292678e Mon Sep 17 00:00:00 2001 From: ginevrar Date: Thu, 6 Feb 2025 18:53:05 +0100 Subject: [PATCH 19/21] Tried to read 4 string state vars in GIB - Not Working --- src/BC/bc_set.f90 | 2 +- src/BC/sponge.f90 | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/BC/bc_set.f90 b/src/BC/bc_set.f90 index e28c973c..96e4d4a1 100644 --- a/src/BC/bc_set.f90 +++ b/src/BC/bc_set.f90 @@ -47,7 +47,7 @@ subroutine bc_set_default(bcs,bcs_namelist) type(bc_set), intent(inout) :: bcs character(len=14), intent(in) :: bcs_namelist integer, parameter :: file_unit = 102 ! 100 for data files, 101 for boundary namelist files, 102 for global namelist - character(len=47) :: bc_string + character(len=54) :: bc_string integer :: i ! open file diff --git a/src/BC/sponge.f90 b/src/BC/sponge.f90 index a912bf5c..28ad9de8 100644 --- a/src/BC/sponge.f90 +++ b/src/BC/sponge.f90 @@ -18,8 +18,8 @@ module sponge_mod ! TO DO: review names character(len=3) :: m_name ! ex: 'gib' integer :: m_n_vars ! BC_mem.f90:94 - character(len=20), allocatable, dimension(:) :: m_var_names ! domrea.f90:161-167 - character(len=20), allocatable, dimension(:) :: m_var_names_data ! bc_gib.f90:113 + character(len=29), allocatable, dimension(:) :: m_var_names ! domrea.f90:161-167 + character(len=29), allocatable, dimension(:) :: m_var_names_data ! bc_gib.f90:113 integer(4), allocatable, dimension(:) :: m_var_names_idx ! tra_matrix_gib double precision, allocatable, dimension(:, :, :) :: m_buffer ! replaces m_aux, now it is a 3D matrix integer(4) :: m_size ! BC_mem.f90:21 @@ -183,6 +183,8 @@ subroutine init_members(self, bc_name, namelist_file) self%m_var_names(i) = vars(i) self%m_var_names_data(i) = trim(self%m_var_names(i)) self%m_var_names_idx(i) = find_index_var(self%m_var_names(i)) + write(*,*) 'self%m_var_names_data(i)', self%m_var_names_data(i) + write(*,*) 'self%m_var_names_idx)', self%m_var_names_idx(i) enddo ! call delegated constructor - related procedures From 0beb2f9b8d32f374203aeed98137026fdeba4e7b Mon Sep 17 00:00:00 2001 From: ginevrar Date: Tue, 25 Feb 2025 15:19:46 +0100 Subject: [PATCH 20/21] Changed required variable name into hg0atm for the input file atm_Hg0_yyyy0107-00:00:00.nc --- src/IO/bc_atm_Hg0.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IO/bc_atm_Hg0.f90 b/src/IO/bc_atm_Hg0.f90 index a1aa47df..c7f6ca7d 100644 --- a/src/IO/bc_atm_Hg0.f90 +++ b/src/IO/bc_atm_Hg0.f90 @@ -109,7 +109,7 @@ SUBROUTINE LOAD_atm_Hg0(datestring) nomefile(1:31)='BC/atm_Hg0_'//datestring//'.nc' if(lwp) write(*,'(A,I4,A,A)') "LOAD_atm_Hg0 --> I am ", myrank, " starting reading forcing fields from ", nomefile(1:31) - call readnc_slice_float_2d(nomefile,'hgconc',buf2,0) + call readnc_slice_float_2d(nomefile,'hg0atm',buf2,0) Hg0_IO(:,:,2) = buf2*tmask(1,:,:) From d5ef88fca5576f994b496854c8d5eb452f9e22ac Mon Sep 17 00:00:00 2001 From: ginRosati Date: Wed, 26 Feb 2025 11:24:13 +0100 Subject: [PATCH 21/21] Added README.md for mercury (Hg) model --- README.md | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 README.md diff --git a/README.md b/README.md new file mode 100644 index 00000000..01e2dbd1 --- /dev/null +++ b/README.md @@ -0,0 +1,21 @@ +## The OGSTM-BFM-Hg model + +OGSTM-BFM-Hg is a transport-reaction model coupling mercury biogeochemistry in the ocean with the plankton-nutrients-detritus dynamics (BFM model) and hydrodynamic transport (OGSTM model). The model simulates the main marine Hg species (HgII, Hg0, MMHg, and DMHg) as well as other key processes of the Hg cycle such as bioaccumulation in phyto- and zooplankton, partitioning to detritus, sinking to the seabed, and outgassing to the atmosphere. The tool is written in Fortran 90 and presented in Rosati et al., (2022). It was developed under the Italian PRIN project ICCC (Impact of Climate Change on the biogeochemistry of Contaminants in the Mediterranean Sea) and improved as part of the NECCTON project. + +**References** + +*Rosati, G., Canu, D., Lazzari, P., & Solidoro, C. (2022). Assessing the spatial and temporal variability of MeHg biogeochemistry and bioaccumulation in the Mediterranean Sea with a coupled 3D model. Biogeosciences, 19(February), 3663–3682. https://doi.org/doi.org/10.5194/bg-2022-14* + + +**Installation** + +The shell script *downloader_ogstm_bfm.sh* automatically downloads the OGSTM and BFM codes from their repositories and switches to the branch containing the mercury dynamics (branch: neccton_WP8). Once the model directories have been downloaded, the *builder_ogstm_bfm.sh* can be executed to compile the model. These file are available at: https://github.com/inogs/ModelBuild/tree/neccton_WP8. +For the procedure to be successful, a request to access the BFM repository must be made in advance to the BFM system team via this link https://docs.google.com/forms/d/e/1FAIpQLScI7N8AcvFxBeCD-EXwMXkQhgMwjhOLz3MYX8Kb47oPCXRv6w/viewform + +**Required data** + +The model namelists, generated during compilation, must be copied in the run directory ‘~/wrkdir/MODEL’: + +cp ${CODEDIR}/ogstm/ready_for_model_namelists/* . + +The run directory must include a FORCINGS directory containing the physical fields for running the model, and RESTARTS and BC directories containing files for initial and boundary conditions of biogeochemical and mercury variables. A file of the domain, called “meshmask.nc”, is required to run the model. All files used are in the netCDF format.