diff --git a/cime_config/__pycache__/buildnmlcpython-39.pyc b/cime_config/__pycache__/buildnmlcpython-39.pyc new file mode 100644 index 0000000..9f8c2af Binary files /dev/null and b/cime_config/__pycache__/buildnmlcpython-39.pyc differ diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 30fe4fb..408c814 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/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 diff --git a/src/riverroute/DommasbMod.F90 b/src/riverroute/DommasbMod.F90 new file mode 100755 index 0000000..ac2fcee --- /dev/null +++ b/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 diff --git a/src/riverroute/MOSART_physics_mod.F90 b/src/riverroute/MOSART_physics_mod.F90 index a2d327f..8b1af83 100644 --- a/src/riverroute/MOSART_physics_mod.F90 +++ b/src/riverroute/MOSART_physics_mod.F90 @@ -13,18 +13,20 @@ MODULE MOSART_physics_mod 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 RtmVar , only : iulog, barrier_timers, nt_rtm, rtm_tracers - use RunoffMod , only : Tctl, TUnit, TRunoff, TPara, rtmCTL + use RtmVar , only : iulog, barrier_timers, nt_rtm, rtm_tracers, nt_rtm_dom, rtm_tracers_dom + use RunoffMod , only : Tctl, TUnit, TRunoff, TPara, rtmCTL,Tdom use RunoffMod , only : SMatP_eroutUp, avsrc_eroutUp, avdst_eroutUp + use RunoffMod , only : SMatP_domRUp, avsrc_domRUp, avdst_domRUp ! Matrix to calculate the DOC coming from all upstream gridcells use RtmSpmd , only : masterproc, mpicom_rof use perf_mod , only: t_startf, t_stopf use mct_mod + use DommasbMod implicit none private real(r8), parameter :: TINYVALUE = 1.0e-14_r8 ! double precision variable has a significance of about 16 decimal digits - integer :: nt ! loop indices + integer :: nt, ntdom ! loop indices real(r8), parameter :: SLOPE1def = 0.1_r8 ! here give it a small value in order to avoid the abrupt change of hydraulic radidus etc. real(r8) :: sinatanSLOPE1defr ! 1.0/sin(atan(slope1)) @@ -49,33 +51,74 @@ subroutine Euler integer :: iunit, m, k, unitUp, cnt, ier !local index real(r8) :: temp_erout, localDeltaT real(r8) :: negchan - + real(r8) :: temp_eroutdom(nt_rtm_dom) + real(r8) :: liqsub(rtmCTL%begr:rtmCTL%endr,nt_rtm),concsub(rtmCTL%begr:rtmCTL%endr,nt_rtm_dom) + real(r8) :: liqhill(rtmCTL%begr:rtmCTL%endr,nt_rtm),conchill(rtmCTL%begr:rtmCTL%endr,nt_rtm_dom) !------------------ ! hillslope !------------------ - + liqsub(:,:)=0._r8 + concsub(:,:)=0._r8 + liqhill(:,:)=0._r8 + conchill(:,:)=0._r8 + call t_startf('mosartr_hillslope') do nt=1,nt_rtm if (TUnit%euler_calc(nt)) then do iunit=rtmCTL%begr,rtmCTL%endr - if(TUnit%mask(iunit) > 0) then + if (TUnit%mask(iunit) > 0) then call hillslopeRouting(iunit,nt,Tctl%DeltaT) TRunoff%wh(iunit,nt) = TRunoff%wh(iunit,nt) + TRunoff%dwh(iunit,nt) * Tctl%DeltaT call UpdateState_hillslope(iunit,nt) TRunoff%etin(iunit,nt) = (-TRunoff%ehout(iunit,nt) + TRunoff%qsub(iunit,nt)) * TUnit%area(iunit) * TUnit%frac(iunit) + !----------------------------------------------------------------------------------------------------------------- + if (nt==1) then ! if LIQ tracer + do ntdom=1,nt_rtm_dom ! loop over DOM tracers + liqhill(iunit,nt)=TRunoff%wh(iunit,nt)-TRunoff%dwh(iunit,nt)*Tctl%DeltaT+TRunoff%qsur(iunit,nt)*Tctl%DeltaT + conchill(iunit,ntdom)=(Tdom%domH(iunit,ntdom) + Tdom%domsur(iunit,ntdom) * Tctl%DeltaT)/liqhill(iunit,nt) + if (liqhill(iunit,nt)>0._r8 .and. conchill(iunit,ntdom)<=0.3_r8 .and. conchill(iunit,ntdom)>0._r8) then + if (TRunoff%wh(iunit,nt)-TRunoff%dwh(iunit,nt)*Tctl%DeltaT < 0._r8) then + Tdom%domResthill(iunit,ntdom) = Tdom%domResthill(iunit,ntdom) + Tdom%domH(iunit,ntdom) + Tdom%domH(iunit,ntdom)=0._r8 + endif + Tdom%domHout(iunit,ntdom) = -TRunoff%ehout(iunit,nt) * conchill(iunit,ntdom) + Tdom%domH(iunit,ntdom) = Tdom%domH(iunit,ntdom) + (Tdom%domsur(iunit,ntdom) - Tdom%domHout(iunit,ntdom)) * Tctl%DeltaT + else if (Tdom%domsur(iunit,ntdom)>0._r8) then + Tdom%domResthill(iunit,ntdom)= Tdom%domResthill(iunit,ntdom)+Tdom%domsur(iunit,ntdom)*Tctl%DeltaT + endif + if (TRunoff%wh(iunit,nt)< 0._r8) then + Tdom%domResthill(iunit,ntdom)= Tdom%domResthill(iunit,ntdom)+Tdom%domH(iunit,ntdom) + Tdom%domH(iunit,ntdom)=0._r8 + endif + if (Tdom%domH(iunit,ntdom) < 0._r8) then + Tdom%domResthill(iunit,ntdom)= Tdom%domResthill(iunit,ntdom)+Tdom%domH(iunit,ntdom) + Tdom%domH(iunit,ntdom)=0._r8 + endif + ! here some checks to make sure the DOM is not at too high or low concentrations + if (Tdom%domH(iunit,ntdom)/TRunoff%wh(iunit,nt) > 0.3_r8 .and. TRunoff%wh(iunit,nt) >= 0._r8) then + Tdom%domResthill(iunit,ntdom)= Tdom%domResthill(iunit,ntdom)+Tdom%domH(iunit,ntdom)-0.3*TRunoff%wh(iunit,nt) + Tdom%domH(iunit,ntdom)=0.3*TRunoff%wh(iunit,nt) + endif + Tdom%domsub(iunit,ntdom) = Tdom%domsub(iunit,ntdom) * TUnit%area(iunit) * TUnit%frac(iunit) ! readjust to correct units kg/m2s --> kg/s + Tdom%domHout(iunit,ntdom) = Tdom%domHout(iunit,ntdom) * TUnit%area(iunit) * TUnit%frac(iunit) ! readjust to correct units kg/m2s --> kg/s + Tdom%domResthill(iunit,ntdom) = Tdom%domResthill(iunit,ntdom) * TUnit%area(iunit) * TUnit%frac(iunit) ! readjust to correct units kg/m2 --> kg + enddo + endif + !-------------------------------------------------------------------------------------------------------------------------- endif end do endif end do call t_stopf('mosartr_hillslope') - TRunoff%flow = 0._r8 + Tdom%domRoutflow = 0._r8 + Tdom%domToutLat2 = 0._r8 + TRunoff%erlateral2 = 0._r8 TRunoff%erout_prev = 0._r8 TRunoff%eroutup_avg = 0._r8 TRunoff%erlat_avg = 0._r8 negchan = 9999.0_r8 do m=1,Tctl%DLevelH2R - !--- accumulate/average erout at prior timestep (used in eroutUp calc) for budget analysis do nt=1,nt_rtm if (TUnit%euler_calc(nt)) then @@ -91,6 +134,7 @@ subroutine Euler call t_startf('mosartr_subnetwork') TRunoff%erlateral(:,:) = 0._r8 + Tdom%domToutLat(:,:) = 0._r8 do nt=1,nt_rtm if (TUnit%euler_calc(nt)) then do iunit=rtmCTL%begr,rtmCTL%endr @@ -101,8 +145,45 @@ subroutine Euler TRunoff%wt(iunit,nt) = TRunoff%wt(iunit,nt) + TRunoff%dwt(iunit,nt) * localDeltaT call UpdateState_subnetwork(iunit,nt) TRunoff%erlateral(iunit,nt) = TRunoff%erlateral(iunit,nt)-TRunoff%etout(iunit,nt) + !---------------------------------------------------------------------------------------------------- + if (nt==1) then ! if LIQ tracer and there is water + do ntdom=1,nt_rtm_dom ! loop over DOM tracers + liqsub(iunit,nt)=TRunoff%wt(iunit,nt)-TRunoff%dwt(iunit,nt)*localDeltaT+TRunoff%etin(iunit,nt)*localDeltaT + concsub(iunit,ntdom)=(Tdom%domT(iunit,ntdom) + (Tdom%domsub(iunit,ntdom)+Tdom%domHout(iunit,ntdom)) * localDeltaT)/liqsub(iunit,nt) + if (liqsub(iunit,nt) > 0._r8 .and. concsub(iunit,ntdom)<=0.3_r8 .and. concsub(iunit,ntdom)>0._r8) then + if (TRunoff%wt(iunit,nt)-TRunoff%dwt(iunit,nt)*localDeltaT< 0._r8) then + Tdom%domRestsubn(iunit,ntdom) = Tdom%domRestsubn(iunit,ntdom) + Tdom%domT(iunit,ntdom) + Tdom%domT(iunit,ntdom)=0._r8 + endif + Tdom%domTout(iunit,ntdom) = -TRunoff%etout(iunit,nt)*concsub(iunit,ntdom) + Tdom%domT(iunit,ntdom) = Tdom%domT(iunit,ntdom) + ( Tdom%domsub(iunit,ntdom) + Tdom%domHout(iunit,ntdom) - Tdom%domTout(iunit,ntdom) ) * localDeltaT + else if ((Tdom%domsub(iunit,ntdom)+Tdom%domHout(iunit,ntdom))>0._r8) then + Tdom%domRestsubn(iunit,ntdom)= Tdom%domRestsubn(iunit,ntdom)+(Tdom%domsub(iunit,ntdom)+Tdom%domHout(iunit,ntdom))*localDeltaT + endif + if (TRunoff%wt(iunit,nt) < 0._r8) then + Tdom%domRestsubn(iunit,ntdom)= Tdom%domRestsubn(iunit,ntdom)+Tdom%domT(iunit,ntdom) + Tdom%domT(iunit,ntdom)=0._r8 + endif + if (Tdom%domT(iunit,ntdom) < 0._r8) then + Tdom%domRestsubn(iunit,ntdom)=Tdom%domRestsubn(iunit,ntdom)+Tdom%domT(iunit,ntdom) + Tdom%domT(iunit,ntdom)=0._r8 + endif + if (Tdom%domT(iunit,ntdom)/TRunoff%wt(iunit,nt) > 0.3_r8 .and. TRunoff%wt(iunit,nt) >= 0._r8) then + Tdom%domRestsubn(iunit,ntdom)=Tdom%domRestsubn(iunit,ntdom)+Tdom%domT(iunit,ntdom)-0.3*TRunoff%wt(iunit,nt) + Tdom%domT(iunit,ntdom)=0.3*TRunoff%wt(iunit,nt) + endif + Tdom%domToutLat(iunit,ntdom) = Tdom%domToutLat(iunit,ntdom) + Tdom%domTout(iunit,ntdom) + enddo + endif end do ! numDT_t TRunoff%erlateral(iunit,nt) = TRunoff%erlateral(iunit,nt) / TUnit%numDT_t(iunit) + TRunoff%erlateral2(iunit,nt) = TRunoff%erlateral2(iunit,nt) + TRunoff%erlateral(iunit,nt) + if (nt==1) then ! if LIQ tracer + do ntdom=1,nt_rtm_dom ! loop over DOM tracers + Tdom%domToutLat(iunit,ntdom) = Tdom%domToutLat(iunit,ntdom) / TUnit%numDT_t(iunit) + Tdom%domToutLat2(iunit,ntdom) = Tdom%domToutLat2(iunit,ntdom) + Tdom%domToutLat(iunit,ntdom) + end do + end if endif end do ! iunit endif ! euler_calc @@ -117,10 +198,11 @@ subroutine Euler call t_startf('mosartr_SMeroutUp_barrier') call mpi_barrier(mpicom_rof,ier) call t_stopf('mosartr_SMeroutUp_barrier') - endif + endif call t_startf('mosartr_SMeroutUp') TRunoff%eroutUp = 0._r8 + Tdom%domRUp = 0._r8 #ifdef NO_MCT do iunit=rtmCTL%begr,rtmCTL%endr do k=1,TUnit%nUp(iunit) @@ -133,16 +215,23 @@ subroutine Euler #else !--- copy erout into avsrc_eroutUp --- call mct_avect_zero(avsrc_eroutUp) + call mct_avect_zero(avsrc_domRUp) cnt = 0 do iunit = rtmCTL%begr,rtmCTL%endr cnt = cnt + 1 do nt = 1,nt_rtm avsrc_eroutUp%rAttr(nt,cnt) = TRunoff%erout(iunit,nt) + if (nt==1) then + do ntdom = 1,nt_rtm_dom + avsrc_domRUp%rAttr(ntdom,cnt) = Tdom%domRout(iunit,ntdom) + end do + endif enddo enddo call mct_avect_zero(avdst_eroutUp) - + call mct_avect_zero(avdst_domRUp) call mct_sMat_avMult(avsrc_eroutUp, sMatP_eroutUp, avdst_eroutUp) + call mct_sMat_avMult(avsrc_domRUp, sMatP_domRUp, avdst_domRUp) !--- add mapped eroutUp to TRunoff --- cnt = 0 @@ -150,6 +239,11 @@ subroutine Euler cnt = cnt + 1 do nt = 1,nt_rtm TRunoff%eroutUp(iunit,nt) = avdst_eroutUp%rAttr(nt,cnt) + if (nt==1) then + do ntdom = 1,nt_rtm_dom + Tdom%domRUp(iunit,ntdom) = avdst_domRUp%rAttr(ntdom,cnt) + end do + endif enddo enddo #endif @@ -157,7 +251,6 @@ subroutine Euler TRunoff%eroutup_avg = TRunoff%eroutup_avg + TRunoff%eroutUp TRunoff%erlat_avg = TRunoff%erlat_avg + TRunoff%erlateral - !------------------ ! channel routing !------------------ @@ -169,6 +262,7 @@ subroutine Euler if(TUnit%mask(iunit) > 0) then localDeltaT = Tctl%DeltaT/Tctl%DLevelH2R/TUnit%numDT_r(iunit) temp_erout = 0._r8 + temp_eroutdom(:) = 0._r8 do k=1,TUnit%numDT_r(iunit) call mainchannelRouting(iunit,nt,localDeltaT) TRunoff%wr(iunit,nt) = TRunoff%wr(iunit,nt) + TRunoff%dwr(iunit,nt) * localDeltaT @@ -179,10 +273,45 @@ subroutine Euler ! end if call UpdateState_mainchannel(iunit,nt) temp_erout = temp_erout + TRunoff%erout(iunit,nt) ! erout here might be inflow to some downstream subbasin, so treat it differently than erlateral + !----------------------------------------------------------------------------------------------------------- + if (nt==1) then ! if LIQ tracer and there is water + do ntdom=1,nt_rtm_dom ! loop over DOM tracers + if (TRunoff%wr(iunit,nt)-TRunoff%dwr(iunit,nt)*localDeltaT+(TRunoff%erlateral(iunit,nt)+TRunoff%erin(iunit,nt))*localDeltaT>0._r8) then + if (TRunoff%wr(iunit,nt)-TRunoff%dwr(iunit,nt)*localDeltaT < 0._r8) then + Tdom%domRestmain(iunit,ntdom)= Tdom%domRestmain(iunit,ntdom)+Tdom%domR(iunit,ntdom) + Tdom%domR(iunit,ntdom)=0._r8 + endif + call mainchannelRoutingDOM(iunit,nt,ntdom,localDeltaT) + temp_eroutdom(ntdom) = temp_eroutdom(ntdom) + Tdom%domRout(iunit,ntdom) + else if ((Tdom%domRUp(iunit,ntdom)+Tdom%domToutLat(iunit,ntdom))>0._r8) then + Tdom%domRestmain(iunit,ntdom)= Tdom%domRestmain(iunit,ntdom)+(Tdom%domRUp(iunit,ntdom)+Tdom%domToutLat(iunit,ntdom))*localDeltaT + endif + if (Tdom%domR(iunit,ntdom)/TRunoff%wr(iunit,nt) > 0.30001_r8) then !.and. Tdom%domR(iunit,ntdom) < 1.e-10*(Tdom%domRout(iunit,ntdom)+Tdom%domToutLat(iunit,ntdom)+Tdom%domRUp(iunit,ntdom))) then + Tdom%domRestmain(iunit,ntdom)=Tdom%domRestmain(iunit,ntdom)+Tdom%domR(iunit,ntdom) + Tdom%domR(iunit,ntdom)=0._r8 + endif + if (Tdom%domR(iunit,ntdom) < 1.e-50_r8) then + Tdom%domRestmain(iunit,ntdom)=Tdom%domRestmain(iunit,ntdom)+Tdom%domR(iunit,ntdom) + Tdom%domR(iunit,ntdom)=0._r8 + endif + !if (Tdom%domR(iunit,ntdom)/TRunoff%wr(iunit,nt) > 0.30001_r8 .or. Tdom%domR(iunit,ntdom) <0._r8) then + ! write(iulog,*)' Concentration in main is too high or too low ',Tdom%domR(iunit,ntdom),TRunoff%wr(iunit,nt) + !call shr_sys_abort('Concentration in main is too high or too low') + !endif + enddo + endif + !---------------------------------------------------------------------------------------------------------------- end do temp_erout = temp_erout / TUnit%numDT_r(iunit) TRunoff%erout(iunit,nt) = temp_erout TRunoff%flow(iunit,nt) = TRunoff%flow(iunit,nt) - TRunoff%erout(iunit,nt) + if (nt==1) then ! if LIQ tracer and there is water + do ntdom=1,nt_rtm_dom ! loop over DOM tracers + temp_eroutdom(ntdom) = temp_eroutdom(ntdom) / TUnit%numDT_r(iunit) + Tdom%domRout(iunit,ntdom) = temp_eroutdom(ntdom) + Tdom%domRoutFlow(iunit,ntdom) = Tdom%domRoutFlow(iunit,ntdom) + Tdom%domRout(iunit,ntdom) + enddo + endif endif end do ! iunit endif ! euler_calc @@ -195,8 +324,11 @@ subroutine Euler ! check for negative channel storage if (negchan < -1.e-10) then write(iulog,*) 'Warning: Negative channel storage found! ',negchan -! call shr_sys_abort('mosart: negative channel storage') + call shr_sys_abort('mosart: negative channel storage') endif + TRunoff%erlateral2= TRunoff%erlateral2/ Tctl%DLevelH2R + Tdom%domRoutFlow = Tdom%domRoutFlow / Tctl%DLevelH2R + Tdom%domToutLat2 = Tdom%domToutLat2 / Tctl%DLevelH2R TRunoff%flow = TRunoff%flow / Tctl%DLevelH2R TRunoff%erout_prev = TRunoff%erout_prev / Tctl%DLevelH2R TRunoff%eroutup_avg = TRunoff%eroutup_avg / Tctl%DLevelH2R diff --git a/src/riverroute/RtmHistFlds.F90 b/src/riverroute/RtmHistFlds.F90 index 8550930..c373317 100644 --- a/src/riverroute/RtmHistFlds.F90 +++ b/src/riverroute/RtmHistFlds.F90 @@ -9,9 +9,9 @@ module RtmHistFlds ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 - use RunoffMod , only : rtmCTL + use RunoffMod , only : rtmCTL, Tdom use RtmHistFile , only : RtmHistAddfld, RtmHistPrintflds - use RtmVar , only : nt_rtm, rtm_tracers + use RtmVar , only : nt_rtm, rtm_tracers, nt_rtm_dom, rtm_tracers_dom implicit none ! @@ -48,7 +48,7 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='RIVER_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART river discharge into ocean: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%runoffocn_nt1, default='inactive') + ptr_rof=rtmCTL%runoffocn_nt1, default='active') call RtmHistAddfld (fname='RIVER_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART river discharge into ocean: '//trim(rtm_tracers(2)), & @@ -72,7 +72,7 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='STORAGE'//'_'//trim(rtm_tracers(1)), units='m3', & avgflag='A', long_name='MOSART storage: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%volr_nt1, default='inactive') + ptr_rof=rtmCTL%volr_nt1, default='active') call RtmHistAddfld (fname='STORAGE'//'_'//trim(rtm_tracers(2)), units='m3', & avgflag='A', long_name='MOSART storage: '//trim(rtm_tracers(2)), & @@ -80,11 +80,11 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='STORAGE_MCH', units='m3', & avgflag='A', long_name='MOSART main channelstorage', & - ptr_rof=rtmCTL%volr_mch, default='inactive') + ptr_rof=rtmCTL%volr_mch, default='active') call RtmHistAddfld (fname='DVOLRDT_LND'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART land change in storage: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%dvolrdtlnd_nt1, default='inactive') + ptr_rof=rtmCTL%dvolrdtlnd_nt1, default='active') call RtmHistAddfld (fname='DVOLRDT_LND'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART land change in storage: '//trim(rtm_tracers(2)), & @@ -92,7 +92,7 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='DVOLRDT_OCN'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART ocean change of storage: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%dvolrdtocn_nt1, default='inactive') + ptr_rof=rtmCTL%dvolrdtocn_nt1, default='active') call RtmHistAddfld (fname='DVOLRDT_OCN'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART ocean change of storage: '//trim(rtm_tracers(2)), & @@ -100,7 +100,7 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='QSUR'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART input surface runoff: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%qsur_nt1, default='inactive') + ptr_rof=rtmCTL%qsur_nt1, default='active') call RtmHistAddfld (fname='QSUR'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART input surface runoff: '//trim(rtm_tracers(2)), & @@ -108,7 +108,7 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='QSUB'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART input subsurface runoff: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%qsub_nt1, default='inactive') + ptr_rof=rtmCTL%qsub_nt1, default='active') call RtmHistAddfld (fname='QSUB'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART input subsurface runoff: '//trim(rtm_tracers(2)), & @@ -116,19 +116,93 @@ subroutine RtmHistFldsInit() call RtmHistAddfld (fname='QGWL'//'_'//trim(rtm_tracers(1)), units='m3/s', & avgflag='A', long_name='MOSART input GWL runoff: '//trim(rtm_tracers(1)), & - ptr_rof=rtmCTL%qgwl_nt1, default='inactive') + ptr_rof=rtmCTL%qgwl_nt1, default='active') call RtmHistAddfld (fname='QGWL'//'_'//trim(rtm_tracers(2)), units='m3/s', & avgflag='A', long_name='MOSART input GWL runoff: '//trim(rtm_tracers(2)), & - ptr_rof=rtmCTL%qgwl_nt2, default='inactive') + ptr_rof=rtmCTL%qgwl_nt2, default='active') call RtmHistAddfld (fname='QIRRIG_FROM_COUPLER', units='m3/s', & avgflag='A', long_name='Amount of water used for irrigation (total flux received from coupler)', & - ptr_rof=rtmCTL%qirrig, default='inactive') + ptr_rof=rtmCTL%qirrig, default='active') call RtmHistAddfld (fname='QIRRIG_ACTUAL', units='m3/s', & avgflag='A', long_name='Actual irrigation (if limited by river storage)', & - ptr_rof=rtmCTL%qirrig_actual, default='inactive') + ptr_rof=rtmCTL%qirrig_actual, default='active') + + call RtmHistAddfld (fname='RIVER_DISCHARGE_OVER_LAND'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART DOM basin flow: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%runofflnddom_ntdom1, default='active') + + call RtmHistAddfld (fname='RIVER_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART DOM discharge into ocean: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%runoffocndom_ntdom1, default='active') + + call RtmHistAddfld (fname='QSUR'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART input surface DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domsur_ntdom1, default='active') + + call RtmHistAddfld (fname='QSUB'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART input subsurface DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domsub_ntdom1, default='active') + + call RtmHistAddfld (fname='STORAGE'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART storage: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%dommas_ntdom1, default='active') + + call RtmHistAddfld (fname='MASS_HILLS'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domH_ntdom1, default='active') + + call RtmHistAddfld (fname='MASS_SUBN'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domT_ntdom1, default='active') + + call RtmHistAddfld (fname='MASS_MAINC'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domR_ntdom1, default='active') + + call RtmHistAddfld (fname='MASS_REST_HILL'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domResthill_ntdom1, default='active') + call RtmHistAddfld (fname='MASS_REST_SUBN'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domRestsubn_ntdom1, default='active') + call RtmHistAddfld (fname='MASS_REST_MAIN'//'_'//trim(rtm_tracers_dom(1)), units='kgC', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domRestmain_ntdom1, default='active') + ! added DOC out of Hills + call RtmHistAddfld (fname='OUT_HILLS'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domHout_ntdom1, default='active') + ! added DOC out of Subn + call RtmHistAddfld (fname='OUT_SUBN'//'_'//trim(rtm_tracers_dom(1)), units='kgC/s', & + avgflag='A', long_name='MOSART DOM: '//trim(rtm_tracers_dom(1)), & + ptr_rof=rtmCTL%domTout_ntdom1, default='active') + ! added water out of Subn + call RtmHistAddfld (fname='OUT_SUBN'//'_'//trim(rtm_tracers(1)), units='m3/s', & + avgflag='A', long_name='MOSART water: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%erlateral2_nt1, default='active') + ! added water out of Hills + call RtmHistAddfld (fname='IN_SUBN'//'_'//trim(rtm_tracers(1)), units='m3/s', & + avgflag='A', long_name='MOSART water: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%etin_nt1, default='active') + + call RtmHistAddfld (fname='OUT_HILLS'//'_'//trim(rtm_tracers(1)), units='m3/s', & + avgflag='A', long_name='MOSART water: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%ehout_nt1, default='active') + + call RtmHistAddfld (fname='MASS_HILLS'//'_'//trim(rtm_tracers(1)), units='m3', & + avgflag='A', long_name='MOSART WATER: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%wh_nt1, default='active') + + call RtmHistAddfld (fname='MASS_SUBN'//'_'//trim(rtm_tracers(1)), units='m3', & + avgflag='A', long_name='MOSART WATER: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%wt_nt1, default='active') + + call RtmHistAddfld (fname='MASS_MAINC'//'_'//trim(rtm_tracers(1)), units='m3', & + avgflag='A', long_name='MOSART WATER: '//trim(rtm_tracers(1)), & + ptr_rof=rtmCTL%wr_nt1, default='active') ! Print masterlist of history fields @@ -180,6 +254,26 @@ subroutine RtmHistFldsSet() rtmCTL%qgwl_nt1(:) = rtmCTL%qgwl(:,1) rtmCTL%qgwl_nt2(:) = rtmCTL%qgwl(:,2) + rtmCTL%domsur_ntdom1(:) = rtmCTL%domsur(:,1) + rtmCTL%domsub_ntdom1(:) = rtmCTL%domsub(:,1) + rtmCTL%dommas_ntdom1(:) = rtmCTL%dommas(:,1) + rtmCTL%runoffocndom_ntdom1(:) = rtmCTL%runoffocndom(:,1) + rtmCTL%runofflnddom_ntdom1(:) = rtmCTL%runofflnddom(:,1) + rtmCTL%domH_ntdom1(:) = rtmCTL%domH(:,1) + rtmCTL%domT_ntdom1(:) = rtmCTL%domT(:,1) + rtmCTL%domR_ntdom1(:) = rtmCTL%domR(:,1) + rtmCTL%domResthill_ntdom1(:) = rtmCTL%domResthill(:,1) + rtmCTL%domRestsubn_ntdom1(:) = rtmCTL%domRestsubn(:,1) + rtmCTL%domRestmain_ntdom1(:) = rtmCTL%domRestmain(:,1) + rtmCTL%domHout_ntdom1(:) = rtmCTL%domHout(:,1) + rtmCTL%domTout_ntdom1(:) = rtmCTL%domTout(:,1) + rtmCTL%wr_nt1(:) = rtmCTL%wr(:,1) + rtmCTL%wt_nt1(:) = rtmCTL%wt(:,1) + rtmCTL%wh_nt1(:) = rtmCTL%wh(:,1) + rtmCTL%etin_nt1(:) = rtmCTL%etin(:,1) + rtmCTL%ehout_nt1(:) = rtmCTL%ehout(:,1) + rtmCTL%erlateral2_nt1(:) = rtmCTL%erlateral2(:,1) + end subroutine RtmHistFldsSet diff --git a/src/riverroute/RtmMod.F90 b/src/riverroute/RtmMod.F90 index 14cf0c6..8925f06 100644 --- a/src/riverroute/RtmMod.F90 +++ b/src/riverroute/RtmMod.F90 @@ -12,7 +12,7 @@ module RtmMod use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_flush use shr_const_mod , only : SHR_CONST_PI, SHR_CONST_CDAY - use RtmVar , only : nt_rtm, rtm_tracers + use RtmVar , only : nt_rtm, rtm_tracers, nt_rtm_dom, rtm_tracers_dom use RtmSpmd , only : masterproc, npes, iam, mpicom_rof, ROFID, mastertask, & MPI_REAL8,MPI_INTEGER,MPI_CHARACTER,MPI_LOGICAL,MPI_MAX use RtmVar , only : re, spval, rtmlon, rtmlat, iulog, ice_runoff, & @@ -24,7 +24,7 @@ module RtmMod barrier_timers use RtmFileUtils , only : getfil, getavu, relavu use RtmTimeManager , only : timemgr_init, get_nstep, get_curr_date - use RtmHistFlds , only : RtmHistFldsInit, RtmHistFldsSet + use RtmHistFlds , only : RtmHistFldsInit, RtmHistFldsSet use RtmHistFile , only : RtmHistUpdateHbuf, RtmHistHtapesWrapup, RtmHistHtapesBuild, & rtmhist_ndens, rtmhist_mfilt, rtmhist_nhtfrq, & rtmhist_avgflag_pertape, rtmhist_avgflag_pertape, & @@ -37,7 +37,9 @@ module RtmMod gsmap_r, & SMatP_dnstrm, avsrc_dnstrm, avdst_dnstrm, & SMatP_direct, avsrc_direct, avdst_direct, & - SMatP_eroutUp, avsrc_eroutUp, avdst_eroutUp + SMatP_eroutUp, avsrc_eroutUp, avdst_eroutUp, & + sMatP_domRUp, avsrc_domRUp, avdst_domRUp, & + Tdom use MOSART_physics_mod, only : Euler use MOSART_physics_mod, only : updatestate_hillslope, updatestate_subnetwork, & updatestate_mainchannel @@ -82,6 +84,12 @@ module RtmMod !local (gdc) real(r8), save, pointer :: evel(:,:) ! effective tracer velocity (m/s) real(r8), save, pointer :: flow(:,:) ! mosart flow (m3/s) + real(r8), save, pointer :: erlateral2(:,:) ! mosart flow (m3/s) + real(r8), save, pointer :: etin(:,:) ! mosart flow (m3/s) + real(r8), save, pointer :: ehout(:,:) ! mosart flow (m3/s) + real(r8), save, pointer :: flowdom(:,:) ! mosart flow (kg/s) + real(r8), save, pointer :: Houtdom(:,:) ! mosart flow (kg/s) + real(r8), save, pointer :: Toutdom(:,:) ! mosart flow (kg/s) real(r8), save, pointer :: erout_prev(:,:) ! erout previous timestep (m3/s) real(r8), save, pointer :: eroutup_avg(:,:)! eroutup average over coupling period (m3/s) real(r8), save, pointer :: erlat_avg(:,:) ! erlateral average over coupling period (m3/s) @@ -136,7 +144,7 @@ subroutine Rtmini(flood_active) real(r8) :: edgee ! East edge of the direction file real(r8) :: edges ! South edge of the direction file real(r8) :: edgew ! West edge of the direction file - integer :: i,j,k,n,ng,g,n2,nt,nn ! loop indices + integer :: i,j,k,n,ng,g,n2,nt,nn,ntdom ! loop indices integer :: i1,j1,i2,j2 integer :: im1,ip1,jm1,jp1,ir,jr,nr ! neighbor indices real(r8) :: deg2rad ! pi/180 @@ -351,6 +359,14 @@ subroutine Rtmini(flood_active) write(iulog,*)'MOSART tracers = ',nt_rtm,trim(rtm_trstr) end if + rtm_trstr = trim(rtm_tracers_dom(1)) + do n = 2,nt_rtm_dom + rtm_trstr = trim(rtm_trstr)//':'//trim(rtm_tracers_dom(n)) + enddo + if (masterproc) then + write(iulog,*)'MOSART tracers dom = ',nt_rtm_dom,trim(rtm_trstr) + end if + !------------------------------------------------------- ! Read input data (river direction file) !------------------------------------------------------- @@ -920,6 +936,12 @@ subroutine Rtmini(flood_active) allocate (evel (rtmCTL%begr:rtmCTL%endr,nt_rtm), & flow (rtmCTL%begr:rtmCTL%endr,nt_rtm), & + flowdom (rtmCTL%begr:rtmCTL%endr,nt_rtm_dom), & + erlateral2 (rtmCTL%begr:rtmCTL%endr,nt_rtm), & + etin (rtmCTL%begr:rtmCTL%endr,nt_rtm), & + ehout (rtmCTL%begr:rtmCTL%endr,nt_rtm), & + Houtdom (rtmCTL%begr:rtmCTL%endr,nt_rtm_dom), & + Toutdom (rtmCTL%begr:rtmCTL%endr,nt_rtm_dom), & erout_prev(rtmCTL%begr:rtmCTL%endr,nt_rtm), & eroutup_avg(rtmCTL%begr:rtmCTL%endr,nt_rtm), & erlat_avg(rtmCTL%begr:rtmCTL%endr,nt_rtm), & @@ -929,6 +951,12 @@ subroutine Rtmini(flood_active) call shr_sys_abort(subname//' Allocationt ERROR flow') end if flow(:,:) = 0._r8 + erlateral2(:,:) = 0._r8 + etin(:,:) = 0._r8 + ehout(:,:) = 0._r8 + flowdom(:,:) = 0._r8 + Houtdom(:,:) = 0._r8 + Toutdom(:,:) = 0._r8 erout_prev(:,:) = 0._r8 eroutup_avg(:,:) = 0._r8 erlat_avg(:,:) = 0._r8 @@ -1317,6 +1345,11 @@ subroutine Rtmini(flood_active) TRunoff%wt = rtmCTL%wt TRunoff%wr = rtmCTL%wr TRunoff%erout= rtmCTL%erout + Tdom%domH = rtmCTL%domH + Tdom%domT = rtmCTL%domT + Tdom%domR = rtmCTL%domR + Tdom%domRout = rtmCTL%domRout + write(iulog,*) 'UPDATED MARIUS' else ! do nt = 1,nt_rtm ! do nr = rtmCTL%begr,rtmCTL%endr @@ -1389,7 +1422,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! ! !LOCAL VARIABLES: !EOP - integer :: i, j, n, nr, ns, nt, n2, nf ! indices + integer :: i, j, n, nr, ns, nt, n2, nf, ntdom! indices real(r8) :: budget_terms(30,nt_rtm) ! BUDGET terms ! BUDGET terms 1-10 are for volumes (m3) ! BUDGET terms 11-30 are for flows (m3/s) @@ -1446,6 +1479,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) budget_terms = 0._r8 flow = 0._r8 + erlateral2 = 0._r8 + etin = 0._r8 + ehout = 0._r8 + flowdom = 0._r8 + Houtdom = 0._r8 + Toutdom = 0._r8 erout_prev = 0._r8 eroutup_avg = 0._r8 erlat_avg = 0._r8 @@ -1491,6 +1530,10 @@ subroutine Rtmrun(rstwr,nlend,rdate) TRunoff%qsub(nr,nt) = rtmCTL%qsub(nr,nt) TRunoff%qgwl(nr,nt) = rtmCTL%qgwl(nr,nt) enddo + do nt = 1,nt_rtm_dom + Tdom%domsur(nr,nt) = rtmCTL%domsur(nr,nt) + Tdom%domsub(nr,nt) = rtmCTL%domsub(nr,nt) + enddo enddo !----------------------------------- @@ -1807,12 +1850,17 @@ subroutine Rtmrun(rstwr,nlend,rdate) call t_stopf('mosartr_budget') ! endif - do nt = 1,nt_rtm + do nr = rtmCTL%begr,rtmCTL%endr + do nt = 1,nt_rtm TRunoff%qsur(nr,nt) = TRunoff%qsur(nr,nt) / rtmCTL%area(nr) TRunoff%qsub(nr,nt) = TRunoff%qsub(nr,nt) / rtmCTL%area(nr) TRunoff%qgwl(nr,nt) = TRunoff%qgwl(nr,nt) / rtmCTL%area(nr) enddo + do nt = 1,nt_rtm_dom + Tdom%domsur(nr,nt) = Tdom%domsur(nr,nt) / rtmCTL%area(nr) + Tdom%domsub(nr,nt) = Tdom%domsub(nr,nt) / rtmCTL%area(nr) + enddo enddo do ns = 1,nsub @@ -1865,9 +1913,19 @@ subroutine Rtmrun(rstwr,nlend,rdate) do nt = 1,nt_rtm do nr = rtmCTL%begr,rtmCTL%endr flow(nr,nt) = flow(nr,nt) + TRunoff%flow(nr,nt) + etin(nr,nt) = etin(nr,nt) + TRunoff%etin(nr,nt) + ehout(nr,nt) = ehout(nr,nt) - TRunoff%ehout(nr,nt)*rtmCTL%area(nr) + erlateral2(nr,nt) = erlateral2(nr,nt) + TRunoff%erlateral2(nr,nt) erout_prev(nr,nt) = erout_prev(nr,nt) + TRunoff%erout_prev(nr,nt) eroutup_avg(nr,nt) = eroutup_avg(nr,nt) + TRunoff%eroutup_avg(nr,nt) erlat_avg(nr,nt) = erlat_avg(nr,nt) + TRunoff%erlat_avg(nr,nt) + if (nt==1) then + do ntdom = 1,nt_rtm_dom + flowdom(nr,ntdom) = flowdom(nr,ntdom) + Tdom%domRoutFlow(nr,ntdom) + Houtdom(nr,ntdom) = Houtdom(nr,ntdom) + Tdom%domHout(nr,ntdom) + Toutdom(nr,ntdom) = Toutdom(nr,ntdom) + Tdom%domToutLat2(nr,ntdom) + enddo + endif enddo enddo @@ -1876,8 +1934,13 @@ subroutine Rtmrun(rstwr,nlend,rdate) !----------------------------------- ! average flow over subcycling !----------------------------------- - + erlateral2 = erlateral2 / float(nsub) + etin = etin / float(nsub) + ehout = ehout / float(nsub) flow = flow / float(nsub) + flowdom = flowdom / float(nsub) + Houtdom = Houtdom / float(nsub) + Toutdom = Toutdom / float(nsub) erout_prev = erout_prev / float(nsub) eroutup_avg = eroutup_avg / float(nsub) erlat_avg = erlat_avg / float(nsub) @@ -1885,28 +1948,62 @@ subroutine Rtmrun(rstwr,nlend,rdate) !----------------------------------- ! update states when subsycling completed !----------------------------------- - - rtmCTL%wh = TRunoff%wh rtmCTL%wt = TRunoff%wt rtmCTL%wr = TRunoff%wr rtmCTL%erout = TRunoff%erout + rtmCTL%domT = Tdom%domT + rtmCTL%domR = Tdom%domR + rtmCTL%domRout = Tdom%domRout + rtmCTL%domResthill = Tdom%domResthill + rtmCTL%domRestsubn = Tdom%domRestsubn + rtmCTL%domRestmain = Tdom%domRestmain do nt = 1,nt_rtm - do nr = rtmCTL%begr,rtmCTL%endr + do nr = rtmCTL%begr,rtmCTL%endr ! both ocean, land and outlet volr_init = rtmCTL%volr(nr,nt) rtmCTL%volr(nr,nt) = (TRunoff%wt(nr,nt) + TRunoff%wr(nr,nt) + & TRunoff%wh(nr,nt)*rtmCTL%area(nr)) + rtmCTL%wh = TRunoff%wh*rtmCTL%area(nr) + if (nt==1) then + do ntdom = 1,nt_rtm_dom + rtmCTL%domH(nr,ntdom) = Tdom%domH(nr,ntdom) * rtmCTL%area(nr) + rtmCTL%dommas(nr,ntdom)=(rtmCTL%area(nr)*Tdom%domH(nr,ntdom) + & + Tdom%domT(nr,ntdom) + & + Tdom%domR(nr,ntdom)) + rtmCTL%domHout(nr,ntdom)=Houtdom(nr,ntdom) + rtmCTL%domTout(nr,ntdom)=Toutdom(nr,ntdom) + enddo + end if rtmCTL%dvolrdt(nr,nt) = (rtmCTL%volr(nr,nt) - volr_init) / delt_coupling rtmCTL%runoff(nr,nt) = flow(nr,nt) - + rtmCTL%runofftot(nr,nt) = rtmCTL%direct(nr,nt) - if (rtmCTL%mask(nr) == 1) then + + rtmCTL%erlateral2(nr,nt)=erlateral2(nr,nt) + rtmCTL%etin(nr,nt)=etin(nr,nt) + rtmCTL%ehout(nr,nt)=ehout(nr,nt) + + if (rtmCTL%mask(nr) == 1) then ! over land rtmCTL%runofflnd(nr,nt) = rtmCTL%runoff(nr,nt) rtmCTL%dvolrdtlnd(nr,nt)= rtmCTL%dvolrdt(nr,nt) - elseif (rtmCTL%mask(nr) >= 2) then + + if (nt==1) then + do ntdom = 1,nt_rtm_dom + rtmCTL%runofflnddom(nr,ntdom)=flowdom(nr,ntdom) + enddo + end if + + elseif (rtmCTL%mask(nr) >= 2) then ! ocean and outlet rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) rtmCTL%runofftot(nr,nt) = rtmCTL%runofftot(nr,nt) + rtmCTL%runoff(nr,nt) rtmCTL%dvolrdtocn(nr,nt)= rtmCTL%dvolrdt(nr,nt) + + if (nt==1) then + do ntdom = 1,nt_rtm_dom + rtmCTL%runoffocndom(nr,ntdom)=flowdom(nr,ntdom) + enddo + end if + endif enddo enddo @@ -2098,7 +2195,7 @@ subroutine RtmFloodInit(frivinp, begr, endr, fthresh, evel ) real(r8) , pointer :: rslope(:) real(r8) , pointer :: max_volr(:) integer, pointer :: compdof(:) ! computational degrees of freedom for pio - integer :: nt,n,cnt ! indices + integer :: nt,ntdom,n,cnt ! indices logical :: readvar ! read variable in or not integer :: ier ! status variable integer :: dids(2) ! variable dimension ids @@ -2195,7 +2292,7 @@ subroutine MOSART_init integer :: dids(2) ! variable dimension ids integer :: dsizes(2) ! variable dimension lengths integer :: ier ! error code - integer :: begr, endr, iunit, nn, n, cnt, nr, nt + integer :: begr, endr, iunit, nn, n, cnt, nr, nt, ntdom integer :: numDT_r, numDT_t integer :: lsize, gsize integer :: igrow, igcol, iwgt @@ -2511,6 +2608,9 @@ subroutine MOSART_init allocate (TRunoff%erlateral(begr:endr,nt_rtm)) TRunoff%erlateral = 0._r8 + allocate (TRunoff%erlateral2(begr:endr,nt_rtm)) + TRunoff%erlateral2 = 0._r8 + allocate (TRunoff%erin(begr:endr,nt_rtm)) TRunoff%erin = 0._r8 @@ -2538,6 +2638,38 @@ subroutine MOSART_init allocate (TPara%c_twid(begr:endr)) TPara%c_twid = 1.0_r8 + !Initialize dom flux variables + allocate (Tdom%domH(begr:endr,nt_rtm_dom)) + Tdom%domH = 0._r8 + allocate (Tdom%domHout(begr:endr,nt_rtm_dom)) + Tdom%domHout = 0._r8 + allocate (Tdom%domT(begr:endr,nt_rtm_dom)) + Tdom%domT = 0._r8 + allocate (Tdom%domTout(begr:endr,nt_rtm_dom)) + Tdom%domTout = 0._r8 + allocate (Tdom%domToutLat(begr:endr,nt_rtm_dom)) + Tdom%domToutLat = 0._r8 + allocate (Tdom%domToutLat2(begr:endr,nt_rtm_dom)) + Tdom%domToutLat2 = 0._r8 + allocate (Tdom%domR(begr:endr,nt_rtm_dom)) + Tdom%domR = 0._r8 + allocate (Tdom%domRout(begr:endr,nt_rtm_dom)) + Tdom%domRout = 0._r8 + allocate (Tdom%domRoutFlow(begr:endr,nt_rtm_dom)) + Tdom%domRoutFlow = 0._r8 + allocate (Tdom%domRUp(begr:endr,nt_rtm_dom)) + Tdom%domRUp = 0._r8 + allocate (Tdom%domResthill(begr:endr,nt_rtm_dom)) + Tdom%domResthill = 0._r8 + allocate (Tdom%domRestsubn(begr:endr,nt_rtm_dom)) + Tdom%domRestsubn = 0._r8 + allocate (Tdom%domRestmain(begr:endr,nt_rtm_dom)) + Tdom%domRestmain = 0._r8 + allocate (Tdom%domsur(begr:endr,nt_rtm_dom)) + Tdom%domsur = 0._r8 + allocate (Tdom%domsub(begr:endr,nt_rtm_dom)) + Tdom%domsub = 0._r8 + call pio_freedecomp(ncid, iodesc_dbl) call pio_freedecomp(ncid, iodesc_int) call pio_closefile(ncid) @@ -2640,6 +2772,7 @@ subroutine MOSART_init enddo call mct_sMatP_Init(sMatP_eroutUp, sMat, gsMap_r, gsMap_r, 0, mpicom_rof, ROFID) + call mct_sMatP_Init(sMatP_domRUp, sMat, gsMap_r, gsMap_r, 0, mpicom_rof, ROFID) elseif (smat_option == 'Xonly' .or. smat_option == 'Yonly') then ! root initialization @@ -2681,6 +2814,7 @@ subroutine MOSART_init call mct_avect_clean(avtmp) call mct_sMatP_Init(sMatP_eroutUp, sMat, gsMap_r, gsMap_r, smat_option, 0, mpicom_rof, ROFID) + call mct_sMatP_Init(sMatP_domRUp, sMat, gsMap_r, gsMap_r, smat_option, 0, mpicom_rof, ROFID) else @@ -2697,8 +2831,11 @@ subroutine MOSART_init if ( masterproc ) write(iulog,*) trim(subname),' MOSART initialize avect ',trim(rList) call mct_aVect_init(avsrc_eroutUp,rList=rList,lsize=rtmCTL%lnumr) call mct_aVect_init(avdst_eroutUp,rList=rList,lsize=rtmCTL%lnumr) + call mct_aVect_init(avsrc_domRUp,rList=rList,lsize=rtmCTL%lnumr) + call mct_aVect_init(avdst_domRUp,rList=rList,lsize=rtmCTL%lnumr) lsize = mct_smat_gNumEl(sMatP_eroutUp%Matrix,mpicom_rof) + lsize = mct_smat_gNumEl(sMatP_domRUp%Matrix,mpicom_rof) if (masterproc) write(iulog,*) subname," Done initializing SmatP_eroutUp, nElements = ",lsize ! keep only sMatP @@ -2732,14 +2869,15 @@ subroutine MOSART_init ! copy avdst to avsrc for next downstream step cnt = 0 call mct_avect_zero(avsrc_eroutUp) + do nr = rtmCTL%begr,rtmCTL%endr cnt = cnt + 1 avsrc_eroutUp%rAttr(1,cnt) = avdst_eroutUp%rAttr(1,cnt) enddo call mct_avect_zero(avdst_eroutUp) - call mct_sMat_avMult(avsrc_eroutUp, sMatP_eroutUp, avdst_eroutUp) + ! add avdst to areatot and compute new global sum cnt = 0 diff --git a/src/riverroute/RtmRestFile.F90 b/src/riverroute/RtmRestFile.F90 index 19c593c..a0c8e58 100644 --- a/src/riverroute/RtmRestFile.F90 +++ b/src/riverroute/RtmRestFile.F90 @@ -17,7 +17,7 @@ module RtmRestFile finidat_rtm, nrevsn_rtm, spval, & nsrContinue, nsrBranch, nsrStartup, & ctitle, version, username, hostname, conventions, source, & - nt_rtm, nt_rtm, rtm_tracers + nt_rtm, nt_rtm, rtm_tracers, nt_rtm_dom use RtmHistFile , only : RtmHistRestart use RtmFileUtils , only : relavu, getavu, opnfil, getfil use RtmTimeManager, only : timemgr_restart, get_nstep, get_curr_date, is_last_step @@ -373,7 +373,7 @@ subroutine RtmRestart(ncid, flag) character(len=*) , intent(in) :: flag ! 'read' or 'write' ! LOCAL VARIABLES: logical :: readvar ! determine if variable is on initial file - integer :: nt,nv,n ! indices + integer :: nt,ntdom,nv,n ! indices real(r8) , pointer :: dfld(:) ! temporary array character(len=32) :: vname,uname character(len=255) :: lname @@ -421,6 +421,7 @@ subroutine RtmRestart(ncid, flag) write(iulog,*) 'Rtm ERROR: illegal nv value a ',nv call shr_sys_abort() endif + if (flag == 'define') then call ncd_defvar(ncid=ncid, varname=trim(vname), & @@ -440,6 +441,54 @@ subroutine RtmRestart(ncid, flag) enddo enddo + + + do nv = 8,11 + do ntdom = 1,nt_rtm_dom + if (nv == 8) then + vname = 'RTM_DOMH_'//trim(rtm_tracers(ntdom)) + lname = 'DOM storage at hillslope in cell' + uname = 'kg' + dfld => rtmCTL%domH(:,ntdom) + elseif (nv == 9) then + vname = 'RTM_DOMT_'//trim(rtm_tracers(ntdom)) + lname = 'DOM storage in tributary channels in cell' + uname = 'kg' + dfld => rtmCTL%domT(:,ntdom) + elseif (nv == 10) then + vname = 'RTM_DOMR_'//trim(rtm_tracers(ntdom)) + lname = 'DOM storage in main channel in cell' + uname = 'kg' + dfld => rtmCTL%domR(:,ntdom) + elseif (nv == 11) then + vname = 'RTM_DOMROUT_'//trim(rtm_tracers(ntdom)) + lname = 'DOM storage in upstream main channels' + uname = 'kg/s' + dfld => rtmCTL%domRout(:,ntdom) + else + write(iulog,*) 'Rtm ERROR: illegal nv value a ',nv + call shr_sys_abort() + endif + + + if (flag == 'define') then + call ncd_defvar(ncid=ncid, varname=trim(vname), & + xtype=ncd_double, dim1name='rtmlon', dim2name='rtmlat', & + long_name=trim(lname), units=trim(uname), fill_value=spval) + else if (flag == 'read' .or. flag == 'write') then + call ncd_io(varname=trim(vname), data=dfld, dim1name='allrof', & + ncid=ncid, flag=flag, readvar=readvar) + if (flag=='read' .and. .not. readvar) then + if (nsrest == nsrContinue) then + call shr_sys_abort() + else + dfld = 0._r8 + end if + end if + end if + + enddo + enddo if (flag == 'read') then do n = rtmCTL%begr,rtmCTL%endr @@ -451,6 +500,14 @@ subroutine RtmRestart(ncid, flag) if (abs(rtmCTL%wt(n,nt)) > 1.e30) rtmCTL%wt(n,nt) = 0. if (abs(rtmCTL%wr(n,nt)) > 1.e30) rtmCTL%wr(n,nt) = 0. if (abs(rtmCTL%erout(n,nt)) > 1.e30) rtmCTL%erout(n,nt) = 0. + if (nt==1) then + do ntdom = 1,nt_rtm_dom + if (abs(rtmCTL%domH(n,ntdom)) > 1.e30) rtmCTL%domH(n,ntdom) = 0. + if (abs(rtmCTL%domT(n,ntdom)) > 1.e30) rtmCTL%domT(n,ntdom) = 0. + if (abs(rtmCTL%domR(n,ntdom)) > 1.e30) rtmCTL%domR(n,ntdom) = 0. + if (abs(rtmCTL%domRout(n,ntdom)) > 1.e30) rtmCTL%domRout(n,ntdom) = 0. + end do + endif end do if (rtmCTL%mask(n) == 1) then do nt = 1,nt_rtm diff --git a/src/riverroute/RtmVar.F90 b/src/riverroute/RtmVar.F90 index 744cf01..b1ec37a 100644 --- a/src/riverroute/RtmVar.F90 +++ b/src/riverroute/RtmVar.F90 @@ -10,6 +10,8 @@ module RtmVar !TODO - nt_rtm and rtm_tracers need to be removed and set by access to the index array integer, parameter, public :: nt_rtm = 2 ! number of tracers character(len=3), parameter, public :: rtm_tracers(nt_rtm) = (/'LIQ','ICE'/) + integer, parameter, public :: nt_rtm_dom = 1 ! number of tracers + character(len=3), parameter, public :: rtm_tracers_dom(nt_rtm_dom) = (/'DOC'/) !DON? ! Constants integer, parameter, private :: iundef = -9999999 diff --git a/src/riverroute/RunoffMod.F90 b/src/riverroute/RunoffMod.F90 index 995be6c..1a5909a 100644 --- a/src/riverroute/RunoffMod.F90 +++ b/src/riverroute/RunoffMod.F90 @@ -10,7 +10,7 @@ module RunoffMod ! ! !USES: use shr_kind_mod, only : r8 => shr_kind_r8 - use RtmVar , only : iulog, spval, nt_rtm + use RtmVar , only : iulog, spval, nt_rtm,nt_rtm_dom use mct_mod ! !PUBLIC TYPES: @@ -30,6 +30,10 @@ module RunoffMod type(mct_sMatP),public :: sMatP_eroutUp ! sparse matrix plus for eroutUp calc type(mct_avect),public :: avsrc_eroutUp ! src avect for SM mult eroutUp calc type(mct_avect),public :: avdst_eroutUp ! dst avect for SM mult eroutUp calc + type(mct_sMatP),public :: sMatP_domRUp ! sparse matrix plus for domRUp calc + type(mct_avect),public :: avsrc_domRUp ! src avect for SM mult domRUp calc + type(mct_avect),public :: avdst_domRUp ! dst avect for SM mult domRUp calc + public :: runoff_flow type runoff_flow @@ -55,13 +59,29 @@ module RunoffMod ! - local real(r8), pointer :: runofflnd(:,:) ! runoff masked for land (m3 H2O/s) real(r8), pointer :: runoffocn(:,:) ! runoff masked for ocn (m3 H2O/s) + real(r8), pointer :: runoffocndom(:,:)! DOM runoff masked for ocn (kgC/s) + real(r8), pointer :: runofflnddom(:,:)! DOM runoff masked for lnd (kgC/s) real(r8), pointer :: runofftot(:,:) ! total runoff masked for ocn (m3 H2O/s) real(r8), pointer :: dvolrdt(:,:) ! RTM change in storage (mm/s) real(r8), pointer :: dvolrdtlnd(:,:) ! dvolrdt masked for land (mm/s) real(r8), pointer :: dvolrdtocn(:,:) ! dvolrdt masked for ocn (mm/s) real(r8), pointer :: volr(:,:) ! RTM storage (m3) + real(r8), pointer :: dommas(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domH(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domT(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domR(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domTout(:,:) ! RTM DOM storage (kgC/s) + real(r8), pointer :: domHout(:,:) ! RTM DOM storage (kgC/s) + real(r8), pointer :: domRoutFlow(:,:) ! RTM DOM storage (kgC/s) + real(r8), pointer :: domRout(:,:) ! RTM DOM storage (kgC/s) + real(r8), pointer :: domResthill(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domRestsubn(:,:) ! RTM DOM storage (kgC) + real(r8), pointer :: domRestmain(:,:) ! RTM DOM storage (kgC) real(r8), pointer :: fthresh(:) ! RTM water flood threshold - + real(r8), pointer :: etin(:,:) + real(r8), pointer :: ehout(:,:) + real(r8), pointer :: erlateral2(:,:) + ! - restarts real(r8), pointer :: wh(:,:) ! MOSART hillslope surface water storage (m) real(r8), pointer :: wt(:,:) ! MOSART sub-network water storage (m3) @@ -72,6 +92,8 @@ module RunoffMod real(r8), pointer :: qsur(:,:) ! coupler surface forcing [m3/s] real(r8), pointer :: qsub(:,:) ! coupler subsurface forcing [m3/s] real(r8), pointer :: qgwl(:,:) ! coupler glacier/wetland/lake forcing [m3/s] + real(r8), pointer :: domsur(:,:) ! surface dom masked for land (kgC/s) + real(r8), pointer :: domsub(:,:) ! subsurface dom masked for land (kgC/s) ! - outputs real(r8), pointer :: flood(:) ! coupler return flood water sent back to clm [m3/s] @@ -103,6 +125,25 @@ module RunoffMod real(r8), pointer :: qsub_nt2(:) real(r8), pointer :: qgwl_nt1(:) real(r8), pointer :: qgwl_nt2(:) + real(r8), pointer :: domResthill_ntdom1(:) + real(r8), pointer :: domRestsubn_ntdom1(:) + real(r8), pointer :: domRestmain_ntdom1(:) + real(r8), pointer :: domsur_ntdom1(:) + real(r8), pointer :: domsub_ntdom1(:) + real(r8), pointer :: dommas_ntdom1(:) + real(r8), pointer :: runoffocndom_ntdom1(:) + real(r8), pointer :: runofflnddom_ntdom1(:) + real(r8), pointer :: domH_ntdom1(:) + real(r8), pointer :: domT_ntdom1(:) + real(r8), pointer :: domR_ntdom1(:) + real(r8), pointer :: domHout_ntdom1(:) + real(r8), pointer :: domTout_ntdom1(:) + real(r8), pointer :: wh_nt1(:) + real(r8), pointer :: wt_nt1(:) + real(r8), pointer :: wr_nt1(:) + real(r8), pointer :: etin_nt1(:) + real(r8), pointer :: ehout_nt1(:) + real(r8), pointer :: erlateral2_nt1(:) end type runoff_flow @@ -242,6 +283,7 @@ module RunoffMod !! exchange fluxes real(r8), pointer :: erlg(:,:) ! evaporation, [m/s] real(r8), pointer :: erlateral(:,:) ! lateral flow from hillslope, including surface and subsurface runoff generation components, [m3/s] + real(r8), pointer :: erlateral2(:,:) ! lateral flow from hillslope, including surface and subsurface runoff generation components, [m3/s] real(r8), pointer :: erin(:,:) ! inflow from upstream links, [m3/s] real(r8), pointer :: erout(:,:) ! outflow into downstream links, [m3/s] real(r8), pointer :: erout_prev(:,:) ! outflow into downstream links from previous timestep, [m3/s] @@ -270,13 +312,37 @@ module RunoffMod real(r8), pointer :: c_nr(:) ! coefficient to adjust the manning's roughness of channels real(r8), pointer :: c_nh(:) ! coefficient to adjust the manning's roughness of overland flow across hillslopes real(r8), pointer :: c_twid(:) ! coefficient to adjust the width of sub-reach channel - end type Tparameter - + end type Tparameter + + ! DOM status and flux variables + public :: Domflux + type Domflux + real(r8), pointer :: domsur(:,:) ! surface DOM flow from land (kgC/s) + real(r8), pointer :: domsub(:,:) ! subsurface DOM flow from land (kgC/s) + !hillslope + real(r8), pointer :: domH(:,:) ! dissolved organic matter in hillslope (kgC) + real(r8), pointer :: domHout(:,:) ! dissolved organic matter generated from hillslope (kgC/s) + !sub-network + real(r8), pointer :: domT(:,:) ! dom mass in sub-network (kgC) + real(r8), pointer :: domTout(:,:) ! dom discharge from sub-network into main reach (kgC/s) + real(r8), pointer :: domToutLat(:,:) + real(r8), pointer :: domToutLat2(:,:) + !main channel upstream interactions + real(r8), pointer :: domR(:,:) ! dom mass in main channel (kgC) + real(r8), pointer :: domRout(:,:) ! dom discharge from main channel into downstream gridcells (kgC/s) + real(r8), pointer :: domRoutFlow(:,:) + real(r8), pointer :: domRUp(:,:) ! outflow sum of upstream gridcells (kgC/m3) + real(r8), pointer :: domResthill(:,:) ! excess DOM in mosart (kgC) + real(r8), pointer :: domRestsubn(:,:) ! excess DOM in mosart (kgC) + real(r8), pointer :: domRestmain(:,:) ! excess DOM in mosart (kgC) + end type Domflux + !== Hongyi type (Tcontrol) , public :: Tctl type (Tspatialunit), public :: TUnit type (TstatusFlux) , public :: TRunoff type (Tparameter) , public :: TPara + type (Domflux) , public :: Tdom !== Hongyi type (runoff_flow) , public :: rtmCTL @@ -333,12 +399,48 @@ subroutine RunoffInit(begr, endr, numr) rtmCTL%wh(begr:endr,nt_rtm), & rtmCTL%wt(begr:endr,nt_rtm), & rtmCTL%wr(begr:endr,nt_rtm), & + rtmCTL%wh_nt1(begr:endr), & + rtmCTL%wt_nt1(begr:endr), & + rtmCTL%wr_nt1(begr:endr), & rtmCTL%erout(begr:endr,nt_rtm), & rtmCTL%qsur(begr:endr,nt_rtm), & rtmCTL%qsub(begr:endr,nt_rtm), & rtmCTL%qgwl(begr:endr,nt_rtm), & rtmCTL%qirrig(begr:endr), & rtmCTL%qirrig_actual(begr:endr), & + rtmCTL%runofflnddom(begr:endr,nt_rtm_dom), & + rtmCTL%runoffocndom(begr:endr,nt_rtm_dom), & + rtmCTL%domsur(begr:endr,nt_rtm_dom), & + rtmCTL%domsub(begr:endr,nt_rtm_dom), & + rtmCTL%dommas(begr:endr,nt_rtm_dom), & + rtmCTL%runofflnddom_ntdom1(begr:endr), & + rtmCTL%runoffocndom_ntdom1(begr:endr), & + rtmCTL%domsur_ntdom1(begr:endr), & + rtmCTL%domsub_ntdom1(begr:endr), & + rtmCTL%dommas_ntdom1(begr:endr), & + rtmCTL%domResthill_ntdom1(begr:endr), & + rtmCTL%domResthill(begr:endr,nt_rtm_dom), & + rtmCTL%domRestsubn_ntdom1(begr:endr), & + rtmCTL%domRestsubn(begr:endr,nt_rtm_dom), & + rtmCTL%domRestmain_ntdom1(begr:endr), & + rtmCTL%domRestmain(begr:endr,nt_rtm_dom), & + rtmCTL%domH_ntdom1(begr:endr), & + rtmCTL%domH(begr:endr,nt_rtm_dom), & + rtmCTL%domT_ntdom1(begr:endr), & + rtmCTL%domT(begr:endr,nt_rtm_dom), & + rtmCTL%domHout_ntdom1(begr:endr), & + rtmCTL%domHout(begr:endr,nt_rtm_dom), & + rtmCTL%domTout_ntdom1(begr:endr), & + rtmCTL%domTout(begr:endr,nt_rtm_dom), & + rtmCTL%domR_ntdom1(begr:endr), & + rtmCTL%domR(begr:endr,nt_rtm_dom), & + rtmCTL%domRout(begr:endr,nt_rtm_dom), & + rtmCTL%erlateral2(begr:endr,nt_rtm), & + rtmCTL%erlateral2_nt1(begr:endr), & + rtmCTL%etin(begr:endr,nt_rtm), & + rtmCTL%etin_nt1(begr:endr), & + rtmCTL%ehout(begr:endr,nt_rtm), & + rtmCTL%ehout_nt1(begr:endr), & stat=ier) if (ier /= 0) then write(iulog,*)'Rtmini ERROR allocation of runoff local arrays' @@ -362,6 +464,17 @@ subroutine RunoffInit(begr, endr, numr) rtmCTL%qsur(:,:) = 0._r8 rtmCTL%qsub(:,:) = 0._r8 rtmCTL%qgwl(:,:) = 0._r8 + rtmCTL%erlateral2(:,:) =0._r8 + rtmCTL%etin(:,:) =0._r8 + rtmCTL%ehout(:,:) =0._r8 + + rtmCTL%runofflnddom(:,:)=spval + rtmCTL%runoffocndom(:,:)=spval + rtmCTL%domsur(:,:) =0._r8 + rtmCTL%domsub(:,:) =0._r8 + rtmCTL%domTout(:,:) =0._r8 + rtmCTL%domHout(:,:) =0._r8 + end subroutine RunoffInit