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() 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. diff --git a/bfmv5/BFMtab.xml b/bfmv5/BFMtab.xml index 8a7053da..f8cb4b59 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,17 @@ + + + + + + + + + + + @@ -197,6 +223,15 @@ + + + + + + + + + @@ -236,6 +271,8 @@ + + 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/rivers.f90 b/src/BC/rivers.f90 index 53808e9b..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 @@ -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) 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 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 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..26ae263e 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 @@ -164,7 +164,11 @@ SUBROUTINE trcbio er(jk,16) = correct_fact * ( gdept(jpk,jj,ji)-gdept(jk,jj,ji) ) /gdept(jpk,jj,ji) enddo #endif - call BFM1D_Input_EcologyDynamics(bottom,a,jtrmax,er) + + 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() 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..af522d2f 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 @@ -43,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/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_mem.f90 b/src/IO/BC_mem.f90 index 0cae6af8..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) @@ -114,7 +114,7 @@ SUBROUTINE alloc_DTATRC jn_gib = 6 ! 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 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/bc_atm_Hg0.f90 b/src/IO/bc_atm_Hg0.f90 new file mode 100644 index 00000000..c7f6ca7d --- /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,'hg0atm',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 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 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