75 add t doc

src/CMakeLists.txt

@@ -20,6 +20,7 @@ add_library(fabm_models_ersem OBJECT
             mesozooplankton.F90
             bacteria.F90
             bacteria_docdyn.F90
+            TDOC.F90
             calcification.F90
             carbonate.F90
             oxygen.F90

src/TDOC.F90

@@ -0,0 +1,373 @@
+#include "fabm_driver.h"
+
+module ersem_TDOC
+
+   use fabm_types
+   use fabm_builtin_models
+   use fabm_particle
+
+   use ersem_shared
+   use ersem_pelagic_base
+
+   implicit none
+
+   private
+
+   type,extends(type_ersem_pelagic_base),public :: type_ersem_TDOC
+      ! Variables
+      type (type_state_variable_id) :: id_RPc,id_T2c, id_O3c
+      type (type_state_variable_id) :: id_RPn, id_RPp, id_T2n, id_T2p, id_N1p, id_N4n
+      type (type_state_variable_id) :: id_TD_older_c,id_TD_older_n,id_TD_older_p
+      type (type_model_id) :: id_T2, id_TD_older
+      type (type_dependency_id)                     :: id_ETW, id_chemEIR
+      type (type_dependency_id) ::   id_shadow_bioflux, id_dz
+      type (type_dependency_id)            :: id_X1X             ! Salinity
+!      type (type_horizontal_dependency_id)          :: id_R1c,id_R1d
+      type (type_horizontal_diagnostic_variable_id) :: id_surface_photolysis, id_surface_photo_aging
+      type (type_diagnostic_variable_id) :: id_photolysis,id_flocc, id_photo_aging, id_bio_aging
+
+   ! Parameters
+
+      real(rk) :: suva, iref,phyref,surf_phyref,phyt,floc,qp,qn,scx,sbx,chemEIR_scaling
+      real(rk) :: age, photoaging,bioaging
+      logical  :: is_photolabile
+
+   contains
+!     Model procedures
+      procedure :: initialize
+      procedure :: do_surface
+      procedure :: do
+   end type
+
+   ! Submodel that is seen by bacteria for microbial degradation
+   ! this is currently needed because the microbial degradation flux (calculated by the bacteria module) is needed here for the aging of tDOC
+   type,extends(type_particle_model),public :: type_ersem_shadow_tDOC
+      type (type_diagnostic_variable_id) :: id_c_shadow   ! shadow state variable
+      type (type_state_variable_id) :: id_parent_c, id_TD_older_parent_c ! ,id_parent_n, id_T1_older_parent_n,id_parent_p, id_T1_older_parent_p
+      type (type_model_id) :: id_parent, id_TD_older_parent
+
+      real (rk)  :: qn, qp, bioaging
+
+   contains
+      procedure :: initialize => shadow_initialize
+      procedure :: do  => shadow_do
+   end type
+
+
+contains
+
+   subroutine initialize(self,configunit)
+!
+! !DESCRIPTION:
+!
+! !INPUT PARAMETERS:
+      class (type_ersem_TDOC),intent(inout),target :: self
+      integer,                         intent(in)           :: configunit
+
+! !LOCAL VARIABLES:
+      real(rk) :: c0,EPS
+      class (type_ersem_shadow_tDOC), pointer :: shadow      
+!
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+!      call self%get_parameter(self%kuw,'kuw', 'm-1',              'extinction of UV in the water')
+!      call self%get_parameter(self%suva,'suva','m2*mmol C-1',            'specific UV absorption at 350 nm')
+      call self%get_parameter(self%qp,    'qp',    'mmol P/mg C','phosphorus to carbon ratio')
+      call self%get_parameter(self%qn,    'qn',    'mmol N/mg C','nitrogen to carbon ratio')
+      call self%get_parameter(EPS,     'EPS',    'm^2/mg C','specific shortwave attenuation', default=4.E-4_rk)
+      call self%get_parameter(c0,'c0','mg C/m^3','background carbon concentration')
+      !call self%get_parameter(self%age, 'age', 'd', 'Characteristic age of the tDOC pool', default=0._rk))  ! for now the parameter is calculated a priori
+      call self%get_parameter(self%photoaging, 'photoaging', '-', 'aging due to photo oxidation')
+      call self%get_parameter(self%bioaging, 'bioaging', '-', 'aging due to microbial oxidation')
+      call self%get_parameter(self%is_photolabile,'photolabile','-',            'definition of the type of tDOC: True=photolabile; False=non-photolabile')
+
+      ! if it is photolabile then set all process parameters
+      if (self%is_photolabile) then
+         call self%get_parameter(self%chemEIR_scaling, 'chemEIR_scaling', '', 'scaling factor for incoming radiation activating photochemical reaction', default=1._rk)
+         call self%get_parameter(self%iref,'iref','W m-2',            'reference irradiance')
+         call self%get_parameter(self%phyref,'phyref','d-1 ',            'reference photooxidation rate',default=0._rk)
+         call self%get_parameter(self%surf_phyref,'surf_phyref','d-1 ',            'reference surface_photooxidation rate',default=0._rk)
+         call self%get_parameter(self%phyt,'phyt','adim',            'photooxidated fraction of T1 going into T2) ')
+         call self%get_parameter(self%floc,'floc','(mmol C-3)-1*d-1 ',            'reference photooxidation rate')
+         call self%get_parameter(self%sbx,   'sbx',   'psu',        'optimal salinity')
+         call self%get_parameter(self%scx,   'scx',   '',        'salinity function parameter')
+      endif
+
+! Allow ERSEM base model to declare our own state variables.
+      call self%initialize_ersem_base(sedimentation=.false.)
+      call self%add_constituent('c',1.e-4_rk,c0,qn=self%qn,qp=self%qp)
+!      call self%add_constituent('c',0.0_rk)
+
+! Register links to nutrient pools.
+!      call self%register_dependency(self%id_EIR,'EIR','W/m^2','downwelling_shortwave_flux', &
+!              standard_variable=standard_variables%downwelling_shortwave_flux,source=source_do_column)
+      !call self%register_dependency(self%id_EIR,standard_variables%downwelling_shortwave_flux)
+
+      if (self%is_photolabile) then
+          call self%register_dependency(self%id_dz, standard_variables%cell_thickness)
+          call self%register_dependency(self%id_X1X,standard_variables%practical_salinity)
+          call self%register_state_dependency(self%id_RPc,'RPc','mg C/m^3',   'particulate organic carbon')
+          call self%register_state_dependency(self%id_RPp,'RPp','mmol P/m^3', 'particulate organic phosphorus')
+          call self%register_state_dependency(self%id_RPn,'RPn','mmol N/m^3', 'particulate organic nitrogen')
+          call self%register_model_dependency(self%id_T2,'T2')
+          call self%register_state_dependency(self%id_T2c,'T2c','mg C/m^3','non photolabile terrigenous DOC')
+          call self%register_state_dependency(self%id_T2n,'T2n','mmol N/m^3','non photolabile terrigenous DON')
+          call self%register_state_dependency(self%id_T2p,'T2p','mmol P/m^3','non photolabile terrigenous DOP')
+          call self%request_coupling_to_model(self%id_T2c,self%id_T2,'c')
+          call self%request_coupling_to_model(self%id_T2n,self%id_T2,standard_variables%total_nitrogen)
+          call self%request_coupling_to_model(self%id_T2p,self%id_T2,standard_variables%total_phosphorus)
+          call self%register_dependency(self%id_chemEIR,'chemEIR','W/m^2','incoming radiation activating photochemical reaction')
+          call self%register_state_dependency(self%id_O3c,'O3c','mmol C/m^3','carbon dioxide sink')
+          call self%register_state_dependency(self%id_N1p,'N1p','mmol P/m^3','phosphate')
+          call self%register_state_dependency(self%id_N4n,'N4n','mmol N/m^3','ammonium')
+          call self%register_diagnostic_variable(self%id_photolysis,'photolysis','mgC/m^3/d','photolysis')
+          call self%register_horizontal_diagnostic_variable(self%id_surface_photolysis,'surface_photolysis','mgC/m^3/d','surface photolysis',source=source_do_surface)
+          call self%register_diagnostic_variable(self%id_flocc,'flocc','mgC/m^3/d','flocculation')
+      end if
+
+      if ((self%bioaging.gt.0._rk).or.(self%photoaging.gt.0._rk)) then
+          call self%register_model_dependency(self%id_TD_older,'TD_older')
+          call self%register_state_dependency(self%id_TD_older_c,'TD_older_c','mg C/m^3','non photolabile terrigenous DOC')
+          call self%register_state_dependency(self%id_TD_older_n,'TD_older_n','mmol N/m^3','non photolabile terrigenous DON')
+          call self%register_state_dependency(self%id_TD_older_p,'TD_older_p','mmol P/m^3','non photolabile terrigenous DOP')
+          call self%request_coupling_to_model(self%id_TD_older_c,self%id_TD_older,'c')
+          call self%request_coupling_to_model(self%id_TD_older_n,self%id_TD_older,standard_variables%total_nitrogen)
+          call self%request_coupling_to_model(self%id_TD_older_p,self%id_TD_older,standard_variables%total_phosphorus)
+          if (self%is_photolabile) then 
+             call self%register_diagnostic_variable(self%id_photo_aging,'photoaging','mgC/m^3/d','aging due to 3D photoloysis')
+             call self%register_horizontal_diagnostic_variable(self%id_surface_photo_aging,'surface_photoaging','mgC/m^3/d','aging due to surface photolysis',source=source_do_surface)
+          end if
+          call self%register_diagnostic_variable(self%id_bio_aging,'bioaging','mgC/m^3/d','aging due to microbial degrdation')      
+
+
+    end if
+
+      ! Register contribution to light extinction
+      call self%add_to_aggregate_variable(standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux, &
+         self%id_c,scale_factor=EPS,include_background=.true.)
+
+
+      allocate(shadow)
+      shadow%qn=self%qn
+      shadow%qp=self%qp
+      shadow%bioaging=self%bioaging
+      call shadow%couplings%set_string('parent',self%name)
+      call self%add_child(shadow,'shadow',configunit=configunit)      
+      if ((self%bioaging.gt.0._rk)) then 
+          call shadow%couplings%set_string('TD_older_parent','../TD_older')
+          !call self%register_dependency(self%id_shadow_bioflux,'B_degradation','','bacterial degradation of tDOC')
+          !call self%request_coupling(self%id_shadow_bioflux,trim(shadow%id_c_shadow%link%target%name//'_sms_tot')) !
+          !call copy_fluxes(self,shadow%id_c_shadow,self%id_c,(1._rk+self%bioaging))
+      endif
+
+   end subroutine initialize
+
+
+   subroutine do(self,_ARGUMENTS_DO_)
+
+      class (type_ersem_TDOC),intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_
+
+   ! !LOCAL VARIABLES:
+    real(rk) :: flocc,photolysis,R8cP,R8c,T2c,T2n,T2p,O3c,O3cP,T2cP,chemEIR,c,Xp,Xn,XXn,T1T2,px,nx, bio_flux
+    real(rk) :: X1X,sal
+
+! Enter spatial loops (if any)
+      _LOOP_BEGIN_
+
+    if (self%is_photolabile) then
+      _GET_(self%id_c,c)
+      _GET_(self%id_T2c,T2cP)
+!      _GET_(self%id_R8c,R8cP)
+      _GET_(self%id_O3c,O3cP)
+      _GET_WITH_BACKGROUND_(self%id_T2c,T2c)
+      _GET_WITH_BACKGROUND_(self%id_T2n,T2n)
+      _GET_WITH_BACKGROUND_(self%id_T2p,T2p)
+!      _GET_WITH_BACKGROUND_(self%id_R8c,R8c)
+      _GET_(self%id_chemEIR,chemEIR)
+      _GET_(self%id_X1X,X1X)
+
+      IF(T2c.gt.0._rk) then
+      Xn=self%qn/(T2n/T2c)
+      Xp=self%qp/(T2p/T2c)
+      nx=self%qn-(T2n/T2c)
+      px=self%qp-(T2p/T2c)
+      else
+      Xn=0._rk
+      Xp=0._rk
+      nx=1._rk
+      px=1._rk
+      endif
+
+      XXn=min(Xn,Xp)
+
+!Flocculation is assumed to be (log-normal) function of salinity(sal)
+
+     X1X=max(X1X,0.01_rk) !to avoid log(0) in sal!!
+     sal=exp(-((log(X1X)-self%sbx)**2._rk)/(2._rk*self%scx**2._rk))
+     flocc=self%floc*sal*c**2
+
+! Photolysis
+     photolysis=self%phyref*(chemEIR*self%chemEIR_scaling/self%iref)*c
+
+! in order to ensure mass conservation, the fraction of carbon going from T1 to T2 after photolysis is down regulated 
+! if T1 has less  N and/or P than T2. Any excess of nutrients is redirected into the dissolved pool (N1p and N4n). (Luca, July 2018)    
+     T1T2=self%phyt*min(1._rk,XXn)
+
+     _SET_ODE_(self%id_c,-(1._rk+self%photoaging)*photolysis-flocc)
+
+     if (self%photoaging.gt.0._rk) then
+        _SET_ODE_(self%id_TD_older_c,self%photoaging*photolysis)
+!         _SET_ODE_(self%id_TD_older_n,self%photoaging*photolysis*self%qn)
+!         _SET_ODE_(self%id_TD_older_p,self%photoaging*photolysis*self%qp)
+     end if
+
+     !if (self%bioaging.gt.0._rk) then
+     !   _GET_(self%id_shadow_bioflux,bio_flux)
+     !   _SET_ODE_(self%id_TD_older_c,self%bioaging*bio_flux)
+     !   _SET_DIAGNOSTIC_(self%id_bio_aging,self%bioaging*bio_flux)
+     !endif
+
+      _SET_ODE_(self%id_RPc,+flocc)
+      _SET_ODE_(self%id_RPn,+flocc*self%qn)
+      _SET_ODE_(self%id_RPp,+flocc*self%qp)
+      _SET_ODE_(self%id_T2c,+photolysis*T1T2)
+      _SET_ODE_(self%id_O3c,+photolysis*(1._rk-T1T2)/CMass)
+      _SET_ODE_(self%id_N1p,+photolysis*(1._rk-T1T2)*self%qp+photolysis*T1T2*px)
+      _SET_ODE_(self%id_N4n,+photolysis*(1._rk-T1T2)*self%qn+photolysis*T1T2*nx)
+      _SET_DIAGNOSTIC_(self%id_photolysis, photolysis)
+      _SET_DIAGNOSTIC_(self%id_flocc, flocc)
+
+   end if !is_photolabile
+   _LOOP_END_
+
+   end subroutine do
+
+   subroutine do_surface(self,_ARGUMENTS_DO_SURFACE_)
+      class (type_ersem_TDOC), intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_SURFACE_
+
+         real (rk) :: chemEIR,photolysis,c,T1T2,T2c,T2n,T2p
+         real (rk) :: Xn,Xp,XXn,nx,px,dz
+
+! Enter horizontal loops (if any)
+      _HORIZONTAL_LOOP_BEGIN_
+
+      if (self%is_photolabile) then
+         _GET_(self%id_chemEIR,chemEIR)
+         _GET_(self%id_c,c)
+         _GET_(self%id_dz,dz)
+         _GET_WITH_BACKGROUND_(self%id_T2c,T2c)
+         _GET_WITH_BACKGROUND_(self%id_T2n,T2n)
+         _GET_WITH_BACKGROUND_(self%id_T2p,T2p)
+
+         ! YA: this BLOCK is commented because does not work as it should
+         ! for now I leave it like that and force behaviour as if T1 stoichiometry is identical to T2 stoichiometry (True for LOCATE 3D runs)
+!          IF(T2c.gt.0._rk) then
+!             Xn=self%qn/(T2n/T2c)
+!             Xp=self%qp/(T2p/T2c)
+!             nx=self%qn-(T2n/T2c)
+!             px=self%qp-(T2p/T2c)
+!          else
+!             Xn=0._rk
+!             Xp=0._rk
+!             nx=1._rk
+!             px=1._rk
+!          endif
+! 
+!          ! the two following checks are needed because due to roundings, N and P are not conserved otherwise
+!          if (abs(nx).lt.1.e-5_rk) then
+!               nx=0._rk
+!               Xn=1._rk
+!          endif
+!          if (abs(px).lt.1.e-5_rk) then 
+!              px=0._rk
+!              Xp=1._rk
+!          endif
+
+              nx=0._rk
+              Xn=1._rk
+              px=0._rk
+              Xp=1._rk
+          XXn=min(Xn,Xp)
+         ! Photolysis
+         ! this is the real photolitic rate in mgC/m3/d
+         photolysis=self%surf_phyref*(chemEIR*self%chemEIR_scaling/self%iref)*c
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_surface_photolysis, photolysis)
+
+         ! because this is a do_surface subroutine only _set_surface_exchange_ can be used, not _set_ode_
+         ! _set_surface_exchange_ requires the flux being vertically integrated (units = 1/m2/d)
+         ! so the flux is multiplied by the cell thickness. _set_surface_exchange_ will re-convert back the unit automatically
+
+         photolysis = photolysis * dz
+
+         ! in order to ensure mass conservation, the fraction of carbon going from T1 to T2 after photolysis is down regulated 
+         ! if T1 has less  N and/or P than T2. Any excess of nutrients is redirected into the dissolved pool (N1p and N4n). (Luca, July 2018)    
+         T1T2=self%phyt*min(1._rk,XXn)
+
+
+         _SET_SURFACE_EXCHANGE_(self%id_c,-photolysis*(1+self%photoaging))
+         _SET_SURFACE_EXCHANGE_(self%id_T2c,+photolysis*T1T2)
+
+         if (self%photoaging.gt.0._rk) then
+            _SET_SURFACE_EXCHANGE_(self%id_TD_older_c,self%photoaging*photolysis)
+!             _SET_SURFACE_EXCHANGE_(self%id_TD_older_n,-self%photoaging*photolysis*self%qn)
+!             _SET_SURFACE_EXCHANGE_(self%id_TD_older_p,-self%photoaging*photolysis*self%qp)
+         end if         
+
+         _SET_SURFACE_EXCHANGE_(self%id_O3c,+photolysis*(1._rk-T1T2)/CMass)
+         if (nx.NE.0._rk) write(6,*) trim(self%name),T2c,T2n,self%qn,Xn
+         _SET_SURFACE_EXCHANGE_(self%id_N1p,+photolysis*(1._rk-T1T2)*self%qp)!+photolysis*T1T2*px)
+         _SET_SURFACE_EXCHANGE_(self%id_N4n,+photolysis*(1._rk-T1T2)*self%qn)!+photolysis*T1T2*nx)
+    end if !is_photolabile
+      _HORIZONTAL_LOOP_END_    
+    end subroutine
+
+   subroutine shadow_initialize(self,configunit)
+
+   ! !INPUT PARAMETERS:
+      class (type_ersem_shadow_TDOC),intent(inout),target :: self
+      integer,                       intent(in)           :: configunit
+
+      self%dt = 3600._rk*24._rk
+
+      call self%register_diagnostic_variable(self%id_c_shadow,'c','mgC/m^3','shadow carbon',act_as_state_variable=.true.,output=output_none)
+      call self%add_to_aggregate_variable(standard_variables%total_carbon,self%id_c_shadow)
+      call self%add_to_aggregate_variable(standard_variables%total_phosphorus,self%id_c_shadow,scale_factor=self%qp)
+      call self%add_to_aggregate_variable(standard_variables%total_nitrogen,self%id_c_shadow,scale_factor=self%qn)
+
+      call self%register_model_dependency(self%id_parent,'parent')
+      call self%register_state_dependency(self%id_parent_c,'parent_c','mg C/m^3','non photolabile terrigenous DOC')
+      call self%request_coupling_to_model(self%id_parent_c,self%id_parent,'c')
+
+      if (self%bioaging.gt.0._rk) then
+          call self%register_model_dependency(self%id_TD_older_parent,'TD_older_parent')
+          call self%register_state_dependency(self%id_TD_older_parent_c,'TD_older_parent_c','mg C/m^3','non photolabile terrigenous DOC')
+          call self%request_coupling_to_model(self%id_TD_older_parent_c,self%id_TD_older_parent,'c')
+          call copy_fluxes(self,self%id_c_shadow,self%id_TD_older_parent_c,-self%bioaging)  ! - because an uptake (negative flux) from bacteria correspond to increase of the older DOC (positive fulx)
+
+      end if
+
+      call copy_fluxes(self,self%id_c_shadow,self%id_parent_c,1._rk+self%bioaging)!scale_factor=(1._rk+self%bioaging))
+
+   end subroutine shadow_initialize
+
+   subroutine shadow_do(self,_ARGUMENTS_DO_)
+
+      class (type_ersem_shadow_TDOC),intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_
+
+         real (rk) :: c_shadow
+
+! Enter horizontal loops (if any)
+      _LOOP_BEGIN_
+
+      _GET_(self%id_parent_c,c_shadow)
+      _SET_DIAGNOSTIC_(self%id_c_shadow,c_shadow)
+
+      _LOOP_END_    
+    end subroutine shadow_do
+
+
+end module