Neccton WP8

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()

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.

bfmv5/BFMtab.xml

@@ -2,6 +2,10 @@
 <root>
 <ModelVars>
 <var name="O2o" maxvalue="1000" highfreq="true" dump="false" ></var>
+<var name="Hg0" maxvalue="50." highfreq="true" dump="false" ></var>
+<var name="Hg2" maxvalue="50." highfreq="true" dump="false" ></var>
+<var name="MMHg" maxvalue="50." highfreq="true" dump="false" ></var>
+<var name="DMHg" maxvalue="50." highfreq="true" dump="false" ></var>
 <var name="N1p" maxvalue= "5."  highfreq="true" dump="false" ></var>
 <var name="N3n" maxvalue="500"  highfreq="true" dump="false" ></var>
 <var name="N4n" maxvalue="20"   highfreq="true" dump="false" ></var>
@@ -16,33 +20,42 @@
 <var name="P1p" maxvalue="0.2"  highfreq="false" dump="false" ></var>
 <var name="P1l" maxvalue="2"    highfreq="true" dump="false" ></var>
 <var name="P1s" maxvalue="2.5"  highfreq="false" dump="false" ></var>
+<var name="P1h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="P2c" maxvalue="125"  highfreq="true" dump="false" ></var>
 <var name="P2n" maxvalue="4"    highfreq="false" dump="false" ></var>
 <var name="P2p" maxvalue="0.5"  highfreq="false" dump="false" ></var>
 <var name="P2l" maxvalue="1.5"  highfreq="true" dump="false" ></var>
+<var name="P2h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="P3c" maxvalue="200"  highfreq="true" dump="false" ></var>
 <var name="P3n" maxvalue="6"    highfreq="false" dump="false" ></var>
 <var name="P3p" maxvalue="0.5"  highfreq="false" dump="false" ></var>
 <var name="P3l" maxvalue="4"    highfreq="true" dump="false" ></var>
+<var name="P3h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="P4c" maxvalue="100"  highfreq="true" dump="false" ></var>
 <var name="P4n" maxvalue="0.6"  highfreq="false" dump="false" ></var>
 <var name="P4p" maxvalue="0.1"  highfreq="false" dump="false" ></var>
 <var name="P4l" maxvalue="1.5"  highfreq="true" dump="false" ></var>
+<var name="P4h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="Z3c" maxvalue="18"   highfreq="false" dump="false" ></var>
 <var name="Z3n" maxvalue="0.3"  highfreq="false" dump="false" ></var>
 <var name="Z3p" maxvalue="0.1"  highfreq="false" dump="false" ></var>
+<var name="Z3h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="Z4c" maxvalue="100"  highfreq="false" dump="false" ></var>
 <var name="Z4n" maxvalue="1.5"  highfreq="false" dump="false" ></var>
 <var name="Z4p" maxvalue="0.3"  highfreq="false" dump="false" ></var>
+<var name="Z4h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="Z5c" maxvalue="80"   highfreq="false" dump="false" ></var>
 <var name="Z5n" maxvalue="3"    highfreq="false" dump="false" ></var>
 <var name="Z5p" maxvalue="0.2"  highfreq="false" dump="false" ></var>
+<var name="Z5h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="Z6c" maxvalue="100"  highfreq="false" dump="false" ></var>
 <var name="Z6n" maxvalue="4"    highfreq="false" dump="false" ></var>
 <var name="Z6p" maxvalue="0.3"  highfreq="false" dump="false" ></var>
+<var name="Z6h" maxvalue="2000" highfreq="false" dump="false" ></var>
 <var name="R1c" maxvalue="60"   highfreq="false" dump="false" ></var>
 <var name="R1n" maxvalue="3"    highfreq="false" dump="false" ></var>
 <var name="R1p" maxvalue="0.3"  highfreq="false" dump="false" ></var>
+
 <var name="R1s" maxvalue="0.1"  highfreq="false" dump="false" ></var>
 <var name="R2c" maxvalue="500"  highfreq="false" dump="false" ></var>
 <var name="R6c" maxvalue="200"  highfreq="false" dump="false" ></var>
@@ -70,6 +83,7 @@
 <var name="ESW"          highfreq="false" dump="false" ></var> 
 <var name="ERHO"         highfreq="false" dump="false" ></var> 
 <var name="EIR"          highfreq="false" dump="false" ></var> 
+<var name="SWR_RT"       highfreq="false" dump="true" ></var>
 <var name="ESS"          highfreq="false" dump="false" ></var> 
 <var name="EPR"          highfreq="false" dump="false" ></var> 
 <var name="BAC_ACT_FACT" highfreq="false" dump="false" ></var> 
@@ -78,6 +92,7 @@
 <var name="DIC"          highfreq="true" dump="true" ></var> 
 <var name="Area"         highfreq="false" dump="false" ></var> 
 <var name="CO2"          highfreq="false" dump="false" ></var> 
+<var name="skvo"         highfreq="false" dump="true" ></var> 
 <var name="pCO2"         highfreq="true" dump="true" ></var> 
 <var name="HCO3"         highfreq="false" dump="false" ></var> 
 <var name="CO3"          highfreq="false" dump="false" ></var> 
@@ -164,6 +179,17 @@
 <var name="ELiPPY_iiP4"  highfreq="false" dump="false" ></var> 
 <var name="ruPPYc"       highfreq="true"  dump="true" comment="deve essere true"></var> 
 <var name="resPPYc"      highfreq="true"  dump="true" comment="deve essere true"></var> 
+<var name="remPOC6"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="remPOC8"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skmet"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skmet2"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skdem"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skdem2"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skox"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skred"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="skphdem"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="pR2c"      highfreq="true"  dump="true" comment="deve essere true"></var>
+<var name="pR3c"      highfreq="true"  dump="true" comment="deve essere true"></var>
 <var name="netPPYc"      highfreq="false" dump="false" ></var> 
 <var name="resPBAc"      highfreq="false" dump="false" ></var> 
 <var name="resZOOc"      highfreq="false" dump="false" ></var> 
@@ -197,6 +223,15 @@
 <var name="qccPPY_iiP2"     highfreq="false" dump="false" ></var>
 <var name="qccPPY_iiP3"     highfreq="false" dump="false" ></var>
 <var name="qccPPY_iiP4"     highfreq="false" dump="false" ></var>
+<var name="qhcPPY_iiP1"     highfreq="false" dump="false" ></var>
+<var name="qhcPPY_iiP2"     highfreq="false" dump="false" ></var>
+<var name="qhcPPY_iiP3"     highfreq="false" dump="false" ></var>
+<var name="qhcPPY_iiP4"     highfreq="false" dump="false" ></var>
+<var name="qhcMEZ_iiZ3"     highfreq="false" dump="false" ></var>
+<var name="qhcMEZ_iiZ4"     highfreq="false" dump="false" ></var>
+<var name="qhcMIZ_iiZ5"     highfreq="false" dump="false" ></var>
+<var name="qhcMIZ_iiZ6"     highfreq="false" dump="false" ></var>
+<var name="qhcPBA_iiB1"     highfreq="false" dump="false" ></var>
 <var name="rcalCARc"     highfreq="false" dump="false" ></var>
 <var name="BFM1D_rsp_iiP1"     highfreq="false" dump="false" ></var>
 <var name="BFM1D_rsp_iiP2"     highfreq="false" dump="false" ></var>
@@ -236,6 +271,8 @@
 <Diagnostics_2D>
 <var name="EPCO2air"   highfreq="false" dump="false" ></var>
 <var name="CO2airflux" highfreq="true" dump="true" ></var>
+<var name="skvo" highfreq="true" dump="true" ></var>
+<var name="Hg0atm" highfreq="true" dump="true" ></var>
 <var name="Area2d"   highfreq="false" dump="false" ></var>
 <var name="ThereIsLight" highfreq="false" dump="false" ></var>
 <var name="SUNQ" highfreq="false" dump="false" ></var>

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

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
@@ -140,6 +140,7 @@ 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
@@ -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))
@@ -170,6 +172,7 @@ subroutine init_members(self, bc_name, namelist_file)
             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 +187,7 @@ subroutine init_members(self, bc_name, namelist_file)
 
         ! deallocation
         deallocate(vars)
+     !   deallocate(var_names_idx)
 
         ! close file
         close(unit=file_unit)

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

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