DOM flow in Mosart

src/cpl/nuopc/rof_import_export.F90

@@ -11,7 +11,7 @@ module rof_import_export
   use shr_sys_mod     , only : shr_sys_abort
   use nuopc_shr_methods , only : chkerr
   use RunoffMod       , only : rtmCTL, TRunoff, TUnit
-  use RtmVar          , only : iulog, nt_rtm, rtm_tracers
+  use RtmVar          , only : iulog, nt_rtm, rtm_tracers,nt_rtm_dom,rtm_tracers_dom
   use RtmSpmd         , only : masterproc, mpicom_rof
   use RtmTimeManager  , only : get_nstep
   use nuopc_shr_methods , only : chkerr
@@ -110,6 +110,8 @@ subroutine advertise_fields(gcomp, flds_scalar_name, do_rtmflood, rc)
 
     call fldlist_add(fldsToRof_num, fldsToRof, trim(flds_scalar_name))
     call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofsur')
+    call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_subdoc')
+    call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_surfdoc')
     call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofgwl')
     call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofsub')
     call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofi')
@@ -244,7 +246,7 @@ subroutine import_fields( gcomp, rc )
     type(ESMF_State) :: importState
     integer          :: n,nt
     integer          :: begr, endr
-    integer          :: nliq, nfrz
+    integer          :: nliq, nfrz, ndoc ! ndoc is carbon tracer, ndon can be added for nitrogen later
     character(len=*), parameter :: subname='(rof_import_export:import_fields)'
     !---------------------------------------------------------------------------
 
@@ -267,6 +269,15 @@ subroutine import_fields( gcomp, rc )
        call shr_sys_abort()
     endif
 
+    ndoc = 0
+    do nt = 1,nt_rtm_dom
+       if (trim(rtm_tracers_dom(nt)) == 'DOC') ndoc = nt ! DOC is currently the only tracer in DOM
+    enddo
+    if (ndoc == 0) then
+       write(iulog,*) trim(subname),': ERROR in rtm_tracers_dom DOC ',ndoc,rtm_tracers_dom
+       call shr_sys_abort()
+    endif
+
     begr = rtmCTL%begr
     endr = rtmCTL%endr
 
@@ -277,6 +288,14 @@ subroutine import_fields( gcomp, rc )
          do_area_correction=.true., rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+    call state_getimport(importState, 'Flrl_surfdoc', begr, endr, rtmCTL%area, output=rtmCTL%domsur(:,ndoc), &
+         do_area_correction=.true., rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    call state_getimport(importState, 'Flrl_subdoc', begr, endr, rtmCTL%area, output=rtmCTL%domsub(:,ndoc), &
+         do_area_correction=.true., rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
     call state_getimport(importState, 'Flrl_rofsub', begr, endr, rtmCTL%area, output=rtmCTL%qsub(:,nliq), &
          do_area_correction=.true., rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return

src/riverroute/DommasbMod.F90

@@ -0,0 +1,79 @@
+MODULE DommasbMod
+  !Description: core code of Dissolved Organic Matter mass balance utilizing river routing models
+  !Developed by Marius Lambert 02-02-2023
+  !This module is currently made interact with MOSART routing model
+  ! USES:
+  use shr_kind_mod  , only : r8 => shr_kind_r8
+  use shr_const_mod , only : SHR_CONST_REARTH, SHR_CONST_PI
+  use shr_sys_mod   , only : shr_sys_abort
+  use RunoffMod     , only : TRunoff, Tdom, TUnit
+  use RtmVar        , only : iulog
+
+  implicit none
+
+  public hillslopeRoutingDOM
+  public subnetworkRoutingDOM
+  public mainchannelRoutingDOM
+  !--------------------------------------------------------------------
+
+  ! ! PUBLIC MEMBER FUNCTIONS:
+  contains
+
+  !----------------------------------------------------------------------
+    subroutine hillslopeRoutingDOM(iunit,nt,ntdom,theDeltaT)
+      ! ! DESCRIPTION: solve the ODEs with Euler algorithm for hillslope routing
+      implicit none
+      integer,  intent(in) :: iunit, nt, ntdom
+      real(r8), intent(in) :: theDeltaT
+      !domsur (kg/m2*s) ,domH  (kg/m2), ehout (m/s), domHout (kg/m2*s), qsur (m/s), wh (m)
+      Tdom%domHout(iunit,ntdom) = -TRunoff%ehout(iunit,nt) * (Tdom%domH(iunit,ntdom) + Tdom%domsur(iunit,ntdom) * theDeltaT)/(TRunoff%wh(iunit,nt)-TRunoff%dwh(iunit,nt)*theDeltaT+TRunoff%qsur(iunit,nt)*theDeltaT)
+
+      !we dont want a too high out
+      !Tdom%domHout(iunit,ntdom) = min(-TRunoff%ehout(iunit,nt) * 0.3_r8, Tdom%domHout(iunit,ntdom))
+      !cannot be less than 0, lower boundary
+      !Tdom%domHout(iunit,ntdom) = max(0._r8,Tdom%domHout(iunit,ntdom))
+      !cannot be more than available carbon, upper boundary
+      !Tdom%domHout(iunit,ntdom) = min((Tdom%domH(iunit,ntdom)+Tdom%domsur(iunit,ntdom)*theDeltaT)/theDeltaT,Tdom%domHout(iunit,ntdom))
+      !Tdom%domH(iunit,ntdom) = max(0._r8,Tdom%domH(iunit,ntdom) + (Tdom%domsur(iunit,ntdom) - Tdom%domHout(iunit,ntdom))* theDeltaT)
+
+      Tdom%domH(iunit,ntdom) = Tdom%domH(iunit,ntdom) + (Tdom%domsur(iunit,ntdom) - Tdom%domHout(iunit,ntdom))* theDeltaT
+
+    end subroutine hillslopeRoutingDOM
+
+    subroutine subnetworkRoutingDOM(iunit,nt,ntdom,theDeltaT)
+      ! solve the ODEs with Euler algorithm for subnetwork routing
+      implicit none
+      integer, intent(in) :: iunit, nt, ntdom
+      real(r8), intent(in) :: theDeltaT
+      ! domTout (kg/s), etout (m3/s), domT (kg), domsub (kg/s), domHout (kg/s), wt (m3), dwt (m3/s), etin (m3/s)
+
+      Tdom%domTout(iunit,ntdom) = -TRunoff%etout(iunit,nt) * (Tdom%domT(iunit,ntdom) + (Tdom%domsub(iunit,ntdom)+Tdom%domHout(iunit,ntdom)) * theDeltaT)/(TRunoff%wt(iunit,nt)-TRunoff%dwt(iunit,nt)*theDeltaT+TRunoff%etin(iunit,nt)*theDeltaT)
+
+      !Tdom%domTout(iunit,ntdom) = min(-TRunoff%etout(iunit,nt) *0.3_r8,Tdom%domTout(iunit,ntdom))
+      !Tdom%domTout(iunit,ntdom) = max(0._r8,Tdom%domTout(iunit,ntdom))
+      !Tdom%domTout(iunit,ntdom) = min((Tdom%domT(iunit,ntdom)+(Tdom%domsub(iunit,ntdom)+Tdom%domHout(iunit,ntdom))* theDeltaT)/theDeltaT,Tdom%domTout(iunit,ntdom))      
+      !Tdom%domT(iunit,ntdom) = max(0._r8,Tdom%domT(iunit,ntdom) + ( Tdom%domsub(iunit,ntdom) + Tdom%domHout(iunit,ntdom) - Tdom%domTout(iunit,ntdom) ) * theDeltaT)
+
+      Tdom%domT(iunit,ntdom) = Tdom%domT(iunit,ntdom) + ( Tdom%domsub(iunit,ntdom) + Tdom%domHout(iunit,ntdom) - Tdom%domTout(iunit,ntdom) ) * theDeltaT
+
+    end subroutine subnetworkRoutingDOM
+
+    subroutine mainchannelRoutingDOM(iunit,nt,ntdom,theDeltaT)
+      ! solve the ODE with Euler algorithm for main-channel routing
+      implicit none
+      integer, intent(in) :: iunit, nt, ntdom
+      real(r8), intent(in) :: theDeltaT
+      real(r8) :: temp_gwl
+      temp_gwl = TRunoff%qgwl(iunit,nt) * TUnit%area(iunit) * TUnit%frac(iunit)
+
+      Tdom%domRout(iunit,ntdom) = -TRunoff%erout(iunit,nt) * (Tdom%domR(iunit,ntdom) + (Tdom%domRUp(iunit,ntdom) + Tdom%domToutLat(iunit,ntdom)) * theDeltaT)/(TRunoff%wr(iunit,nt)-TRunoff%dwr(iunit,nt)*theDeltaT+(TRunoff%erlateral(iunit,nt)+TRunoff%erin(iunit,nt)+temp_gwl)*theDeltaT)
+
+      !Tdom%domRout(iunit,ntdom) = min(-TRunoff%erout(iunit,nt) *0.3_r8,Tdom%domRout(iunit,ntdom))
+      !Tdom%domRout(iunit,ntdom) = max(0._r8,Tdom%domRout(iunit,ntdom))
+      !Tdom%domRout(iunit,ntdom) = min((Tdom%domR(iunit,ntdom)+(Tdom%domRUp(iunit,ntdom) + Tdom%domToutLat(iunit,ntdom))* theDeltaT)/theDeltaT,Tdom%domRout(iunit,ntdom))
+      !Tdom%domR(iunit,ntdom) = max(0._r8,Tdom%domR(iunit,ntdom) + (Tdom%domRUp(iunit,ntdom) + Tdom%domToutLat(iunit,ntdom) - Tdom%domRout(iunit,ntdom)) * theDeltaT)
+
+      Tdom%domR(iunit,ntdom) = Tdom%domR(iunit,ntdom) + (Tdom%domRUp(iunit,ntdom) + Tdom%domToutLat(iunit,ntdom) - Tdom%domRout(iunit,ntdom)) * theDeltaT
+    end subroutine mainchannelRoutingDOM
+!-------------------------------------------------------------------------
+end MODULE DommasbMod