diff --git a/models/rof/mosart/bld/namelist_files/namelist_defaults_mosart.xml b/models/rof/mosart/bld/namelist_files/namelist_defaults_mosart.xml
index 558ea97aa92f..a521cdfe2676 100644
--- a/models/rof/mosart/bld/namelist_files/namelist_defaults_mosart.xml
+++ b/models/rof/mosart/bld/namelist_files/namelist_defaults_mosart.xml
@@ -28,7 +28,7 @@ for the CLM data in the CESM distribution
-rof/mosart/MOSART_Global_half_20151015a.nc
+rof/mosart/MOSART_Global_half_20151130a.nc
rof/mosart/MOSART_NLDAS_8th.nc
diff --git a/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90 b/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
index a23d690058ad..465b64e0e542 100644
--- a/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
+++ b/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
@@ -54,18 +54,20 @@ subroutine Euler
! hillslope
!------------------
- call t_startf('mosart_hillslope')
+ 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 .and. TUnit%areaTotal(iunit) > 0._r8) 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)
endif
end do
+ endif
end do
- call t_stopf('mosart_hillslope')
+ call t_stopf('mosartr_hillslope')
TRunoff%flow = 0._r8
TRunoff%erout_prev = 0._r8
@@ -76,20 +78,23 @@ subroutine Euler
!--- accumulate/average erout at prior timestep (used in eroutUp calc) for budget analysis
do nt=1,nt_rtm
+ if (TUnit%euler_calc(nt)) then
do iunit=rtmCTL%begr,rtmCTL%endr
TRunoff%erout_prev(iunit,nt) = TRunoff%erout_prev(iunit,nt) + TRunoff%erout(iunit,nt)
- enddo
- enddo
+ end do
+ endif
+ end do
!------------------
! subnetwork
!------------------
- call t_startf('mosart_subnetwork')
+ call t_startf('mosartr_subnetwork')
TRunoff%erlateral(:,:) = 0._r8
do nt=1,nt_rtm
+ if (TUnit%euler_calc(nt)) then
do iunit=rtmCTL%begr,rtmCTL%endr
- if(TUnit%mask(iunit) > 0 .and. TUnit%areaTotal(iunit) > 0._r8) then
+ if(TUnit%mask(iunit) > 0) then
localDeltaT = Tctl%DeltaT/Tctl%DLevelH2R/TUnit%numDT_t(iunit)
do k=1,TUnit%numDT_t(iunit)
call subnetworkRouting(iunit,nt,localDeltaT)
@@ -100,20 +105,21 @@ subroutine Euler
TRunoff%erlateral(iunit,nt) = TRunoff%erlateral(iunit,nt) / TUnit%numDT_t(iunit)
endif
end do ! iunit
+ endif ! euler_calc
end do ! nt
- call t_stopf('mosart_subnetwork')
+ call t_stopf('mosartr_subnetwork')
!------------------
! upstream interactions
!------------------
if (barrier_timers) then
- call t_startf('mosart_SMeroutUp_barrier')
+ call t_startf('mosartr_SMeroutUp_barrier')
call mpi_barrier(mpicom_rof,ier)
- call t_stopf('mosart_SMeroutUp_barrier')
+ call t_stopf('mosartr_SMeroutUp_barrier')
endif
- call t_startf('mosart_SMeroutUp')
+ call t_startf('mosartr_SMeroutUp')
TRunoff%eroutUp = 0._r8
#ifdef NO_MCT
do iunit=rtmCTL%begr,rtmCTL%endr
@@ -147,7 +153,7 @@ subroutine Euler
enddo
enddo
#endif
- call t_stopf('mosart_SMeroutUp')
+ call t_stopf('mosartr_SMeroutUp')
TRunoff%eroutup_avg = TRunoff%eroutup_avg + TRunoff%eroutUp
TRunoff%erlat_avg = TRunoff%erlat_avg + TRunoff%erlateral
@@ -156,10 +162,11 @@ subroutine Euler
! channel routing
!------------------
- call t_startf('mosart_chanroute')
+ call t_startf('mosartr_chanroute')
do nt=1,nt_rtm
+ if (TUnit%euler_calc(nt)) then
do iunit=rtmCTL%begr,rtmCTL%endr
- if(TUnit%mask(iunit) > 0 .and. TUnit%areaTotal(iunit) > TINYVALUE) then
+ if(TUnit%mask(iunit) > 0) then
localDeltaT = Tctl%DeltaT/Tctl%DLevelH2R/TUnit%numDT_r(iunit)
temp_erout = 0._r8
do k=1,TUnit%numDT_r(iunit)
@@ -177,11 +184,12 @@ subroutine Euler
TRunoff%erout(iunit,nt) = temp_erout
TRunoff%flow(iunit,nt) = TRunoff%flow(iunit,nt) - TRunoff%erout(iunit,nt)
endif
- end do
- end do
+ end do ! iunit
+ endif ! euler_calc
+ end do ! nt
negchan = min(negchan, minval(TRunoff%wr(:,:)))
- call t_stopf('mosart_chanroute')
+ call t_stopf('mosartr_chanroute')
end do
! check for negative channel storage
@@ -294,7 +302,7 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
TRunoff%vr(iunit,nt) = 0._r8
TRunoff%erout(iunit,nt) = -TRunoff%erin(iunit,nt)-TRunoff%erlateral(iunit,nt)
else
- if(TUnit%areaTotal(iunit)/TUnit%rwidth(iunit)/TUnit%rlen(iunit) > 1e6_r8) then
+ if(TUnit%areaTotal2(iunit)/TUnit%rwidth(iunit)/TUnit%rlen(iunit) > 1e6_r8) then
TRunoff%erout(iunit,nt) = -TRunoff%erin(iunit,nt)-TRunoff%erlateral(iunit,nt)
else
! !TRunoff%vr(iunit,nt) = CRVRMAN(TUnit%rslp(iunit), TUnit%nr(iunit), TRunoff%rr(iunit,nt))
diff --git a/models/rof/mosart/src/riverroute/RtmHistFile.F90 b/models/rof/mosart/src/riverroute/RtmHistFile.F90
index a84b455508b6..79e731e6bb15 100644
--- a/models/rof/mosart/src/riverroute/RtmHistFile.F90
+++ b/models/rof/mosart/src/riverroute/RtmHistFile.F90
@@ -8,7 +8,7 @@ module RtmHistFile
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_sys_mod , only : shr_sys_flush, shr_sys_abort
- use RunoffMod , only : rtmCTL
+ use RunoffMod , only : rtmCTL, Tunit
use RtmVar , only : rtmlon, rtmlat, spval, ispval, secspday, frivinp_rtm, &
iulog, nsrest, caseid, inst_suffix, nsrStartup, nsrBranch, &
ctitle, version, hostname, username, conventions, source
@@ -809,8 +809,14 @@ subroutine htape_timeconst(t, mode)
long_name='runoff coordinate longitude', units='degrees_east', ncid=nfid(t))
call ncd_defvar(varname='lat', xtype=tape(t)%ncprec, dim1name='lat', &
long_name='runoff coordinate latitude', units='degrees_north', ncid=nfid(t))
-! call ncd_defvar(varname='mask', xtype=ncd_int, dim1name='lon', dim2name='lat', &
-! long_name='runoff mask', units='unitless', ncid=nfid(t))
+ call ncd_defvar(varname='mask', xtype=ncd_int, dim1name='lon', dim2name='lat', &
+ long_name='runoff mask', units='unitless', ncid=nfid(t))
+ call ncd_defvar(varname='area', xtype=tape(t)%ncprec, dim1name='lon', dim2name='lat', &
+ long_name='runoff grid area', units='m2', ncid=nfid(t))
+ call ncd_defvar(varname='areatotal', xtype=tape(t)%ncprec, dim1name='lon', dim2name='lat', &
+ long_name='basin upstream areatotal', units='m2', ncid=nfid(t))
+ call ncd_defvar(varname='areatotal2', xtype=tape(t)%ncprec, dim1name='lon', dim2name='lat', &
+ long_name='computed basin upstream areatotal', units='m2', ncid=nfid(t))
else if (mode == 'write') then
@@ -839,7 +845,14 @@ subroutine htape_timeconst(t, mode)
call ncd_io(varname='lon', data=rtmCTL%rlon, ncid=nfid(t), flag='write')
call ncd_io(varname='lat', data=rtmCTL%rlat, ncid=nfid(t), flag='write')
-! call ncd_io(varname='mask', data=rtmCTL%mask, ncid=nfid(t), flag='write')
+ call ncd_io(flag='write', varname='mask', dim1name='allrof', &
+ data=rtmCTL%mask, ncid=nfid(t))
+ call ncd_io(flag='write', varname='area', dim1name='allrof', &
+ data=rtmCTL%area, ncid=nfid(t))
+ call ncd_io(flag='write', varname='areatotal', dim1name='allrof', &
+ data=Tunit%areatotal, ncid=nfid(t))
+ call ncd_io(flag='write', varname='areatotal2', dim1name='allrof', &
+ data=Tunit%areatotal2, ncid=nfid(t))
endif
diff --git a/models/rof/mosart/src/riverroute/RtmHistFlds.F90 b/models/rof/mosart/src/riverroute/RtmHistFlds.F90
index 8867e66f2df5..8f770367f61c 100644
--- a/models/rof/mosart/src/riverroute/RtmHistFlds.F90
+++ b/models/rof/mosart/src/riverroute/RtmHistFlds.F90
@@ -38,27 +38,43 @@ subroutine RtmHistFldsInit()
implicit none
!-------------------------------------------------------
- call RtmHistAddfld (fname='QCHANR'//'_'//trim(rtm_tracers(1)), units='m3/s', &
- avgflag='A', long_name='MOSART river flow: '//trim(rtm_tracers(1)), &
+ call RtmHistAddfld (fname='RIVER_DISCHARGE_OVER_LAND'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ avgflag='A', long_name='MOSART river basin flow: '//trim(rtm_tracers(1)), &
ptr_rof=rtmCTL%runofflnd_nt1, default='active')
- call RtmHistAddfld (fname='QCHANR'//'_'//trim(rtm_tracers(2)), units='m3/s', &
- avgflag='A', long_name='MOSART river flow: '//trim(rtm_tracers(2)), &
+ call RtmHistAddfld (fname='RIVER_DISCHARGE_OVER_LAND'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART river basin flow: '//trim(rtm_tracers(2)), &
ptr_rof=rtmCTL%runofflnd_nt2, default='active')
- call RtmHistAddfld (fname='QCHOCNR'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ 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='active')
- call RtmHistAddfld (fname='QCHOCNR'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ 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)), &
ptr_rof=rtmCTL%runoffocn_nt2, default='active')
- call RtmHistAddfld (fname='VOLR'//'_'//trim(rtm_tracers(1)), units='m3', &
+ call RtmHistAddfld (fname='TOTAL_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ avgflag='A', long_name='MOSART total discharge into ocean: '//trim(rtm_tracers(1)), &
+ ptr_rof=rtmCTL%runofftot_nt1, default='active')
+
+ call RtmHistAddfld (fname='TOTAL_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART total discharge into ocean: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%runofftot_nt2, default='active')
+
+ call RtmHistAddfld (fname='DIRECT_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ avgflag='A', long_name='MOSART direct discharge into ocean: '//trim(rtm_tracers(1)), &
+ ptr_rof=rtmCTL%runoffdir_nt1, default='active')
+
+ call RtmHistAddfld (fname='DIRECT_DISCHARGE_TO_OCEAN'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART direct discharge into ocean: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%runoffdir_nt2, default='active')
+
+ 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='active')
- call RtmHistAddfld (fname='VOLR'//'_'//trim(rtm_tracers(2)), units='m3', &
+ call RtmHistAddfld (fname='STORAGE'//'_'//trim(rtm_tracers(2)), units='m3', &
avgflag='A', long_name='MOSART storage: '//trim(rtm_tracers(2)), &
ptr_rof=rtmCTL%volr_nt2, default='active')
@@ -135,6 +151,12 @@ subroutine RtmHistFldsSet()
rtmCTL%runoffocn_nt1(:) = rtmCTL%runoffocn(:,1)
rtmCTL%runoffocn_nt2(:) = rtmCTL%runoffocn(:,2)
+ rtmCTL%runofftot_nt1(:) = rtmCTL%runofftot(:,1)
+ rtmCTL%runofftot_nt2(:) = rtmCTL%runofftot(:,2)
+
+ rtmCTL%runoffdir_nt1(:) = rtmCTL%direct(:,1)
+ rtmCTL%runoffdir_nt2(:) = rtmCTL%direct(:,2)
+
rtmCTL%dvolrdtlnd_nt1(:) = rtmCTL%dvolrdtlnd(:,1)
rtmCTL%dvolrdtlnd_nt2(:) = rtmCTL%dvolrdtlnd(:,2)
diff --git a/models/rof/mosart/src/riverroute/RtmMod.F90 b/models/rof/mosart/src/riverroute/RtmMod.F90
index e0dea28ce7ee..a761d411da1c 100644
--- a/models/rof/mosart/src/riverroute/RtmMod.F90
+++ b/models/rof/mosart/src/riverroute/RtmMod.F90
@@ -39,11 +39,11 @@ module RtmMod
use MOSART_physics_mod, only : Euler
use MOSART_physics_mod, only : updatestate_hillslope, updatestate_subnetwork, &
updatestate_mainchannel
- use RtmIO
#ifdef INCLUDE_WRM
use WRM_subw_IO_mod , only : WRM_init
use WRM_modules , only : Euler_WRM
#endif
+ use RtmIO
use mct_mod
use perf_mod
use pio
@@ -1016,18 +1016,6 @@ subroutine Rtmini(rtm_active,flood_active)
rtmCTL%lonc(nr) = rtmCTL%rlon(i)
rtmCTL%latc(nr) = rtmCTL%rlat(j)
- if (rtmCTL%mask(nr) == 1) then
- do nt = 1,nt_rtm
- rtmCTL%runofflnd(nr,nt) = rtmCTL%runoff(nr,nt)
- rtmCTL%dvolrdtlnd(nr,nt)= rtmCTL%dvolrdt(nr,nt)
- enddo
- elseif (rtmCTL%mask(nr) >= 2) then
- do nt = 1,nt_rtm
- rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt)
- rtmCTL%dvolrdtocn(nr,nt)= rtmCTL%dvolrdt(nr,nt)
- enddo
- endif
-
rtmCTL%outletg(nr) = idxocn(n)
rtmCTL%area(nr) = area_global(n)
lrtmarea = lrtmarea + rtmCTL%area(nr)
@@ -1176,10 +1164,7 @@ subroutine Rtmini(rtm_active,flood_active)
! would be in the full matrix. Nrows= size of output vector=nb.
! Ncols = size of input vector = na.
- cnt = 0
- do nr=rtmCTL%begr,rtmCTL%endr
- if(rtmCTL%outletg(nr) > 0) cnt = cnt + 1
- enddo
+ cnt = rtmCTL%endr - rtmCTL%begr + 1
call mct_sMat_init(sMat, gsize, gsize, cnt)
igrow = mct_sMat_indexIA(sMat,'grow')
@@ -1192,8 +1177,17 @@ subroutine Rtmini(rtm_active,flood_active)
sMat%data%rAttr(iwgt ,cnt) = 1.0_r8
sMat%data%iAttr(igrow,cnt) = rtmCTL%outletg(nr)
sMat%data%iAttr(igcol,cnt) = rtmCTL%gindex(nr)
+ else
+ cnt = cnt + 1
+ sMat%data%rAttr(iwgt ,cnt) = 1.0_r8
+ sMat%data%iAttr(igrow,cnt) = rtmCTL%gindex(nr)
+ sMat%data%iAttr(igcol,cnt) = rtmCTL%gindex(nr)
endif
enddo
+ if (cnt /= rtmCTL%endr - rtmCTL%begr + 1) then
+ write(iulog,*) trim(subname),' MOSART ERROR: smat cnt1 ',cnt,rtmCTL%endr-rtmCTL%begr+1
+ call shr_sys_abort(trim(subname)//' ERROR smat cnt1')
+ endif
call mct_sMatP_Init(sMatP_direct, sMat, gsMap_r, gsMap_r, 0, mpicom_rof, ROFID)
@@ -1211,12 +1205,8 @@ subroutine Rtmini(rtm_active,flood_active)
enddo
call mct_avect_gather(avtmp,avtmpG,gsmap_r,mastertask,mpicom_rof)
if (masterproc) then
- cnt = 0
- do n = 1,rtmlon*rtmlat
- if (avtmpG%rAttr(2,n) > 0) then
- cnt = cnt + 1
- endif
- enddo
+
+ cnt = rtmlon*rtmlat
call mct_sMat_init(sMat, gsize, gsize, cnt)
igrow = mct_sMat_indexIA(sMat,'grow')
@@ -1230,8 +1220,17 @@ subroutine Rtmini(rtm_active,flood_active)
sMat%data%rAttr(iwgt ,cnt) = 1.0_r8
sMat%data%iAttr(igrow,cnt) = avtmpG%rAttr(2,n)
sMat%data%iAttr(igcol,cnt) = avtmpG%rAttr(1,n)
+ else
+ cnt = cnt + 1
+ sMat%data%rAttr(iwgt ,cnt) = 1.0_r8
+ sMat%data%iAttr(igrow,cnt) = avtmpG%rAttr(1,n)
+ sMat%data%iAttr(igcol,cnt) = avtmpG%rAttr(1,n)
endif
enddo
+ if (cnt /= rtmlon*rtmlat) then
+ write(iulog,*) trim(subname),' MOSART ERROR: smat cnt2 ',cnt,rtmlon*rtmlat
+ call shr_sys_abort(trim(subname)//' ERROR smat cnt2')
+ endif
call mct_avect_clean(avtmpG)
else
call mct_sMat_init(sMat,1,1,1)
@@ -1415,12 +1414,14 @@ subroutine Rtmrun(rstwr,nlend,rdate)
! !LOCAL VARIABLES:
!EOP
integer :: i, j, n, nr, ns, nt, n2, nf ! indices
- real(r8) :: budget_terms(30,nt_rtm) ! BUDGET terms
+ 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)
real(r8) :: budget_input, budget_output, budget_volume, budget_total, &
- budget_euler, budget_eroutlag ! BUDGET terms
- real(r8) :: budget_global(30,nt_rtm) ! global budget sum
+ budget_euler, budget_eroutlag
+ real(r8),save :: budget_accum(nt_rtm) ! BUDGET accumulator over run
+ integer ,save :: budget_accum_cnt ! counter for budget_accum
+ real(r8) :: budget_global(30,nt_rtm) ! global budget sum
logical :: budget_check ! do global budget check
real(r8) :: volr_init ! temporary storage to compute dvolrdt
real(r8),parameter :: budget_tolerance = 1.0e-6 ! budget tolerance, m3/day
@@ -1457,6 +1458,8 @@ subroutine Rtmrun(rstwr,nlend,rdate)
delt_coupling = coupling_period*1.0_r8
if (first_call) then
+ budget_accum = 0._r8
+ budget_accum_cnt = 0
if (masterproc) write(iulog,'(2a,g20.12)') trim(subname),' MOSART coupling period ',delt_coupling
end if
@@ -1481,7 +1484,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
! BUDGET
! BUDGET terms 1-10 are for volumes (m3)
! BUDGET terms 11-30 are for flows (m3/s)
- if (budget_check) then
+! if (budget_check) then
call t_startf('mosartr_budget')
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1498,7 +1501,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
enddo
enddo
call t_stopf('mosartr_budget')
- endif
+! endif
! data for euler solver, in m3/s here
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1563,8 +1566,78 @@ subroutine Rtmrun(rstwr,nlend,rdate)
do nr = rtmCTL%begr,rtmCTL%endr
cnt = cnt + 1
do nt = 1,nt_rtm
- avsrc_direct%rAttr(nt,cnt) = TRunoff%qgwl(nr,nt)
+
+ !-------------------------------
+ ! This water is routed directly to the outlet and passed out
+ ! Turn on and off terms as desired via commenting them out
+ !-------------------------------
+
+ !---- all dto water, none was going to TRunoff ---
+ !---- *** DO NOT TURN THIS ONE OFF, conservation will fail *** ---
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + rtmCTL%qdto(nr,nt)
+
+ !---- negative gwl water less than channel volume, remove from TRunoff ---
+ !---- scs
+ qgwl_volume = TRunoff%qgwl(nr,nt) * delt_mosart
+ river_volume_minimum = river_depth_minimum * rtmCTL%area(nr)
+ !---- if qgwl is negative, and adding it to the main channel would bring
+ !---- main channel storage below a threshold, send qgwl directly to ocean
+ if (((qgwl_volume + TRunoff%wr(nr,nt)) < river_volume_minimum) &
+ .and. (TRunoff%qgwl(nr,nt) < 0._r8)) then
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + TRunoff%qgwl(nr,nt)
+ TRunoff%qgwl(nr,nt) = 0._r8
+ endif
+ !---- scs
+
+ !---- negative qgwl water, remove from TRunoff ---
+ if (TRunoff%qgwl(nr,nt) < 0._r8) then
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + TRunoff%qgwl(nr,nt)
+ TRunoff%qgwl(nr,nt) = 0._r8
+ endif
+
+ !---- all gwl water, remove from TRunoff ---
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + TRunoff%qgwl(nr,nt)
TRunoff%qgwl(nr,nt) = 0._r8
+
+ !---- negative qsub water, remove from TRunoff ---
+ if (TRunoff%qsub(nr,nt) < 0._r8) then
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + TRunoff%qsub(nr,nt)
+ TRunoff%qsub(nr,nt) = 0._r8
+ endif
+
+ !---- negative qsur water, remove from TRunoff ---
+ if (TRunoff%qsur(nr,nt) < 0._r8) then
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + TRunoff%qsur(nr,nt)
+ TRunoff%qsur(nr,nt) = 0._r8
+ endif
+
+ !---- water outside the basin ---
+ !---- *** DO NOT TURN THIS ONE OFF, conservation will fail *** ---
+ if (TUnit%mask(nr) > 0) then
+ ! mosart euler
+ else
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + &
+ TRunoff%qsub(nr,nt) + &
+ TRunoff%qsur(nr,nt) + &
+ TRunoff%qgwl(nr,nt)
+ TRunoff%qsub(nr,nt) = 0._r8
+ TRunoff%qsur(nr,nt) = 0._r8
+ TRunoff%qgwl(nr,nt) = 0._r8
+ endif
+
+ !---- all nt=2 water ---
+ !---- can turn off euler_calc for this tracer ----
+ if (nt == 2) then
+ TUnit%euler_calc(nt) = .false.
+ avsrc_direct%rAttr(nt,cnt) = avsrc_direct%rAttr(nt,cnt) + &
+ TRunoff%qsub(nr,nt) + &
+ TRunoff%qsur(nr,nt) + &
+ TRunoff%qgwl(nr,nt)
+ TRunoff%qsub(nr,nt) = 0._r8
+ TRunoff%qsur(nr,nt) = 0._r8
+ TRunoff%qgwl(nr,nt) = 0._r8
+ endif
+
enddo
enddo
call mct_avect_zero(avdst_direct)
@@ -1576,70 +1649,23 @@ subroutine Rtmrun(rstwr,nlend,rdate)
do nr = rtmCTL%begr,rtmCTL%endr
cnt = cnt + 1
do nt = 1,nt_rtm
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + avdst_direct%rAttr(nt,cnt)
+ rtmCTL%direct(nr,nt) = avdst_direct%rAttr(nt,cnt)
enddo
enddo
call t_stopf('mosartr_SMdirect')
- !-----------------------------------
- !--- add other direct terms
- !-----------------------------------
-
- do nt = 1,nt_rtm
- do nr = rtmCTL%begr,rtmCTL%endr
- ! --- dto water
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + rtmCTL%qdto(nr,nt)
-
- ! --- gwl negative volume
-!scs
- qgwl_volume = TRunoff%qgwl(nr,nt) * delt_mosart
- river_volume_minimum = river_depth_minimum * rtmCTL%area(nr)
-! if qgwl is negative, and adding it to the main channel would bring
-! main channel storage below a threshold, send qgwl directly to ocean
- if (((qgwl_volume + TRunoff%wr(nr,nt)) < river_volume_minimum) &
- .and. (TRunoff%qgwl(nr,nt) < 0._r8)) then
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + TRunoff%qgwl(nr,nt)
- TRunoff%qgwl(nr,nt) = 0._r8
- endif
-!scs
-
- if (TRunoff%qgwl(nr,nt) < 0._r8) then
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + TRunoff%qgwl(nr,nt)
- TRunoff%qgwl(nr,nt) = 0._r8
- endif
-
- if (TRunoff%qsub(nr,nt) < 0._r8) then
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + TRunoff%qsub(nr,nt)
- TRunoff%qsub(nr,nt) = 0._r8
- endif
-
- if (TRunoff%qsur(nr,nt) < 0._r8) then
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + TRunoff%qsur(nr,nt)
- TRunoff%qsur(nr,nt) = 0._r8
- endif
-
-! if (TUnit%mask(nr) > 0 .and. TUnit%areaTotal(nr) > 0._r8) then
- if (TUnit%mask(nr) > 0) then
- ! mosart euler
- else
- rtmCTL%direct(nr,nt) = rtmCTL%direct(nr,nt) + &
- TRunoff%qsub(nr,nt) + &
- TRunoff%qsur(nr,nt) + &
- TRunoff%qgwl(nr,nt)
- TRunoff%qsub(nr,nt) = 0._r8
- TRunoff%qsur(nr,nt) = 0._r8
- TRunoff%qgwl(nr,nt) = 0._r8
- endif
-
- enddo
- enddo
-
!-----------------------------------
! MOSART Subcycling
!-----------------------------------
call t_startf('mosartr_subcycling')
+ if (first_call) then
+ do nt = 1,nt_rtm
+ write(iulog,'(2a,i6,l4)') trim(subname),' euler_calc for nt = ',nt,TUnit%euler_calc(nt)
+ enddo
+ endif
+
nsub = coupling_period/delt_mosart
if (nsub*delt_mosart < coupling_period) then
nsub = nsub + 1
@@ -1660,7 +1686,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
! --- convert TRunoff fields from m3/s to m/s before calling Euler
!-----------------------------------
- if (budget_check) then
+! if (budget_check) then
call t_startf('mosartr_budget')
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1670,7 +1696,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
enddo
enddo
call t_stopf('mosartr_budget')
- endif
+! endif
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1774,11 +1800,13 @@ subroutine Rtmrun(rstwr,nlend,rdate)
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%runofflnd(nr,nt) = rtmCTL%runoff(nr,nt)
rtmCTL%dvolrdtlnd(nr,nt)= rtmCTL%dvolrdt(nr,nt)
elseif (rtmCTL%mask(nr) >= 2) then
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)
endif
enddo
@@ -1794,7 +1822,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
! BUDGET terms 1-10 are for volumes (m3)
! BUDGET terms 11-30 are for flows (m3/s)
! BUDGET only ocean runoff and direct gets out of the system
- if (budget_check) then
+! if (budget_check) then
call t_startf('mosartr_budget')
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1822,6 +1850,22 @@ subroutine Rtmrun(rstwr,nlend,rdate)
budget_terms(22,nt) = budget_terms(22,nt) + rtmCTL%flood(nr)
enddo
+ ! accumulate the budget total over the run to make sure it's decreasing on avg
+ budget_accum_cnt = budget_accum_cnt + 1
+ do nt = 1,nt_rtm
+ budget_volume = (budget_terms( 2,nt) - budget_terms( 1,nt)) / delt_coupling
+ budget_input = (budget_terms(13,nt) + budget_terms(14,nt) + &
+ budget_terms(15,nt) + budget_terms(16,nt))
+ budget_output = (budget_terms(18,nt) + budget_terms(19,nt) + &
+ budget_terms(21,nt))
+ budget_total = budget_volume - budget_input + budget_output
+ budget_accum(nt) = budget_accum(nt) + budget_total
+ budget_terms(30,nt) = budget_accum(nt)/budget_accum_cnt
+ enddo
+ call t_stopf('mosartr_budget')
+
+ if (budget_check) then
+ call t_startf('mosartr_budget')
!--- check budget
! convert fluxes from m3/s to m3 by mult by coupling_period
@@ -1876,6 +1920,7 @@ subroutine Rtmrun(rstwr,nlend,rdate)
write(iulog,'(2a,i4,f22.6)') trim(subname),' net (dv-i+o) = ',nt,budget_total
!write(iulog,'(2a,i4,f22.6)') trim(subname),' net euler = ',nt,budget_euler
write(iulog,'(2a,i4,f22.6)') trim(subname),' eul erout lag = ',nt,budget_eroutlag
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' accum (dv-i+o)= ',nt,budget_global(30,nt)
!write(iulog,'(2a)') trim(subname),'----------------'
!write(iulog,'(2a,i4,f22.6)') trim(subname),' erout_prev no= ',nt,budget_global(23,nt)
!write(iulog,'(2a,i4,f22.6)') trim(subname),' erout no= ',nt,budget_global(24,nt)
@@ -1887,9 +1932,12 @@ subroutine Rtmrun(rstwr,nlend,rdate)
!write(iulog,'(2a,i4,f22.6)') trim(subname),' net main chan = ',nt,budget_global(6,nt)-budget_global(5,nt)+budget_global(24,nt)-budget_global(23,nt)+budget_global(27,nt)+budget_global(28,nt)+budget_global(29,nt)
!write(iulog,'(2a)') trim(subname),'----------------'
- if ((budget_total+budget_eroutlag) /= 0._r8) then
- if ((budget_total-budget_eroutlag)/(budget_total+budget_eroutlag) > 0.01_r8) then
- write(iulog,'(2a,i4)') trim(subname),' ***** BUDGET WARNING for nt = ',nt
+ if ((budget_total-budget_eroutlag) > 1.0e-6) then
+ write(iulog,'(2a,i4)') trim(subname),' ***** BUDGET WARNING error gt 1. m3 for nt = ',nt
+ endif
+ if ((budget_total+budget_eroutlag) >= 1.0e-6) then
+ if ((budget_total-budget_eroutlag)/(budget_total+budget_eroutlag) > 0.001_r8) then
+ write(iulog,'(2a,i4)') trim(subname),' ***** BUDGET WARNING out of balance for nt = ',nt
endif
endif
enddo
@@ -2054,6 +2102,8 @@ subroutine MOSART_init
integer :: igrow, igcol, iwgt
type(mct_avect) :: avtmp, avtmpG ! temporary avects
type(mct_sMat) :: sMat ! temporary sparse matrix, needed for sMatP
+ real(r8):: areatot_prev, areatot_tmp, areatot_new
+ integer :: tcnt
character(len=16384) :: rList ! list of fields for SM multiply
character(len=1000) :: fname
character(len=*),parameter :: subname = '(MOSART_init)'
@@ -2084,6 +2134,9 @@ subroutine MOSART_init
deallocate(compdof)
call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
+ allocate(TUnit%euler_calc(nt_rtm))
+ Tunit%euler_calc = .true.
+
allocate(TUnit%frac(begr:endr))
ier = pio_inq_varid(ncid, name='frac', vardesc=vardesc)
call pio_read_darray(ncid, vardesc, iodesc_dbl, TUnit%frac, ier)
@@ -2534,6 +2587,71 @@ subroutine MOSART_init
end if ! endr >= begr
+ !--- compute areatot from area using dnID ---
+ !--- this basically advects upstream areas downstream and
+ !--- adds them up as it goes until all upstream areas are accounted for
+
+ allocate(Tunit%areatotal2(rtmCTL%begr:rtmCTL%endr))
+ Tunit%areatotal2 = 0._r8
+
+ ! initialize avdst to local area and add that to areatotal2
+ cnt = 0
+ call mct_avect_zero(avdst_eroutUp)
+ do nr = rtmCTL%begr,rtmCTL%endr
+ cnt = cnt + 1
+ avdst_eroutUp%rAttr(1,cnt) = rtmCTL%area(nr)
+ Tunit%areatotal2(nr) = avdst_eroutUp%rAttr(1,cnt)
+ enddo
+
+ tcnt = 0
+ areatot_prev = -99._r8
+ areatot_new = -50._r8
+ do while (areatot_new /= areatot_prev .and. tcnt < rtmlon*rtmlat)
+
+ tcnt = tcnt + 1
+
+ ! 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
+ areatot_prev = areatot_new
+ areatot_tmp = 0._r8
+ do nr = rtmCTL%begr,rtmCTL%endr
+ cnt = cnt + 1
+ Tunit%areatotal2(nr) = Tunit%areatotal2(nr) + avdst_eroutUp%rAttr(1,cnt)
+ areatot_tmp = areatot_tmp + Tunit%areatotal2(nr)
+ enddo
+ call shr_mpi_sum(areatot_tmp, areatot_new, mpicom_rof, 'areatot_new', all=.true.)
+
+ if (masterproc) then
+ write(iulog,*) trim(subname),' areatot calc ',tcnt,areatot_new
+ endif
+
+ enddo
+
+ if (areatot_new /= areatot_prev) then
+ write(iulog,*) trim(subname),' MOSART ERROR: areatot incorrect ',areatot_new, areatot_prev
+ call shr_sys_abort(trim(subname)//' ERROR areatot incorrect')
+ endif
+
+! do nr = rtmCTL%begr,rtmCTL%endr
+! if (TUnit%areatotal(nr) > 0._r8 .and. Tunit%areatotal2(nr) /= TUnit%areatotal(nr)) then
+! write(iulog,'(2a,i12,2e16.4,f16.4)') trim(subname),' areatot diff ',nr,TUnit%areatotal(nr),Tunit%areatota!l2(nr),&
+! abs(TUnit%areatotal(nr)-Tunit%areatotal2(nr))/(TUnit%areatotal(nr))
+! endif
+! enddo
+
+
! control parameters
Tctl%RoutingMethod = 1
!Tctl%DATAH = rtm_nsteps*get_step_size()
@@ -2579,7 +2697,7 @@ subroutine SubTimestep
do iunit=rtmCTL%begr,rtmCTL%endr
if(TUnit%mask(iunit) > 0 .and. TUnit%rlen(iunit) > 0._r8) then
- TUnit%phi_r(iunit) = TUnit%areaTotal(iunit)*sqrt(TUnit%rslp(iunit))/(TUnit%rlen(iunit)*TUnit%rwidth(iunit))
+ TUnit%phi_r(iunit) = TUnit%areaTotal2(iunit)*sqrt(TUnit%rslp(iunit))/(TUnit%rlen(iunit)*TUnit%rwidth(iunit))
if(TUnit%phi_r(iunit) >= 10._r8) then
TUnit%numDT_r(iunit) = (TUnit%numDT_r(iunit)*log10(TUnit%phi_r(iunit))*Tctl%DLevelR) + 1
else
diff --git a/models/rof/mosart/src/riverroute/RunoffMod.F90 b/models/rof/mosart/src/riverroute/RunoffMod.F90
index 5ddebd00e99f..fcbcdc492e62 100644
--- a/models/rof/mosart/src/riverroute/RunoffMod.F90
+++ b/models/rof/mosart/src/riverroute/RunoffMod.F90
@@ -56,6 +56,7 @@ 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 :: 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)
@@ -84,6 +85,10 @@ module RunoffMod
real(r8), pointer :: runofflnd_nt2(:)
real(r8), pointer :: runoffocn_nt1(:)
real(r8), pointer :: runoffocn_nt2(:)
+ real(r8), pointer :: runofftot_nt1(:)
+ real(r8), pointer :: runofftot_nt2(:)
+ real(r8), pointer :: runoffdir_nt1(:)
+ real(r8), pointer :: runoffdir_nt2(:)
real(r8), pointer :: dvolrdtlnd_nt1(:)
real(r8), pointer :: dvolrdtlnd_nt2(:)
real(r8), pointer :: dvolrdtocn_nt1(:)
@@ -146,9 +151,11 @@ module RunoffMod
real(r8), pointer :: lon(:) ! longitude of the centroid of the cell
real(r8), pointer :: area(:) ! area of local cell, [m2]
real(r8), pointer :: areaTotal(:) ! total upstream drainage area, [m2]
+ real(r8), pointer :: areaTotal2(:)! computed total upstream drainage area, [m2]
real(r8), pointer :: rlenTotal(:) ! length of all reaches, [m]
real(r8), pointer :: Gxr(:) ! drainage density within the cell, [1/m]
real(r8), pointer :: frac(:) ! fraction of cell included in the study area, [-]
+ logical , pointer :: euler_calc(:) ! flag for calculating tracers in euler
! hillslope properties
real(r8), pointer :: nh(:) ! manning's roughness of the hillslope (channel network excluded)
@@ -289,6 +296,7 @@ subroutine RunoffInit(begr, endr, numr)
rtmCTL%dvolrdtlnd(begr:endr,nt_rtm), &
rtmCTL%runoffocn(begr:endr,nt_rtm), &
rtmCTL%dvolrdtocn(begr:endr,nt_rtm), &
+ rtmCTL%runofftot(begr:endr,nt_rtm), &
rtmCTL%area(begr:endr), &
rtmCTL%volr(begr:endr,nt_rtm), &
rtmCTL%lonc(begr:endr), &
@@ -299,6 +307,10 @@ subroutine RunoffInit(begr, endr, numr)
rtmCTL%runofflnd_nt2(begr:endr), &
rtmCTL%runoffocn_nt1(begr:endr), &
rtmCTL%runoffocn_nt2(begr:endr), &
+ rtmCTL%runofftot_nt1(begr:endr), &
+ rtmCTL%runofftot_nt2(begr:endr), &
+ rtmCTL%runoffdir_nt1(begr:endr), &
+ rtmCTL%runoffdir_nt2(begr:endr), &
rtmCTL%volr_nt1(begr:endr), &
rtmCTL%volr_nt2(begr:endr), &
rtmCTL%dvolrdtlnd_nt1(begr:endr), &
@@ -335,6 +347,7 @@ subroutine RunoffInit(begr, endr, numr)
rtmCTL%runoff(:,:) = 0._r8
rtmCTL%runofflnd(:,:) = spval
rtmCTL%runoffocn(:,:) = spval
+ rtmCTL%runofftot(:,:) = spval
rtmCTL%dvolrdt(:,:) = 0._r8
rtmCTL%dvolrdtlnd(:,:) = spval
rtmCTL%dvolrdtocn(:,:) = spval