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 d68cc18278b7..558ea97aa92f 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_20130604a.nc
+rof/mosart/MOSART_Global_half_20151015a.nc
rof/mosart/MOSART_NLDAS_8th.nc
diff --git a/models/rof/mosart/doc/ChangeLog b/models/rof/mosart/doc/ChangeLog
index eafd10d56f08..b68d9bfaf2bc 100644
--- a/models/rof/mosart/doc/ChangeLog
+++ b/models/rof/mosart/doc/ChangeLog
@@ -1,3 +1,26 @@
+===============================================================
+Tag name: mosart1_0_09
+Originator(s): tcraig
+Date: Nov 29 2015
+One-line Summary: Code cleanup, add budget diagnostics, history files
+
+M src/riverroute/RtmMod.F90
+M src/riverroute/RtmHistFlds.F90
+M src/riverroute/MOSART_physics_mod.F90
+M src/riverroute/RtmHistFile.F90
+M src/riverroute/RunoffMod.F90
+M src/riverroute/RtmRestFile.F90
+M src/cpl/rof_cpl_indices.F90
+M src/cpl/rof_comp_mct.F90
+M src/cpl/rof_comp_esmf.F90
+===============================================================
+Tag name: mosart1_0_08
+Originator(s): tcraig
+Date: Nov 24 2015
+One-line Summary: Fix exact restart in atan slope calc
+
+M src/riverroute/MOSART_physics_mod.F90
+
===============================================================
Tag name: mosart1_0_07
Originator(s): tcraig
diff --git a/models/rof/mosart/src/cpl/rof_comp_esmf.F90 b/models/rof/mosart/src/cpl/rof_comp_esmf.F90
index 3756e94987e4..ee20f74a5fff 100644
--- a/models/rof/mosart/src/cpl/rof_comp_esmf.F90
+++ b/models/rof/mosart/src/cpl/rof_comp_esmf.F90
@@ -2,7 +2,7 @@ module rof_comp_esmf
#ifdef ESMF_INTERFACE
!===============================================================================
-! Interface of the active river runoff componetn (RTM) model component of CESM
+! Interface of the active river runoff componetn (mosart) model component of CESM
! with the main CESM driver. This is a thin interface taking CESM driver information
! in MCT (Model Coupling Toolkit) format and converting it to use by RTM and outputing
! if in ESMF (Earth System Modelling Framework) format.
@@ -65,11 +65,6 @@ module rof_comp_esmf
private :: rof_export_esmf ! Export the river runoff model data to the CESM coupler
!
! ! PRIVATE DATA MEMBERS:
- real(r8), pointer :: totrunin(:,:) ! total runoff on rtm grid (mm/s)
- real(r8), pointer :: surrunin(:,:) ! surface runoff on rtm grid (mm/s)
- real(r8), pointer :: subrunin(:,:) ! subsurface runoff on rtm grid (mm/s)
- real(r8), pointer :: gwlrunin(:,:) ! glacier, wetlands and lakes water balance residual on rtm grid (mm/s)
- real(r8), pointer :: dtorunin(:,:) ! direct to ocean runoff on rtm grid (mm/s)
!
!===============================================================================
@@ -111,13 +106,13 @@ subroutine rof_init_esmf(comp, import_state, export_state, EClock, rc)
!---------------------------------------------------------------------------
! !DESCRIPTION:
- ! Initialize runoff model (rtm)
+ ! Initialize runoff model (mosart)
!
! !ARGUMENTS:
implicit none
- type(ESMF_GridComp) :: comp ! RTM gridded component
- type(ESMF_State) :: import_state ! RTM import state
- type(ESMF_State) :: export_state ! RTM export state
+ type(ESMF_GridComp) :: comp ! MOSART gridded component
+ type(ESMF_State) :: import_state ! MOSART import state
+ type(ESMF_State) :: export_state ! MOSART export state
type(ESMF_Clock) :: EClock ! ESMF synchronization clock
integer, intent(out) :: rc ! Return code
!
@@ -171,7 +166,7 @@ subroutine rof_init_esmf(comp, import_state, export_state, EClock, rc)
call MPI_Comm_dup(mpicom_vm, mpicom_rof, rc)
if(rc /= 0) call ESMF_Finalize(rc=rc, endflag=ESMF_END_ABORT)
- ! Initialize rtm MPI communicator
+ ! Initialize mosart MPI communicator
call ESMF_AttributeGet(export_state, name="ID", value=ROFID, rc=rc)
if (rc /= ESMF_SUCCESS) call ESMF_Finalize(rc=rc, endflag=ESMF_END_ABORT)
@@ -207,15 +202,15 @@ subroutine rof_init_esmf(comp, import_state, export_state, EClock, rc)
call shr_file_setLogUnit (iulog)
if (masterproc) then
- write(iulog,*) ' rtm npes = ',npes
- write(iulog,*) ' rtm iam = ',iam
- write(iulog,*) ' rtm ROFID = ',ROFID
- write(iulog,*) ' rtm name = ',trim(inst_name)
- write(iulog,*) ' rtm inst = ',inst_index
- write(iulog,*) ' rtm suffix = ',trim(inst_suffix)
+ write(iulog,*) ' mosart npes = ',npes
+ write(iulog,*) ' mosart iam = ',iam
+ write(iulog,*) ' mosart ROFID = ',ROFID
+ write(iulog,*) ' mosart name = ',trim(inst_name)
+ write(iulog,*) ' mosart inst = ',inst_index
+ write(iulog,*) ' mosart suffix = ',trim(inst_suffix)
endif
- ! Initialize RTM
+ ! Initialize MOSART
call seq_timemgr_EClockGetData(EClock, &
start_ymd=start_ymd, &
@@ -266,18 +261,13 @@ subroutine rof_init_esmf(comp, import_state, export_state, EClock, rc)
nsrest_in=nsrest, version_in=version, &
hostname_in=hostname, username_in=username)
- ! Initialize rtm and determine if rtm will be active
+ ! Initialize mosart and determine if mosart will be active
call Rtmini(rtm_active=rof_prognostic,flood_active=flood_present)
if ( rof_prognostic ) then
begr = rtmCTL%begr
endr = rtmCTL%endr
- allocate(totrunin(begr:endr,nt_rtm))
- allocate(surrunin(begr:endr,nt_rtm))
- allocate(subrunin(begr:endr,nt_rtm))
- allocate(gwlrunin(begr:endr,nt_rtm))
- allocate(dtorunin(begr:endr,nt_rtm))
! Initialize rof distgrid and domain
@@ -353,17 +343,17 @@ subroutine rof_init_esmf(comp, import_state, export_state, EClock, rc)
call ESMF_AttributeAdd(comp, &
convention=convCIM, purpose=purpComp, rc=rc)
- call ESMF_AttributeSet(comp, "ShortName", "RTM", &
+ call ESMF_AttributeSet(comp, "ShortName", "MOSART", &
convention=convCIM, purpose=purpComp, rc=rc)
call ESMF_AttributeSet(comp, "LongName", &
- "River Runoff from RTM", &
+ "River Runoff from MOSART", &
convention=convCIM, purpose=purpComp, rc=rc)
call ESMF_AttributeSet(comp, "Description", &
- "The RTM component is " // &
+ "The MOSART component is " // &
"the river runoff model used in the CESM1.1. " // &
- "More information on the RTM project " // &
- "and access to previous RTM model versions and " // &
- "documentation can be found via the RTM Web Page.", &
+ "More information on the MOSART project " // &
+ "and access to previous MOSART model versions and " // &
+ "documentation can be found via the MOSART Web Page.", &
convention=convCIM, purpose=purpComp, rc=rc)
call ESMF_AttributeSet(comp, "ReleaseDate", "2012", &
convention=convCIM, purpose=purpComp, rc=rc)
@@ -387,13 +377,13 @@ subroutine rof_run_esmf(comp, import_state, export_state, EClock, rc)
!---------------------------------------------------------------------------
! !DESCRIPTION:
- ! Run rtm model
+ ! Run mosart model
!
! !ARGUMENTS:
implicit none
- type(ESMF_GridComp) :: comp ! RTM gridded component
- type(ESMF_State) :: import_state ! RTM import state
- type(ESMF_State) :: export_state ! RTM export state
+ type(ESMF_GridComp) :: comp ! MOSART gridded component
+ type(ESMF_State) :: import_state ! MOSART import state
+ type(ESMF_State) :: export_state ! MOSART export state
type(ESMF_Clock) :: EClock ! ESMF synchronization clock
integer, intent(out) :: rc ! Return code
!
@@ -404,11 +394,11 @@ subroutine rof_run_esmf(comp, import_state, export_state, EClock, rc)
integer :: mon_sync ! Sync current month
integer :: day_sync ! Sync current day
integer :: tod_sync ! Sync current time of day (sec)
- integer :: ymd ! RTM current date (YYYYMMDD)
- integer :: yr ! RTM current year
- integer :: mon ! RTM current month
- integer :: day ! RTM current day
- integer :: tod ! RTM current time of day (sec)
+ integer :: ymd ! MOSART current date (YYYYMMDD)
+ integer :: yr ! MOSART current year
+ integer :: mon ! MOSART current month
+ integer :: day ! MOSART current day
+ integer :: tod ! MOSART current time of day (sec)
logical :: rstwr_sync ! .true. ==> write restart file before returning
logical :: rstwr ! .true. ==> write restart file before returning
logical :: nlend_sync ! Flag signaling last time-step
@@ -441,7 +431,7 @@ subroutine rof_run_esmf(comp, import_state, export_state, EClock, rc)
curr_ymd=ymd, curr_tod=tod_sync, &
curr_yr=yr_sync, curr_mon=mon_sync, curr_day=day_sync)
- ! Map ESMF to rtm input (rof) data type
+ ! Map ESMF to mosart input (rof) data type
call t_startf ('lc_rof_import')
call ESMF_StateGet(import_state, itemName="x2d", array=x2r, rc=rc)
@@ -449,14 +439,14 @@ subroutine rof_run_esmf(comp, import_state, export_state, EClock, rc)
call rof_import_esmf( x2r, rc=rc )
call t_stopf ('lc_rof_import')
- ! Run rtm (input is *runin, output is rtmCTL%runoff)
- ! First advance rtm time step
+ ! Run mosart (input is *runin, output is rtmCTL%runoff)
+ ! First advance mosart time step
write(rdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr_sync,mon_sync,day_sync,tod_sync
nlend = seq_timemgr_StopAlarmIsOn( EClock )
rstwr = seq_timemgr_RestartAlarmIsOn( EClock )
call advance_timestep()
- call Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate)
+ call Rtmrun(rstwr,nlend,rdate)
! Map roff data to MCT datatype (input is rtmCTL%runoff, output is r2x_r)
@@ -472,9 +462,9 @@ subroutine rof_run_esmf(comp, import_state, export_state, EClock, rc)
tod = tod
if ( .not. seq_timemgr_EClockDateInSync( EClock, ymd, tod ) )then
call seq_timemgr_EclockGetData( EClock, curr_ymd=ymd_sync, curr_tod=tod_sync )
- write(iulog,*)' rtm ymd=',ymd ,' rtm tod= ',tod
+ write(iulog,*)' mosart ymd=',ymd ,' mosart tod= ',tod
write(iulog,*)'sync ymd=',ymd_sync,' sync tod= ',tod_sync
- call shr_sys_abort( sub//":: RTM clock not in sync with Master Sync clock" )
+ call shr_sys_abort( sub//":: MOSART clock not in sync with Master Sync clock" )
end if
! Reset shr logging to my original values
@@ -498,13 +488,13 @@ subroutine rof_final_esmf(comp, import_state, export_state, EClock, rc)
!------------------------------------------------------------------------------
! !DESCRIPTION:
- ! Finalize rtm
+ ! Finalize mosart
!
! !ARGUMENTS:
implicit none
- type(ESMF_GridComp) :: comp ! RTM gridded component
- type(ESMF_State) :: import_state ! RTM import state
- type(ESMF_State) :: export_state ! RTM export state
+ type(ESMF_GridComp) :: comp ! MOSART gridded component
+ type(ESMF_State) :: import_state ! MOSART import state
+ type(ESMF_State) :: export_state ! MOSART export state
type(ESMF_Clock) :: EClock ! ESMF synchronization clock
integer, intent(out) :: rc ! Return code
!---------------------------------------------------------------------------
@@ -655,6 +645,7 @@ subroutine rof_import_esmf( x2r_array, rc)
!---------------------------------------------------------------------------
! DESCRIPTION:
! Obtain the runoff input from the coupler
+ ! convert from kg/m2s to m3/s
!
! ARGUMENTS:
implicit none
@@ -691,34 +682,21 @@ subroutine rof_import_esmf( x2r_array, rc)
endr = rtmCTL%endr
do n = begr,endr
n2 = n - begr + 1
- totrunin(n,nliq) = fptr(index_x2r_Flrl_rofsur,n2) + &
- fptr(index_x2r_Flrl_rofgwl,n2) + &
- fptr(index_x2r_Flrl_rofsub,n2)
- totrunin(n,nfrz) = fptr(index_x2r_Flrl_rofi,n2)
- surrunin(n,nliq) = fptr(index_x2r_Flrl_rofsur,n2)
- surrunin(n,nfrz) = fptr(index_x2r_Flrl_rofi,n2)
-
- subrunin(n,nliq) = fptr(index_x2r_Flrl_rofsub,n2)
- subrunin(n,nfrz) = 0.0_r8
-
- gwlrunin(n,nliq) = fptr(index_x2r_Flrl_rofgwl,n2)
- gwlrunin(n,nfrz) = 0.0_r8
-
- dtorunin(n,nliq) = 0.0_r8
- dtorunin(n,nfrz) = 0.0_r8
-
- rtmCTL%qsur(n) = fptr(index_x2r_Flrl_rofsur,n2)
- rtmCTL%qsub(n) = fptr(index_x2r_Flrl_rofsub,n2)
- rtmCTL%qgwl(n) = fptr(index_x2r_Flrl_rofgwl,n2)
- rtmCTL%qdto(n) = 0.0_r8
+ rtmCTL%qsur(n,nliq) = fptr(index_x2r_Flrl_rofsur,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qsub(n,nliq) = fptr(index_x2r_Flrl_rofsub,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qgwl(n,nliq) = fptr(index_x2r_Flrl_rofgwl,n2) * (rtmCTL%area(n)*0.001_r8)
if (index_x2r_Flrl_rofdto > 0) then
- totrunin(n,nliq) = totrunin(n,nliq) + &
- fptr(index_x2r_Flrl_rofdto,n2)
- dtorunin(n,nliq) = fptr(index_x2r_Flrl_rofdto,n2)
- rtmCTL%qdto(n) = ftpr(index_x2r_Flrl_rofdto,n2)
+ rtmCTL%qdto(n,nliq) = fptr(index_x2r_Flrl_rofdto,n2) * (rtmCTL%area(n)*0.001_r8)
+ else
+ rtmCTL%qdto(n,nliq) = 0.0_r8
endif
+ rtmCTL%qsur(n,nfrz) = fptr(index_x2r_Flrl_rofi,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qsub(n,nfrz) = 0.0_r8
+ rtmCTL%qgwl(n,nfrz) = 0.0_r8
+ rtmCTL%qdto(n,nfrz) = 0.0_r8
+
enddo
end subroutine rof_import_esmf
@@ -729,7 +707,8 @@ subroutine rof_export_esmf( r2x_array, rc)
!---------------------------------------------------------------------------
! !DESCRIPTION:
- ! Send the rtm export state to the CESM coupler
+ ! Send the mosart export state to the CESM coupler
+ ! convert from kg/m2s to m3/s
!
! !ARGUMENTS:
implicit none
@@ -739,6 +718,7 @@ subroutine rof_export_esmf( r2x_array, rc)
! Local variables
integer :: ni, n, nt, nliq, nfrz
real(R8), pointer :: fptr(:, :)
+ logical,save :: first_time = .true.
character(len=*), parameter :: sub = 'rof_export_esmf'
!---------------------------------------------------------------------------
@@ -763,16 +743,29 @@ subroutine rof_export_esmf( r2x_array, rc)
if (rc /= ESMF_SUCCESS) call ESMF_Finalize(rc=rc, endflag=ESMF_END_ABORT)
fptr(:,:) = 0._r8
+ if (first_time) then
+ if (masterproc) then
+ if ( ice_runoff )then
+ write(iulog,*)'Snow capping will flow out in frozen river runoff'
+ else
+ write(iulog,*)'Snow capping will flow out in liquid river runoff'
+ endif
+ endif
+ first_time = .false.
+ end if
+
ni = 0
if ( ice_runoff )then
+ ! separate liquid and ice runoff
do n = rtmCTL%begr,rtmCTL%endr
ni = ni + 1
+ fptr(index_r2x_Forr_rofl,ni) = rtmCTL%direct(n,nliq) / (rtmCTL%area(n)*0.001_r8)
+ fptr(index_r2x_Forr_rofi,ni) = rtmCTL%direct(n,nfrz) / (rtmCTL%area(n)*0.001_r8)
if (rtmCTL%mask(n) >= 2) then
- ! liquid and ice runoff are treated separately
- fptr(index_r2x_Forr_rofl,ni) = &
- rtmCTL%runoffall(n,nliq)/(rtmCTL%area(n)*0.001_r8)
- fptr(index_r2x_Forr_rofi,ni) = &
- rtmCTL%runoffall(n,nfrz)/(rtmCTL%area(n)*0.001_r8)
+ fptr(index_r2x_Forr_rofl,ni) = fptr(index_r2x_Forr_rofl,ni) + &
+ rtmCTL%runoff(n,nliq) / (rtmCTL%area(n)*0.001_r8)
+ fptr(index_r2x_Forr_rofi,ni) = fptr(index_r2x_Forr_rofi,ni) + &
+ rtmCTL%runoff(n,nfrz) / (rtmCTL%area(n)*0.001_r8)
if (ni > rtmCTL%lnumr) then
write(iulog,*) sub, ' : ERROR runoff count',n,ni
call shr_sys_abort( sub//' : ERROR runoff > expected' )
@@ -780,13 +773,14 @@ subroutine rof_export_esmf( r2x_array, rc)
endif
end do
else
+ ! liquid and ice runoff added to liquid runoff, ice runoff is zero
do n = rtmCTL%begr,rtmCTL%endr
ni = ni + 1
+ fptr(index_r2x_Forr_rofl,ni) = &
+ (rtmCTL%direct(n,nfrz)+rtmCTL%direct(n,nliq)) / (rtmCTL%area(n)*0.001_r8)
if (rtmCTL%mask(n) >= 2) then
- ! liquid and ice runoff are bundled together to liquid runoff, and then ice runoff set to zero
- fptr(index_r2x_Forr_rofl,ni) = &
- (rtmCTL%runoffall(n,nfrz)+rtmCTL%runoffall(n,nliq))/(rtmCTL%area(n)*0.001_r8)
- fptr(index_r2x_Forr_rofi,ni) = 0.0_r8
+ fptr(index_r2x_Forr_rofl,ni) = fptr(index_r2x_Forr_rofl,ni) + &
+ (rtmCTL%runoff(n,nfrz)+rtmCTL%runoff(n,nliq))/(rtmCTL%area(n)*0.001_r8)
if (ni > rtmCTL%lnumr) then
write(iulog,*) sub, ' : ERROR runoff count',n,ni
call shr_sys_abort( sub//' : ERROR runoff > expected' )
@@ -795,15 +789,6 @@ subroutine rof_export_esmf( r2x_array, rc)
end do
end if
-! Add direct to ocean runoff prior to passing runoff to coupler
-! tcraig, This is now done in rtmrun to include in the budget
-! ni = 0
-! do n = rtmCTL%begr,rtmCTL%endr
-! ni = ni + 1
-! fptr(index_r2x_Forr_rofl,ni) &
-! = fptr(index_r2x_Forr_rofl,ni) + rtmCTL%qdto(n)
-! end do
-
! Flooding back to land, sign convention is positive in land->rof direction
! so if water is sent from rof to land, the flux must be negative.
ni = 0
diff --git a/models/rof/mosart/src/cpl/rof_comp_mct.F90 b/models/rof/mosart/src/cpl/rof_comp_mct.F90
index afeb7c6abc53..63e2cbf210d5 100644
--- a/models/rof/mosart/src/cpl/rof_comp_mct.F90
+++ b/models/rof/mosart/src/cpl/rof_comp_mct.F90
@@ -4,7 +4,7 @@ module rof_comp_mct
! DESCRIPTION:
! Interface of the active runoff component of CESM
! with the main CESM driver. This is a thin interface taking CESM driver information
-! in MCT (Model Coupling Toolkit) format and converting it to use by RTM
+! in MCT (Model Coupling Toolkit) format and converting it to use by MOSART
use seq_flds_mod
use shr_kind_mod , only : r8 => shr_kind_r8
@@ -56,11 +56,6 @@ module rof_comp_mct
private :: rof_export_mct ! Export the river runoff model data to the CESM coupler
!
! PRIVATE DATA MEMBERS:
- real(r8), pointer :: totrunin(:,:) ! total runoff on rtm grid (mm/s)
- real(r8), pointer :: subrunin(:,:) ! subsurface runoff on rtm grid (mm/s)
- real(r8), pointer :: surrunin(:,:) ! surface runoff on rtm grid (mm/s)
- real(r8), pointer :: gwlrunin(:,:) ! water residual from glacier, wetlands and lakes water balance on rtm grid (mm/s)
- real(r8), pointer :: dtorunin(:,:) ! direct to ocean runoff on rtm grid (mm/s)
! REVISION HISTORY:
! Author: Mariana Vertenstein
@@ -124,7 +119,7 @@ subroutine rof_init_mct( EClock, cdata_r, x2r_r, r2x_r, NLFilename)
! Determine attriute vector indices
call rof_cpl_indices_set()
- ! Initialize rtm MPI communicator
+ ! Initialize mosart MPI communicator
call RtmSpmdInit(mpicom_loc)
#if (defined _MEMTRACE)
@@ -155,12 +150,12 @@ subroutine rof_init_mct( EClock, cdata_r, x2r_r, r2x_r, NLFilename)
call shr_file_setLogUnit (iulog)
if (masterproc) then
- write(iulog,*) ' rtm npes = ',npes
- write(iulog,*) ' rtm iam = ',iam
+ write(iulog,*) ' mosart npes = ',npes
+ write(iulog,*) ' mosart iam = ',iam
write(iulog,*) ' inst_name = ',trim(inst_name)
endif
- ! Initialize rtm
+ ! Initialize mosart
call seq_timemgr_EClockGetData(EClock, &
start_ymd=start_ymd, &
start_tod=start_tod, ref_ymd=ref_ymd, &
@@ -201,11 +196,6 @@ subroutine rof_init_mct( EClock, cdata_r, x2r_r, r2x_r, NLFilename)
! Initialize memory for input state
begr = rtmCTL%begr
endr = rtmCTL%endr
- allocate (totrunin(begr:endr,nt_rtm))
- allocate (subrunin(begr:endr,nt_rtm))
- allocate (surrunin(begr:endr,nt_rtm))
- allocate (gwlrunin(begr:endr,nt_rtm))
- allocate (dtorunin(begr:endr,nt_rtm))
! Initialize rof gsMap for ocean rof and land rof
call rof_SetgsMap_mct( mpicom_rof, ROFID, gsMap_rof)
@@ -214,11 +204,11 @@ subroutine rof_init_mct( EClock, cdata_r, x2r_r, r2x_r, NLFilename)
lsize = mct_gsMap_lsize(gsMap_rof, mpicom_rof)
call rof_domain_mct( lsize, gsMap_rof, dom_r )
- ! Initialize lnd -> rtm attribute vector
+ ! Initialize lnd -> mosart attribute vector
call mct_aVect_init(x2r_r, rList=seq_flds_x2r_fields, lsize=lsize)
call mct_aVect_zero(x2r_r)
- ! Initialize rtm -> ocn attribute vector
+ ! Initialize mosart -> ocn attribute vector
call mct_aVect_init(r2x_r, rList=seq_flds_r2x_fields, lsize=lsize)
call mct_aVect_zero(r2x_r)
@@ -302,13 +292,13 @@ subroutine rof_run_mct( EClock, cdata_r, x2r_r, r2x_r)
call rof_import_mct( x2r_r)
call t_stopf ('lc_rof_import')
- ! Run rtm (input is *runin, output is rtmCTL%runoff)
- ! First advance rtm time step
+ ! Run mosart (input is *runin, output is rtmCTL%runoff)
+ ! First advance mosart time step
write(rdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr_sync,mon_sync,day_sync,tod_sync
nlend = seq_timemgr_StopAlarmIsOn( EClock )
rstwr = seq_timemgr_RestartAlarmIsOn( EClock )
call advance_timestep()
- call Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate)
+ call Rtmrun(rstwr,nlend,rdate)
! Map roff data to MCT datatype (input is rtmCTL%runoff, output is r2x_r)
call t_startf ('lc_rof_export')
@@ -321,9 +311,9 @@ subroutine rof_run_mct( EClock, cdata_r, x2r_r, r2x_r)
tod = tod
if ( .not. seq_timemgr_EClockDateInSync( EClock, ymd, tod ) )then
call seq_timemgr_EclockGetData( EClock, curr_ymd=ymd_sync, curr_tod=tod_sync )
- write(iulog,*)' rtm ymd=',ymd ,' rtm tod= ',tod
+ write(iulog,*)' mosart ymd=',ymd ,' mosart tod= ',tod
write(iulog,*)'sync ymd=',ymd_sync,' sync tod= ',tod_sync
- call shr_sys_abort( sub//":: RTM clock is not in sync with Master Sync clock" )
+ call shr_sys_abort( sub//":: MOSART clock is not in sync with Master Sync clock" )
end if
! Reset shr logging to my original values
@@ -514,6 +504,7 @@ subroutine rof_import_mct( x2r_r)
!---------------------------------------------------------------------------
! DESCRIPTION:
! Obtain the runoff input from the coupler
+ ! convert from kg/m2s to m3/s
!
! ARGUMENTS:
implicit none
@@ -545,33 +536,21 @@ subroutine rof_import_mct( x2r_r)
endr = rtmCTL%endr
do n = begr,endr
n2 = n - begr + 1
- totrunin(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofsur,n2) + &
- x2r_r%rAttr(index_x2r_Flrl_rofgwl,n2) + &
- x2r_r%rAttr(index_x2r_Flrl_rofsub,n2)
- totrunin(n,nfrz) = x2r_r%rAttr(index_x2r_Flrl_rofi,n2)
- surrunin(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofsur,n2)
- surrunin(n,nfrz) = x2r_r%rAttr(index_x2r_Flrl_rofi,n2)
-
- subrunin(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofsub,n2)
- subrunin(n,nfrz) = 0.0_r8
-
- gwlrunin(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofgwl,n2)
- gwlrunin(n,nfrz) = 0.0_r8
-
- dtorunin(n,nliq) = 0.0_r8
- dtorunin(n,nfrz) = 0.0_r8
-
- rtmCTL%qsur(n) = x2r_r%rAttr(index_x2r_Flrl_rofsur,n2)
- rtmCTL%qsub(n) = x2r_r%rAttr(index_x2r_Flrl_rofsub,n2)
- rtmCTL%qgwl(n) = x2r_r%rAttr(index_x2r_Flrl_rofgwl,n2)
- rtmCTL%qdto(n) = 0.0_r8
+ rtmCTL%qsur(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofsur,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qsub(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofsub,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qgwl(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofgwl,n2) * (rtmCTL%area(n)*0.001_r8)
if (index_x2r_Flrl_rofdto > 0) then
- totrunin(n,nliq) = totrunin(n,nliq) + &
- x2r_r%rAttr(index_x2r_Flrl_rofdto,n2)
- dtorunin(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofdto,n2)
- rtmCTL%qdto(n) = x2r_r%rAttr(index_x2r_Flrl_rofdto,n2)
+ rtmCTL%qdto(n,nliq) = x2r_r%rAttr(index_x2r_Flrl_rofdto,n2) * (rtmCTL%area(n)*0.001_r8)
+ else
+ rtmCTL%qdto(n,nliq) = 0.0_r8
endif
+
+ rtmCTL%qsur(n,nfrz) = x2r_r%rAttr(index_x2r_Flrl_rofi,n2) * (rtmCTL%area(n)*0.001_r8)
+ rtmCTL%qsub(n,nfrz) = 0.0_r8
+ rtmCTL%qgwl(n,nfrz) = 0.0_r8
+ rtmCTL%qdto(n,nfrz) = 0.0_r8
+
enddo
end subroutine rof_import_mct
@@ -583,6 +562,7 @@ subroutine rof_export_mct( r2x_r )
!---------------------------------------------------------------------------
! DESCRIPTION:
! Send the runoff model export state to the coupler
+ ! convert from m3/s to kg/m2s
!
! ARGUMENTS:
implicit none
@@ -590,7 +570,7 @@ subroutine rof_export_mct( r2x_r )
!
! LOCAL VARIABLES
integer :: ni, n, nt, nliq, nfrz
- logical :: first_time = .true.
+ logical,save :: first_time = .true.
character(len=32), parameter :: sub = 'rof_export_mct'
!---------------------------------------------------------------------------
@@ -624,14 +604,17 @@ subroutine rof_export_mct( r2x_r )
ni = 0
if ( ice_runoff )then
+ ! separate liquid and ice runoff
do n = rtmCTL%begr,rtmCTL%endr
ni = ni + 1
+ r2x_r%rAttr(index_r2x_Forr_rofl,ni) = rtmCTL%direct(n,nliq) / (rtmCTL%area(n)*0.001_r8)
+ r2x_r%rAttr(index_r2x_Forr_rofi,ni) = rtmCTL%direct(n,nfrz) / (rtmCTL%area(n)*0.001_r8)
if (rtmCTL%mask(n) >= 2) then
! liquid and ice runoff are treated separately - this is what goes to the ocean
- r2x_r%rAttr(index_r2x_Forr_rofl,ni) = &
- rtmCTL%runoffall(n,nliq)/(rtmCTL%area(n)*0.001_r8)
- r2x_r%rAttr(index_r2x_Forr_rofi,ni) = &
- rtmCTL%runoffall(n,nfrz)/(rtmCTL%area(n)*0.001_r8)
+ r2x_r%rAttr(index_r2x_Forr_rofl,ni) = r2x_r%rAttr(index_r2x_Forr_rofl,ni) + &
+ rtmCTL%runoff(n,nliq) / (rtmCTL%area(n)*0.001_r8)
+ r2x_r%rAttr(index_r2x_Forr_rofi,ni) = r2x_r%rAttr(index_r2x_Forr_rofi,ni) + &
+ rtmCTL%runoff(n,nfrz) / (rtmCTL%area(n)*0.001_r8)
if (ni > rtmCTL%lnumr) then
write(iulog,*) sub, ' : ERROR runoff count',n,ni
call shr_sys_abort( sub//' : ERROR runoff > expected' )
@@ -639,13 +622,14 @@ subroutine rof_export_mct( r2x_r )
endif
end do
else
+ ! liquid and ice runoff added to liquid runoff, ice runoff is zero
do n = rtmCTL%begr,rtmCTL%endr
ni = ni + 1
+ r2x_r%rAttr(index_r2x_Forr_rofl,ni) = &
+ (rtmCTL%direct(n,nfrz)+rtmCTL%direct(n,nliq)) / (rtmCTL%area(n)*0.001_r8)
if (rtmCTL%mask(n) >= 2) then
- ! liquid and ice runoff are bundled together to liquid runoff, and then ice runoff set to zero
- r2x_r%rAttr(index_r2x_Forr_rofl,ni) = &
- (rtmCTL%runoffall(n,nfrz)+rtmCTL%runoffall(n,nliq))/(rtmCTL%area(n)*0.001_r8)
- r2x_r%rAttr(index_r2x_Forr_rofi,ni) = 0._r8
+ r2x_r%rAttr(index_r2x_Forr_rofl,ni) = r2x_r%rAttr(index_r2x_Forr_rofl,ni) + &
+ (rtmCTL%runoff(n,nfrz)+rtmCTL%runoff(n,nliq)) / (rtmCTL%area(n)*0.001_r8)
if (ni > rtmCTL%lnumr) then
write(iulog,*) sub, ' : ERROR runoff count',n,ni
call shr_sys_abort( sub//' : ERROR runoff > expected' )
@@ -654,21 +638,12 @@ subroutine rof_export_mct( r2x_r )
end do
end if
-! Add direct to ocean runoff prior to passing runoff to coupler
-! tcraig, This is now done in rtmrun to include in the budget
-! ni = 0
-! do n = rtmCTL%begr,rtmCTL%endr
-! ni = ni + 1
-! r2x_r%rAttr(index_r2x_Forr_rofl,ni) &
-! = r2x_r%rAttr(index_r2x_Forr_rofl,ni) + rtmCTL%qdto(n)
-! end do
-
! Flooding back to land, sign convention is positive in land->rof direction
! so if water is sent from rof to land, the flux must be negative.
ni = 0
do n = rtmCTL%begr, rtmCTL%endr
ni = ni + 1
- r2x_r%rattr(index_r2x_Flrr_flood,ni) = -rtmCTL%flood(n)/(rtmCTL%area(n)*0.001_r8)
+ r2x_r%rattr(index_r2x_Flrr_flood,ni) = -rtmCTL%flood(n) / (rtmCTL%area(n)*0.001_r8)
r2x_r%rattr(index_r2x_Flrr_volr,ni) = (Trunoff%wr(n,nliq) + Trunoff%wt(n,nliq)) / rtmCTL%area(n)
r2x_r%rattr(index_r2x_Flrr_volrmch,ni) = Trunoff%wr(n,nliq) / rtmCTL%area(n)
end do
diff --git a/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90 b/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
index c89c58be2bef..a23d690058ad 100644
--- a/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
+++ b/models/rof/mosart/src/riverroute/MOSART_physics_mod.F90
@@ -49,12 +49,6 @@ subroutine Euler
integer :: iunit, m, k, unitUp, cnt, ier !local index
real(r8) :: temp_erout, localDeltaT
real(r8) :: negchan
- logical, save :: first_call = .true.
-
- if (first_call) then
- sinatanSLOPE1defr = 1.0_r8/(sin(atan(SLOPE1def)))
- endif
- first_call = .false.
!------------------
! hillslope
@@ -74,9 +68,19 @@ subroutine Euler
call t_stopf('mosart_hillslope')
TRunoff%flow = 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
+ do iunit=rtmCTL%begr,rtmCTL%endr
+ TRunoff%erout_prev(iunit,nt) = TRunoff%erout_prev(iunit,nt) + TRunoff%erout(iunit,nt)
+ enddo
+ enddo
+
!------------------
! subnetwork
!------------------
@@ -110,7 +114,7 @@ subroutine Euler
endif
call t_startf('mosart_SMeroutUp')
- Trunoff%eroutUp = 0._r8
+ TRunoff%eroutUp = 0._r8
#ifdef NO_MCT
do iunit=rtmCTL%begr,rtmCTL%endr
do k=1,TUnit%nUp(iunit)
@@ -122,28 +126,32 @@ subroutine Euler
end do
#else
!--- copy erout into avsrc_eroutUp ---
+ call mct_avect_zero(avsrc_eroutUp)
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)
+ avsrc_eroutUp%rAttr(nt,cnt) = TRunoff%erout(iunit,nt)
enddo
enddo
call mct_avect_zero(avdst_eroutUp)
call mct_sMat_avMult(avsrc_eroutUp, sMatP_eroutUp, avdst_eroutUp)
- !--- add mapped eroutUp to Trunoff ---
+ !--- add mapped eroutUp to TRunoff ---
cnt = 0
do iunit = rtmCTL%begr,rtmCTL%endr
cnt = cnt + 1
do nt = 1,nt_rtm
- Trunoff%eroutUp(iunit,nt) = avdst_eroutUp%rAttr(nt,cnt)
+ TRunoff%eroutUp(iunit,nt) = avdst_eroutUp%rAttr(nt,cnt)
enddo
enddo
#endif
call t_stopf('mosart_SMeroutUp')
+ TRunoff%eroutup_avg = TRunoff%eroutup_avg + TRunoff%eroutUp
+ TRunoff%erlat_avg = TRunoff%erlat_avg + TRunoff%erlateral
+
!------------------
! channel routing
!------------------
@@ -171,7 +179,7 @@ subroutine Euler
endif
end do
end do
- negchan = min(negchan, minval(Trunoff%wr(:,:)))
+ negchan = min(negchan, minval(TRunoff%wr(:,:)))
call t_stopf('mosart_chanroute')
end do
@@ -182,6 +190,9 @@ subroutine Euler
! call shr_sys_abort('mosart: negative channel storage')
endif
TRunoff%flow = TRunoff%flow / Tctl%DLevelH2R
+ TRunoff%erout_prev = TRunoff%erout_prev / Tctl%DLevelH2R
+ TRunoff%eroutup_avg = TRunoff%eroutup_avg / Tctl%DLevelH2R
+ TRunoff%erlat_avg = TRunoff%erlat_avg / Tctl%DLevelH2R
end subroutine Euler
@@ -200,16 +211,7 @@ subroutine hillslopeRouting(iunit, nt, theDeltaT)
TRunoff%wh(iunit,nt) + (TRunoff%qsur(iunit,nt) + TRunoff%ehout(iunit,nt)) * theDeltaT < TINYVALUE) then
TRunoff%ehout(iunit,nt) = -(TRunoff%qsur(iunit,nt) + TRunoff%wh(iunit,nt) / theDeltaT)
end if
- TRunoff%dwh(iunit,nt) = (TRunoff%qsur(iunit,nt) + TRunoff%ehout(iunit,nt)) !* TUnit%area(iunit) * TUnit%frac(iunit)
- TRunoff%etin(iunit,nt) = (-TRunoff%ehout(iunit,nt) + TRunoff%qsub(iunit,nt)) * TUnit%area(iunit) * TUnit%frac(iunit)
- !if(TRunoff%etin(iunit,1) < 0) then
- ! call hillslopeRouting(iunit, Tctl%DeltaT)
- !end if
-
-! check stability
-! if(TRunoff%yh(iunit,nt) < -TINYVALUE .or. TRunoff%yh(iunit,nt) > 10) then
-! write(iulog,*) "Numerical error in hillslopeRouting, ", iunit,nt,TRunoff%yh(iunit,nt)
-! end if
+ TRunoff%dwh(iunit,nt) = (TRunoff%qsur(iunit,nt) + TRunoff%ehout(iunit,nt))
end subroutine hillslopeRouting
@@ -239,7 +241,7 @@ subroutine subnetworkRouting(iunit,nt,theDeltaT)
! check stability
! if(TRunoff%vt(iunit,nt) < -TINYVALUE .or. TRunoff%vt(iunit,nt) > 30) then
-! write(iulog,*) "Numerical error in subnetworkRouting, ", iunit,nt,Trunoff%vt(iunit,nt)
+! write(iulog,*) "Numerical error in subnetworkRouting, ", iunit,nt,TRunoff%vt(iunit,nt)
! end if
end subroutine subnetworkRouting
@@ -281,11 +283,11 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
TRunoff%erin(iunit,nt) = 0._r8
! tcraig, moved this out of the inner main channel loop to before main channel call
-! now it's precomputed as Trunoff%eroutUp
+! now it's precomputed as TRunoff%eroutUp
! do k=1,TUnit%nUp(iunit)
! TRunoff%erin(iunit,nt) = TRunoff%erin(iunit,nt) - TRunoff%erout(TUnit%iUp(iunit,k),nt)
! end do
- TRunoff%erin(iunit,nt) = TRunoff%erin(iunit,nt) - Trunoff%eroutUp(iunit,nt)
+ TRunoff%erin(iunit,nt) = TRunoff%erin(iunit,nt) - TRunoff%eroutUp(iunit,nt)
! estimate the outflow
if(TUnit%rlen(iunit) <= 0._r8) then ! no river network, no channel routing
@@ -311,9 +313,13 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
temp_gwl = TRunoff%qgwl(iunit,nt) * TUnit%area(iunit) * TUnit%frac(iunit)
temp_gwl0 = temp_gwl
if(abs(temp_gwl) <= TINYVALUE) then
+! write(iulog,*) 'mosart: ERROR dropping temp_gwl too small'
+! call shr_sys_abort('mosart: ERROR temp_gwl too small')
temp_gwl = 0._r8
end if
if(temp_gwl < -TINYVALUE) then
+ write(iulog,*) 'mosart: ERROR temp_gwl negative',iunit,nt,TRunoff%qgwl(iunit,nt)
+ call shr_sys_abort('mosart: ERROR temp_gwl negative ')
if(TRunoff%wr(iunit,nt) < TINYVALUE) then
temp_gwl = 0._r8
else
@@ -328,7 +334,7 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
! check for stability
! if(TRunoff%vr(iunit,nt) < -TINYVALUE .or. TRunoff%vr(iunit,nt) > 30) then
-! write(iulog,*) "Numerical error inRouting_KW, ", iunit,nt,Trunoff%vr(iunit,nt)
+! write(iulog,*) "Numerical error inRouting_KW, ", iunit,nt,TRunoff%vr(iunit,nt)
! end if
! check for negative wr
@@ -615,8 +621,13 @@ function GRPR(hr_, rwidth_, rwidth0_,rdepth_) result(pr_)
real(r8) :: SLOPE1 ! slope of flood plain, TO DO
real(r8) :: deltahr_
+ logical, save :: first_call = .true.
SLOPE1 = SLOPE1def
+ if (first_call) then
+ sinatanSLOPE1defr = 1.0_r8/(sin(atan(SLOPE1def)))
+ endif
+ first_call = .false.
if(hr_ < TINYVALUE) then
pr_ = 0._r8
diff --git a/models/rof/mosart/src/riverroute/RtmHistFile.F90 b/models/rof/mosart/src/riverroute/RtmHistFile.F90
index 41b91df5bcb1..a84b455508b6 100644
--- a/models/rof/mosart/src/riverroute/RtmHistFile.F90
+++ b/models/rof/mosart/src/riverroute/RtmHistFile.F90
@@ -809,6 +809,8 @@ 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))
else if (mode == 'write') then
@@ -837,6 +839,7 @@ 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')
endif
@@ -1178,7 +1181,7 @@ subroutine RtmHistRestart (ncid, flag, rdate)
! Create the restart history filename and open it
!
write(hnum,'(i1.1)') t-1
- locfnhr(t) = "./" // trim(caseid) //".rtm"// trim(inst_suffix) &
+ locfnhr(t) = "./" // trim(caseid) //".mosart"// trim(inst_suffix) &
// ".rh" // hnum //"."// trim(rdate) //".nc"
call htape_create( t, histrest=.true. )
!
@@ -1617,7 +1620,7 @@ character(len=256) function set_hist_filename (hist_freq, rtmhist_mfilt, hist_fi
write(cdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr,mon,day,sec
endif
write(hist_index,'(i1.1)') hist_file - 1
- set_hist_filename = "./"//trim(caseid)//".rtm"//trim(inst_suffix)//&
+ set_hist_filename = "./"//trim(caseid)//".mosart"//trim(inst_suffix)//&
".h"//hist_index//"."//trim(cdate)//".nc"
end function set_hist_filename
@@ -1689,21 +1692,21 @@ subroutine RtmHistAddfld (fname, units, avgflag, long_name, ptr_rof, default)
if (present(default)) then
if (trim(default) == 'inactive') return
- else
- ! Look through master list for input field name.
- ! When found, set active flag for that tape to true.
- found = .false.
- do f = 1,nfmaster
- if (trim(fname) == trim(masterlist(f)%field%name)) then
- masterlist(f)%actflag(1) = .true.
- found = .true.
- exit
- end if
- end do
- if (.not. found) then
- write(iulog,*) trim(subname),' ERROR: field=', fname, ' not found'
- call shr_sys_abort()
+ endif
+
+ ! Look through master list for input field name.
+ ! When found, set active flag for that tape to true.
+ found = .false.
+ do f = 1,nfmaster
+ if (trim(fname) == trim(masterlist(f)%field%name)) then
+ masterlist(f)%actflag(1) = .true.
+ found = .true.
+ exit
end if
+ end do
+ if (.not. found) then
+ write(iulog,*) trim(subname),' ERROR: field=', fname, ' not found'
+ call shr_sys_abort()
end if
end subroutine RtmHistAddfld
diff --git a/models/rof/mosart/src/riverroute/RtmHistFlds.F90 b/models/rof/mosart/src/riverroute/RtmHistFlds.F90
index f8ba3590e6e8..8867e66f2df5 100644
--- a/models/rof/mosart/src/riverroute/RtmHistFlds.F90
+++ b/models/rof/mosart/src/riverroute/RtmHistFlds.F90
@@ -38,61 +38,77 @@ subroutine RtmHistFldsInit()
implicit none
!-------------------------------------------------------
- call RtmHistAddfld (fname='QCHANR', units='m3/s', &
- avgflag='A', long_name='RTM river flow: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%runofflnd_nt1)
+ call RtmHistAddfld (fname='QCHANR'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ avgflag='A', long_name='MOSART river 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='RTM river flow: '//trim(rtm_tracers(2)), &
- ptr_rof=rtmCTL%runofflnd_nt2)
+ avgflag='A', long_name='MOSART river flow: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%runofflnd_nt2, default='active')
- call RtmHistAddfld (fname='QCHOCNR', units='m3/s', &
- avgflag='A', long_name='RTM river discharge into ocean: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%runoffocn_nt1)
+ call RtmHistAddfld (fname='QCHOCNR'//'_'//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', &
- avgflag='A', long_name='RTM river discharge into ocean: '//trim(rtm_tracers(2)), &
- ptr_rof=rtmCTL%runoffocn_nt2)
+ avgflag='A', long_name='MOSART river discharge into ocean: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%runoffocn_nt2, default='active')
- call RtmHistAddfld (fname='VOLR', units='m3', &
- avgflag='A', long_name='RTM storage: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%volr_nt1)
+ call RtmHistAddfld (fname='VOLR'//'_'//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', &
- avgflag='A', long_name='RTM storage: '//trim(rtm_tracers(2)), &
- ptr_rof=rtmCTL%volr_nt2, default='inactive')
+ avgflag='A', long_name='MOSART storage: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%volr_nt2, default='active')
- call RtmHistAddfld (fname='DVOLRDT_LND', units='mm/s', &
- avgflag='A', long_name='RTM land change in storage: '//trim(rtm_tracers(1)), &
+ 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')
- call RtmHistAddfld (fname='DVOLRDT_LND'//'_'//trim(rtm_tracers(2)), units='mm/s', &
- avgflag='A', long_name='RTM land change in storage: '//trim(rtm_tracers(2)), &
+ 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)), &
ptr_rof=rtmCTL%dvolrdtlnd_nt2, default='inactive')
- call RtmHistAddfld (fname='DVOLRDT_OCN', units='mm/s', &
- avgflag='A', long_name='RTM ocean change of storage: '//trim(rtm_tracers(1)), &
+ 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')
- call RtmHistAddfld (fname='DVOLRDT_OCN'//'_'//trim(rtm_tracers(2)), units='mm/s', &
- avgflag='A', long_name='RTM ocean change of storage: '//trim(rtm_tracers(2)), &
+ 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)), &
ptr_rof=rtmCTL%dvolrdtocn_nt2, default='inactive')
- call RtmHistAddfld (fname='QSUR'//'_'//trim(rtm_tracers(1)), units='mm/s', &
- avgflag='A', long_name='RTM input surface runoff: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%qsur, default='inactive')
+ 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='active')
- call RtmHistAddfld (fname='QSUB'//'_'//trim(rtm_tracers(1)), units='mm/s', &
- avgflag='A', long_name='RTM input subsurface runoff: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%qsub, default='inactive')
+ call RtmHistAddfld (fname='QSUR'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART input surface runoff: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%qsur_nt2, default='active')
- call RtmHistAddfld (fname='QGWL'//'_'//trim(rtm_tracers(1)), units='mm/s', &
- avgflag='A', long_name='RTM input GWL runoff: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%qgwl, default='inactive')
+ 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='active')
- call RtmHistAddfld (fname='QDTO'//'_'//trim(rtm_tracers(1)), units='mm/s', &
- avgflag='A', long_name='RTM input direct to ocean runoff: '//trim(rtm_tracers(1)), &
- ptr_rof=rtmCTL%qdto, default='inactive')
+ call RtmHistAddfld (fname='QSUB'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART input subsurface runoff: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%qsub_nt2, default='active')
+
+ 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='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='active')
+
+ call RtmHistAddfld (fname='QDTO'//'_'//trim(rtm_tracers(1)), units='m3/s', &
+ avgflag='A', long_name='MOSART input direct to ocean runoff: '//trim(rtm_tracers(1)), &
+ ptr_rof=rtmCTL%qdto_nt1, default='active')
+
+ call RtmHistAddfld (fname='QDTO'//'_'//trim(rtm_tracers(2)), units='m3/s', &
+ avgflag='A', long_name='MOSART input direct to ocean runoff: '//trim(rtm_tracers(2)), &
+ ptr_rof=rtmCTL%qdto_nt2, default='active')
! Print masterlist of history fields
@@ -106,7 +122,7 @@ subroutine RtmHistFldsSet()
!-----------------------------------------------------------------------
! !DESCRIPTION:
- ! Set rtm history fields as 1d poitner arrays
+ ! Set mosart history fields as 1d poitner arrays
!
implicit none
!-----------------------------------------------------------------------
@@ -125,8 +141,20 @@ subroutine RtmHistFldsSet()
rtmCTL%dvolrdtocn_nt1(:) = rtmCTL%dvolrdtocn(:,1)
rtmCTL%dvolrdtocn_nt2(:) = rtmCTL%dvolrdtocn(:,2)
- rtmCTL%volr_nt1(:) = rtmCTL%volrlnd(:,1)
- rtmCTL%volr_nt2(:) = rtmCTL%volrlnd(:,2)
+ rtmCTL%volr_nt1(:) = rtmCTL%volr(:,1)
+ rtmCTL%volr_nt2(:) = rtmCTL%volr(:,2)
+
+ rtmCTL%qsub_nt1(:) = rtmCTL%qsub(:,1)
+ rtmCTL%qsub_nt2(:) = rtmCTL%qsub(:,2)
+
+ rtmCTL%qsur_nt1(:) = rtmCTL%qsur(:,1)
+ rtmCTL%qsur_nt2(:) = rtmCTL%qsur(:,2)
+
+ rtmCTL%qgwl_nt1(:) = rtmCTL%qgwl(:,1)
+ rtmCTL%qgwl_nt2(:) = rtmCTL%qgwl(:,2)
+
+ rtmCTL%qdto_nt1(:) = rtmCTL%qdto(:,1)
+ rtmCTL%qdto_nt2(:) = rtmCTL%qdto(:,2)
end subroutine RtmHistFldsSet
diff --git a/models/rof/mosart/src/riverroute/RtmMod.F90 b/models/rof/mosart/src/riverroute/RtmMod.F90
index 9096f9bbe242..e0dea28ce7ee 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
@@ -73,6 +73,7 @@ module RtmMod
! MOSART constants
real(r8) :: cfl_scale = 1.0_r8 ! cfl scale factor, must be <= 1.0
+ real(r8) :: river_depth_minimum = 1.e-4 ! gridcell average minimum river depth [m]
!global (glo)
integer , pointer :: ID0_global(:) ! local ID index
@@ -82,10 +83,10 @@ module RtmMod
!local (gdc)
real(r8), save, pointer :: evel(:,:) ! effective tracer velocity (m/s)
- real(r8), save, pointer :: sfluxin(:,:) ! cell tracer influx (m3/s)
- real(r8), save, pointer :: fluxout(:,:) ! cell tracer outlflux (m3/s)
- real(r8), save, pointer :: direct(:,:) ! direct to outlet flow (m3/s)
real(r8), save, pointer :: flow(:,:) ! mosart flow (m3/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)
! global MOSART grid
real(r8),pointer :: rlatc(:) ! latitude of 1d grid cell (deg)
@@ -211,12 +212,12 @@ subroutine Rtmini(rtm_active,flood_active)
do_rtm = .true.
do_rtmflood = .false.
ice_runoff = .true.
+ wrmflag = .false.
finidat_rtm = ' '
nrevsn_rtm = ' '
coupling_period = -1
delt_mosart = 3600
decomp_option = 'basin'
- wrmflag = .false.
smat_option = 'opt'
nlfilename_rof = "mosart_in" // trim(inst_suffix)
@@ -933,20 +934,20 @@ subroutine Rtmini(rtm_active,flood_active)
call t_startf('mosarti_vars')
- allocate (fluxout (rtmCTL%begr:rtmCTL%endr,nt_rtm), &
- evel (rtmCTL%begr:rtmCTL%endr,nt_rtm), &
- direct (rtmCTL%begr:rtmCTL%endr,nt_rtm), &
+ allocate (evel (rtmCTL%begr:rtmCTL%endr,nt_rtm), &
flow (rtmCTL%begr:rtmCTL%endr,nt_rtm), &
- sfluxin (rtmCTL%begr:rtmCTL%endr,nt_rtm), stat=ier)
+ erout_prev(rtmCTL%begr:rtmCTL%endr,nt_rtm), &
+ eroutup_avg(rtmCTL%begr:rtmCTL%endr,nt_rtm), &
+ erlat_avg(rtmCTL%begr:rtmCTL%endr,nt_rtm), &
+ stat=ier)
if (ier /= 0) then
- write(iulog,*) subname,' Allocation ERROR for ',&
- 'fluxout'
- call shr_sys_abort(subname//' Allocationt ERROR volr')
+ write(iulog,*) subname,' Allocation ERROR for flow'
+ call shr_sys_abort(subname//' Allocationt ERROR flow')
end if
- fluxout(:,:) = 0._r8
- sfluxin(:,:) = 0._r8
- direct(:,:) = 0._r8
flow(:,:) = 0._r8
+ erout_prev(:,:) = 0._r8
+ eroutup_avg(:,:) = 0._r8
+ erlat_avg(:,:) = 0._r8
!-------------------------------------------------------
! Allocate runoff datatype
@@ -990,8 +991,7 @@ subroutine Rtmini(rtm_active,flood_active)
nr = rglo2gdc(n)
if (nr > 0) then
numr = numr + 1
- rgdc2glo(nr) = n
- rtmCTL%mask(nr) = gmask(n)
+ rgdc2glo(nr) = n
endif
end do
end do
@@ -999,17 +999,13 @@ subroutine Rtmini(rtm_active,flood_active)
write(iulog,*) subname,'ERROR numr and rtmCTL%numr are different ',numr,rtmCTL%numr
call shr_sys_abort(subname//' ERROR numr')
endif
- if (minval(rtmCTL%mask) < 1) then
- write(iulog,*) subname,'ERROR rtmCTL mask lt 1 ',minval(rtmCTL%mask),maxval(rtmCTL%mask)
- call shr_sys_abort(subname//' ERROR rtmCTL mask')
- endif
- deallocate(gmask)
! Determine runoff datatype variables
lrtmarea = 0.0_r8
cnt = 0
do nr = rtmCTL%begr,rtmCTL%endr
rtmCTL%gindex(nr) = rgdc2glo(nr)
+ rtmCTL%mask(nr) = gmask(rgdc2glo(nr))
n = rgdc2glo(nr)
i = mod(n-1,rtmlon) + 1
j = (n-1)/rtmlon + 1
@@ -1024,7 +1020,6 @@ subroutine Rtmini(rtm_active,flood_active)
do nt = 1,nt_rtm
rtmCTL%runofflnd(nr,nt) = rtmCTL%runoff(nr,nt)
rtmCTL%dvolrdtlnd(nr,nt)= rtmCTL%dvolrdt(nr,nt)
- rtmCTL%volrlnd(nr,nt)= rtmCTL%volr(nr,nt)
enddo
elseif (rtmCTL%mask(nr) >= 2) then
do nt = 1,nt_rtm
@@ -1048,6 +1043,7 @@ subroutine Rtmini(rtm_active,flood_active)
rtmCTL%dsig(nr) = dnID_global(n)
endif
enddo
+ deallocate(gmask)
deallocate(rglo2gdc)
deallocate(rgdc2glo)
deallocate (dnID_global,area_global)
@@ -1055,6 +1051,11 @@ subroutine Rtmini(rtm_active,flood_active)
call shr_mpi_sum(lrtmarea,rtmCTL%totarea,mpicom_rof,'mosart totarea',all=.true.)
if (masterproc) write(iulog,*) subname,' earth area ',4.0_r8*shr_const_pi*1.0e6_r8*re*re
if (masterproc) write(iulog,*) subname,' MOSART area ',rtmCTL%totarea
+ if (minval(rtmCTL%mask) < 1) then
+ write(iulog,*) subname,'ERROR rtmCTL mask lt 1 ',minval(rtmCTL%mask),maxval(rtmCTL%mask)
+ call shr_sys_abort(subname//' ERROR rtmCTL mask')
+ endif
+
!-------------------------------------------------------
! Compute Sparse Matrix for downstream advection
@@ -1335,27 +1336,31 @@ subroutine Rtmini(rtm_active,flood_active)
if ((nsrest == nsrStartup .and. finidat_rtm /= ' ') .or. &
(nsrest == nsrContinue) .or. &
(nsrest == nsrBranch )) then
- call RtmRestFileRead( file=fnamer )
- fluxout(:,:) = rtmCTL%fluxout(:,:)
- !write(iulog,*) ' MOSART init file is read'
- TRunoff%wh = rtmCTL%wh
- TRunoff%wt = rtmCTL%wt
- TRunoff%wr = rtmCTL%wr
- TRunoff%erout= rtmCTL%erout
- do nt = 1,nt_rtm
- do nr = rtmCTL%begr,rtmCTL%endr
- call UpdateState_hillslope(nr,nt)
- call UpdateState_subnetwork(nr,nt)
- call UpdateState_mainchannel(nr,nt)
- enddo
- enddo
- do nt = 1,nt_rtm
- do nr = rtmCTL%begr,rtmCTL%endr
- rtmCTL%volr(nr,nt) = (Trunoff%wt(nr,nt) + Trunoff%wr(nr,nt) + &
- Trunoff%wh(nr,nt)*rtmCTL%area(nr))
- enddo
- enddo
- end if
+ call RtmRestFileRead( file=fnamer )
+ !write(iulog,*) ' MOSART init file is read'
+ TRunoff%wh = rtmCTL%wh
+ TRunoff%wt = rtmCTL%wt
+ TRunoff%wr = rtmCTL%wr
+ TRunoff%erout= rtmCTL%erout
+ else
+! do nt = 1,nt_rtm
+! do nr = rtmCTL%begr,rtmCTL%endr
+! TRunoff%wh(nr,nt) = rtmCTL%area(nr) * river_depth_minimum * 1.e-10_r8
+! TRunoff%wt(nr,nt) = rtmCTL%area(nr) * river_depth_minimum * 1.e-8_r8
+! TRunoff%wr(nr,nt) = rtmCTL%area(nr) * river_depth_minimum * 10._r8
+! enddo
+! enddo
+ endif
+
+ do nt = 1,nt_rtm
+ do nr = rtmCTL%begr,rtmCTL%endr
+ call UpdateState_hillslope(nr,nt)
+ call UpdateState_subnetwork(nr,nt)
+ call UpdateState_mainchannel(nr,nt)
+ rtmCTL%volr(nr,nt) = (TRunoff%wt(nr,nt) + TRunoff%wr(nr,nt) + &
+ TRunoff%wh(nr,nt)*rtmCTL%area(nr))
+ enddo
+ enddo
call t_stopf('mosarti_restart')
@@ -1387,7 +1392,7 @@ end subroutine Rtmini
! !IROUTINE: Rtmrun
!
! !INTERFACE:
- subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate)
+ subroutine Rtmrun(rstwr,nlend,rdate)
!
! !DESCRIPTION:
! River routing model
@@ -1396,11 +1401,6 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
!
! !ARGUMENTS:
implicit none
- real(r8), pointer :: totrunin(:,:) ! total lnd forcing on mosart grid (mm/s)
- real(r8), pointer :: surrunin(:,:) ! surface forcing on mosart grid (mm/s)
- real(r8), pointer :: subrunin(:,:) ! subsurface forcing on mosart grid (mm/s)
- real(r8), pointer :: gwlrunin(:,:) ! glacier/wetland/lake forcing on mosart grid (mm/s)
- real(r8), pointer :: dtorunin(:,:) ! direct-to-ocean forcing on mosart grid (mm/s)
logical , intent(in) :: rstwr ! true => write restart file this step)
logical , intent(in) :: nlend ! true => end of run on this step
character(len=*), intent(in) :: rdate ! restart file time stamp for name
@@ -1415,11 +1415,12 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
! !LOCAL VARIABLES:
!EOP
integer :: i, j, n, nr, ns, nt, n2, nf ! indices
- real(r8) :: budget_terms(20,nt_rtm) ! BUDGET terms
- ! BUDGET terms 1 and 2 are initial and final volumes in m3
- ! BUDGET terms 3,4,5,6,7,8,9,11 are sur,sub,gwl,dto,tot,flow,flood,direct in m3/s
- real(r8) :: budget_input, budget_output, budget_volume, budget_total ! BUDGET terms
- real(r8) :: budget_global(20,nt_rtm) ! global budget sum
+ 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
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
@@ -1436,6 +1437,11 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
integer :: cnt ! counter for gridcells
integer :: ier ! error code
integer,parameter :: dbug = 1 ! local debug flag
+!scs
+! parameters used in negative runoff partitioning algorithm
+ real(r8) :: river_volume_minimum ! gridcell area multiplied by average river_depth_minimum [m3]
+ real(r8) :: qgwl_volume ! volume of runoff during time step [m3]
+!scs
character(len=*),parameter :: subname = '(Rtmrun) '
!-----------------------------------------------------------------------
@@ -1460,7 +1466,12 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
budget_terms = 0._r8
flow = 0._r8
- rtmCTL%runoffall = 0._r8
+ erout_prev = 0._r8
+ eroutup_avg = 0._r8
+ erlat_avg = 0._r8
+ rtmCTL%runoff = 0._r8
+ rtmCTL%direct = 0._r8
+ rtmCTL%flood = 0._r8
rtmCTL%runofflnd = spval
rtmCTL%runoffocn = spval
rtmCTL%dvolrdt = 0._r8
@@ -1468,29 +1479,33 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
rtmCTL%dvolrdtocn = spval
! BUDGET
- ! BUDGET terms 1 and 2 are initial and final volumes in m3
- ! BUDGET terms 3,4,5,6,7,8,9,11 are sur,sub,gwl,dto,tot,flow,flood,direct in m3/s
+ ! BUDGET terms 1-10 are for volumes (m3)
+ ! BUDGET terms 11-30 are for flows (m3/s)
if (budget_check) then
call t_startf('mosartr_budget')
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
- budget_terms(1,nt) = budget_terms(1,nt) + rtmCTL%volr(nr,nt)
- budget_terms(3,nt) = budget_terms(3,nt) + surrunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
- budget_terms(4,nt) = budget_terms(4,nt) + subrunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
- budget_terms(5,nt) = budget_terms(5,nt) + gwlrunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
- budget_terms(6,nt) = budget_terms(6,nt) + dtorunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
- budget_terms(7,nt) = budget_terms(7,nt) + totrunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
+ budget_terms( 1,nt) = budget_terms( 1,nt) + rtmCTL%volr(nr,nt)
+ budget_terms( 3,nt) = budget_terms( 3,nt) + TRunoff%wt(nr,nt)
+ budget_terms( 5,nt) = budget_terms( 5,nt) + TRunoff%wr(nr,nt)
+ budget_terms( 7,nt) = budget_terms( 7,nt) + TRunoff%wh(nr,nt)*rtmCTL%area(nr)
+ budget_terms(13,nt) = budget_terms(13,nt) + rtmCTL%qsur(nr,nt)
+ budget_terms(14,nt) = budget_terms(14,nt) + rtmCTL%qsub(nr,nt)
+ budget_terms(15,nt) = budget_terms(15,nt) + rtmCTL%qgwl(nr,nt)
+ budget_terms(16,nt) = budget_terms(16,nt) + rtmCTL%qdto(nr,nt)
+ budget_terms(17,nt) = budget_terms(17,nt) + rtmCTL%qsur(nr,nt) + rtmCTL%qsub(nr,nt) + &
+ rtmCTL%qgwl(nr,nt) + rtmCTL%qdto(nr,nt)
enddo
enddo
call t_stopf('mosartr_budget')
endif
- ! convert from kg/m2-s (mm/s) to m/s
+ ! data for euler solver, in m3/s here
do nr = rtmCTL%begr,rtmCTL%endr
do nt = 1,nt_rtm
- TRunoff%qsur(nr,nt) = surrunin(nr,nt) * 0.001_r8
- TRunoff%qsub(nr,nt) = subrunin(nr,nt) * 0.001_r8
- TRunoff%qgwl(nr,nt) = gwlrunin(nr,nt) * 0.001_r8
+ TRunoff%qsur(nr,nt) = rtmCTL%qsur(nr,nt)
+ TRunoff%qsub(nr,nt) = rtmCTL%qsub(nr,nt)
+ TRunoff%qgwl(nr,nt) = rtmCTL%qgwl(nr,nt)
enddo
enddo
@@ -1522,15 +1537,15 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
! it at the end or even during the run loop as the
! new volume is computed. fluxout depends on volr, so
! how this is implemented does impact the solution.
- Trunoff%qsur(nr,nt) = Trunoff%qsur(nr,nt) - (rtmCTL%flood(nr)/rtmCTL%area(nr))
+ TRunoff%qsur(nr,nt) = TRunoff%qsur(nr,nt) - rtmCTL%flood(nr)
endif
endif
enddo
call t_stopf('mosartr_flood')
!-----------------------------------
- ! DIRECT transfer to outlet point using sMat
- ! Remember to subract water from Trunoff forcing
+ ! DIRECT sMAT transfer to outlet point using sMat
+ ! Remember to subract water from TRunoff forcing
!-----------------------------------
if (barrier_timers) then
@@ -1548,10 +1563,8 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
do nr = rtmCTL%begr,rtmCTL%endr
cnt = cnt + 1
do nt = 1,nt_rtm
-! if (Trunoff%qgwl(nr,nt) < 0._r8) then
-! avsrc_direct%rAttr(nt,cnt) = TRunoff%qgwl(nr,nt) * TUnit%area(nr) * 0.001_r8
-! TRunoff%qgwl(nr,nt) = 0._r8
-! endif
+ avsrc_direct%rAttr(nt,cnt) = TRunoff%qgwl(nr,nt)
+ TRunoff%qgwl(nr,nt) = 0._r8
enddo
enddo
call mct_avect_zero(avdst_direct)
@@ -1560,21 +1573,64 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
!--- copy direct transfer water from AV to output field ---
cnt = 0
- direct(:,:) = 0._r8
do nr = rtmCTL%begr,rtmCTL%endr
cnt = cnt + 1
do nt = 1,nt_rtm
- direct(nr,nt) = avdst_direct%rAttr(nt,cnt)
+ rtmCTL%direct(nr,nt) = 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
- direct(nr,nt) = direct(nr,nt) + &
- dtorunin(nr,nt)*rtmCTL%area(nr)*0.001_r8
+ ! --- 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
@@ -1599,11 +1655,32 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
delt_save = delt
Tctl%DeltaT = delt
- do ns = 1,nsub
+ !-----------------------------------
+ ! mosart euler solver
+ ! --- convert TRunoff fields from m3/s to m/s before calling Euler
+ !-----------------------------------
- !-----------------------------------
- ! mosart euler solver
- !-----------------------------------
+ if (budget_check) then
+ call t_startf('mosartr_budget')
+ do nt = 1,nt_rtm
+ do nr = rtmCTL%begr,rtmCTL%endr
+ budget_terms(20,nt) = budget_terms(20,nt) + TRunoff%qsur(nr,nt) + &
+ TRunoff%qsub(nr,nt) + TRunoff%qgwl(nr,nt)
+ budget_terms(29,nt) = budget_terms(29,nt) + TRunoff%qgwl(nr,nt)
+ enddo
+ enddo
+ call t_stopf('mosartr_budget')
+ endif
+
+ do nt = 1,nt_rtm
+ do nr = rtmCTL%begr,rtmCTL%endr
+ 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
+ enddo
+
+ do ns = 1,nsub
if (wrmflag) then
#ifdef INCLUDE_WRM
@@ -1663,6 +1740,9 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
flow(nr,nt) = flow(nr,nt) + TRunoff%flow(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)
enddo
enddo
@@ -1672,7 +1752,10 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
! average flow over subcycling
!-----------------------------------
- flow = flow / float(nsub)
+ flow = flow / float(nsub)
+ erout_prev = erout_prev / float(nsub)
+ eroutup_avg = eroutup_avg / float(nsub)
+ erlat_avg = erlat_avg / float(nsub)
!-----------------------------------
! update states when subsycling completed
@@ -1686,15 +1769,16 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
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%dvolrdt = (rtmCTL%volr(nr,nt) - volr_init) / delt_coupling
- rtmCTL%runoffall(nr,nt) = flow(nr,nt) + direct(nr,nt)
- if (rtmCTL%mask(n) == 1) then
- rtmCTL%runofflnd(nr,nt) = rtmCTL%runoffall(nr,nt)
+ rtmCTL%volr(nr,nt) = (TRunoff%wt(nr,nt) + TRunoff%wr(nr,nt) + &
+ TRunoff%wh(nr,nt)*rtmCTL%area(nr))
+ rtmCTL%dvolrdt(nr,nt) = (rtmCTL%volr(nr,nt) - volr_init) / delt_coupling
+ rtmCTL%runoff(nr,nt) = flow(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(n) >= 2) then
- rtmCTL%runoffocn(nr,nt) = rtmCTL%runoffall(nr,nt)
+ elseif (rtmCTL%mask(nr) >= 2) then
+ rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt)
rtmCTL%dvolrdtocn(nr,nt)= rtmCTL%dvolrdt(nr,nt)
endif
enddo
@@ -1707,27 +1791,41 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
!-----------------------------------
! BUDGET
- ! BUDGET terms 1 and 2 are initial and final volumes in m3
- ! BUDGET terms 3,4,5,6,7,8,9,11 are sur,sub,gwl,dto,tot,flow,flood,direct in m3/s
+ ! 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
call t_startf('mosartr_budget')
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
budget_terms( 2,nt) = budget_terms( 2,nt) + rtmCTL%volr(nr,nt)
- budget_terms( 8,nt) = budget_terms( 8,nt) + flow(nr,nt)
- budget_terms(11,nt) = budget_terms(11,nt) + direct(nr,nt)
- budget_terms(12,nt) = budget_terms(12,nt) + rtmCTL%runoffall(nr,nt)
+ budget_terms( 4,nt) = budget_terms( 4,nt) + TRunoff%wt(nr,nt)
+ budget_terms( 6,nt) = budget_terms( 6,nt) + TRunoff%wr(nr,nt)
+ budget_terms( 8,nt) = budget_terms( 8,nt) + TRunoff%wh(nr,nt)*rtmCTL%area(nr)
+ budget_terms(21,nt) = budget_terms(21,nt) + rtmCTL%direct(nr,nt)
+ if (rtmCTL%mask(nr) >= 2) then
+ budget_terms(18,nt) = budget_terms(18,nt) + rtmCTL%runoff(nr,nt)
+ budget_terms(26,nt) = budget_terms(26,nt) - erout_prev(nr,nt)
+ budget_terms(27,nt) = budget_terms(27,nt) + flow(nr,nt)
+ else
+ budget_terms(23,nt) = budget_terms(23,nt) - erout_prev(nr,nt)
+ budget_terms(24,nt) = budget_terms(24,nt) + flow(nr,nt)
+ endif
+ budget_terms(25,nt) = budget_terms(25,nt) - eroutup_avg(nr,nt)
+ budget_terms(28,nt) = budget_terms(28,nt) - erlat_avg(nr,nt)
+ budget_terms(22,nt) = budget_terms(22,nt) + rtmCTL%runoff(nr,nt) + rtmCTL%direct(nr,nt) + eroutup_avg(nr,nt)
enddo
enddo
nt = 1
do nr = rtmCTL%begr,rtmCTL%endr
- budget_terms( 9,nt) = budget_terms( 9,nt) + rtmCTL%flood(nr)
+ budget_terms(19,nt) = budget_terms(19,nt) + rtmCTL%flood(nr)
+ budget_terms(22,nt) = budget_terms(22,nt) + rtmCTL%flood(nr)
enddo
!--- check budget
! convert fluxes from m3/s to m3 by mult by coupling_period
- budget_terms(3:20,:) = budget_terms(3:20,:) * delt_coupling
+ budget_terms(11:30,:) = budget_terms(11:30,:) * delt_coupling
! convert terms from m3 to million m3
budget_terms(:,:) = budget_terms(:,:) * 1.0e-6_r8
@@ -1740,31 +1838,60 @@ subroutine Rtmrun(totrunin,surrunin,subrunin,gwlrunin,dtorunin,rstwr,nlend,rdate
write(iulog,'(2a,i10,i6)') trim(subname),' MOSART BUDGET diagnostics (million m3) for ',ymd,tod
do nt = 1,nt_rtm
budget_volume = (budget_global( 2,nt) - budget_global( 1,nt))
- budget_input = (budget_global( 3,nt) + budget_global( 4,nt) + &
- budget_global( 5,nt) + budget_global( 6,nt))
- budget_output = (budget_global( 8,nt) + budget_global( 9,nt) + &
- budget_global(11,nt))
+ budget_input = (budget_global(13,nt) + budget_global(14,nt) + &
+ budget_global(15,nt) + budget_global(16,nt))
+ budget_output = (budget_global(18,nt) + budget_global(19,nt) + &
+ budget_global(21,nt))
budget_total = budget_volume - budget_input + budget_output
+ budget_euler = budget_volume - budget_global(20,nt) + budget_global(18,nt)
+ budget_eroutlag = budget_global(23,nt) - budget_global(24,nt)
write(iulog,'(2a,i4)') trim(subname),' tracer = ',nt
write(iulog,'(2a,i4,f22.6)') trim(subname),' volume init = ',nt,budget_global(1,nt)
write(iulog,'(2a,i4,f22.6)') trim(subname),' volume final = ',nt,budget_global(2,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input surface = ',nt,budget_global(3,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input subsurf = ',nt,budget_global(4,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input gwl = ',nt,budget_global(5,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input dto = ',nt,budget_global(6,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input total = ',nt,budget_global(7,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' input check = ',nt,budget_input - budget_global(7,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' output flow = ',nt,budget_global(8,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' output direct = ',nt,budget_global(11,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' output rofall = ',nt,budget_global(12,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' output flood = ',nt,budget_global(9,nt)
- write(iulog,'(2a,i4,f22.6)') trim(subname),' output check = ',nt,budget_output - budget_global(12,nt) - budget_global(9,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumeh init = ',nt,budget_global(7,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumeh final = ',nt,budget_global(8,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumet init = ',nt,budget_global(3,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumet final = ',nt,budget_global(4,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumer init = ',nt,budget_global(5,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' volumer final = ',nt,budget_global(6,nt)
+ !write(iulog,'(2a)') trim(subname),'----------------'
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' input surface = ',nt,budget_global(13,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' input subsurf = ',nt,budget_global(14,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' input gwl = ',nt,budget_global(15,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' input dto = ',nt,budget_global(16,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' input total = ',nt,budget_global(17,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' input check = ',nt,budget_input - budget_global(17,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' input euler = ',nt,budget_global(20,nt)
+ !write(iulog,'(2a)') trim(subname),'----------------'
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' output flow = ',nt,budget_global(18,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' output direct = ',nt,budget_global(21,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' output flood = ',nt,budget_global(19,nt)
+ write(iulog,'(2a,i4,f22.6)') trim(subname),' output total = ',nt,budget_global(22,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' output check = ',nt,budget_output - budget_global(22,nt)
!write(iulog,'(2a)') trim(subname),'----------------'
write(iulog,'(2a,i4,f22.6)') trim(subname),' sum input = ',nt,budget_input
write(iulog,'(2a,i4,f22.6)') trim(subname),' sum dvolume = ',nt,budget_volume
write(iulog,'(2a,i4,f22.6)') trim(subname),' sum output = ',nt,budget_output
- !write(iulog,'(a,16x,a)') trim(subname),'----------------'
+ !write(iulog,'(2a)') trim(subname),'----------------'
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)') 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)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' eroutup_avg = ',nt,budget_global(25,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' erout_prev out= ',nt,budget_global(26,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' erout out= ',nt,budget_global(27,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' erlateral = ',nt,budget_global(28,nt)
+ !write(iulog,'(2a,i4,f22.6)') trim(subname),' euler gwl = ',nt,budget_global(29,nt)
+ !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
+ endif
+ endif
enddo
write(iulog,'(a)') '----------------------------------- '
endif
@@ -1978,8 +2105,16 @@ subroutine MOSART_init
Tunit%mask(n) = 0
elseif (Tunit%mask(n) == 0) then
Tunit%mask(n) = 2
+ if (abs(Tunit%frac(n)-1.0_r8)>1.0e-9) then
+ write(iulog,*) subname,' ERROR frac ne 1.0',n,Tunit%frac(n)
+ call shr_sys_abort(subname//' ERROR frac ne 1.0')
+ endif
elseif (Tunit%mask(n) > 0) then
Tunit%mask(n) = 1
+ if (abs(Tunit%frac(n)-1.0_r8)>1.0e-9) then
+ write(iulog,*) subname,' ERROR frac ne 1.0',n,Tunit%frac(n)
+ call shr_sys_abort(subname//' ERROR frac ne 1.0')
+ endif
else
call shr_sys_abort(subname//' Tunit mask error')
endif
@@ -2019,6 +2154,14 @@ subroutine MOSART_init
if (masterproc) write(iulog,FORMR) trim(subname),' read area ',minval(Tunit%area),maxval(Tunit%area)
call shr_sys_flush(iulog)
+ do n=rtmCtl%begr, rtmCTL%endr
+ if (TUnit%area(n) < 0._r8) TUnit%area(n) = rtmCTL%area(n)
+ if (TUnit%area(n) /= rtmCTL%area(n)) then
+ write(iulog,*) subname,' ERROR area mismatch',TUnit%area(n),rtmCTL%area(n)
+ call shr_sys_abort(subname//' ERROR area mismatch')
+ endif
+ enddo
+
allocate(TUnit%areaTotal(begr:endr))
ier = pio_inq_varid(ncid, name='areaTotal', vardesc=vardesc)
call pio_read_darray(ncid, vardesc, iodesc_dbl, TUnit%areaTotal, ier)
@@ -2218,9 +2361,18 @@ subroutine MOSART_init
allocate (TRunoff%erout(begr:endr,nt_rtm))
TRunoff%erout = 0._r8
+ allocate (TRunoff%erout_prev(begr:endr,nt_rtm))
+ TRunoff%erout_prev = 0._r8
+
allocate (TRunoff%eroutUp(begr:endr,nt_rtm))
TRunoff%eroutUp = 0._r8
+ allocate (TRunoff%eroutUp_avg(begr:endr,nt_rtm))
+ TRunoff%eroutUp_avg = 0._r8
+
+ allocate (TRunoff%erlat_avg(begr:endr,nt_rtm))
+ TRunoff%erlat_avg = 0._r8
+
allocate (TRunoff%ergwl(begr:endr,nt_rtm))
TRunoff%ergwl = 0._r8
diff --git a/models/rof/mosart/src/riverroute/RtmRestFile.F90 b/models/rof/mosart/src/riverroute/RtmRestFile.F90
index 54c5266f91eb..97795818b1da 100644
--- a/models/rof/mosart/src/riverroute/RtmRestFile.F90
+++ b/models/rof/mosart/src/riverroute/RtmRestFile.F90
@@ -211,7 +211,7 @@ subroutine RtmRestGetfile( file, path )
! Check case name consistency (case name must be different
! for branch run, unless brnch_retain_casename is set)
- ctest = 'xx.'//trim(caseid)//'.rtm'
+ ctest = 'xx.'//trim(caseid)//'.mosart'
ftest = 'xx.'//trim(file)
status = index(trim(ftest),trim(ctest))
if (status /= 0 .and. .not.(brnch_retain_casename)) then
@@ -308,7 +308,7 @@ character(len=256) function RtmRestFileName( rdate )
implicit none
character(len=*), intent(in) :: rdate ! input date for restart file name
- RtmRestFileName = "./"//trim(caseid)//".rtm"//trim(inst_suffix)//".r."//trim(rdate)//".nc"
+ RtmRestFileName = "./"//trim(caseid)//".mosart"//trim(inst_suffix)//".r."//trim(rdate)//".nc"
if (masterproc) then
write(iulog,*)'writing restart file ',trim(RtmRestFileName),' for model date = ',rdate
end if
@@ -381,7 +381,7 @@ subroutine RtmRestart(ncid, flag)
character(len=128) :: lname
!-----------------------------------------------------------------------
- do nv = 1,8
+ do nv = 1,7
do nt = 1,nt_rtm
if (nv == 1) then
@@ -390,36 +390,31 @@ subroutine RtmRestart(ncid, flag)
uname = 'm3'
dfld => rtmCTL%volr(:,nt)
elseif (nv == 2) then
- vname = 'RTM_FLUXOUT_'//trim(rtm_tracers(nt))
- lname = 'water fluxout in cell (fluxout)'
- uname = 'm3/s'
- dfld => rtmCTL%fluxout(:,nt)
- elseif (nv == 3) then
vname = 'RTM_RUNOFF_'//trim(rtm_tracers(nt))
lname = 'runoff (runoff)'
uname = 'm3/s'
dfld => rtmCTL%runoff(:,nt)
- elseif (nv == 4) then
+ elseif (nv == 3) then
vname = 'RTM_DVOLRDT_'//trim(rtm_tracers(nt))
lname = 'water volume change in cell (dvolrdt)'
uname = 'mm/s'
dfld => rtmCTL%dvolrdt(:,nt)
- elseif (nv == 5) then
+ elseif (nv == 4) then
vname = 'RTM_WH_'//trim(rtm_tracers(nt))
lname = 'surface water storage at hillslopes in cell'
uname = 'm'
dfld => rtmCTL%wh(:,nt)
- elseif (nv == 6) then
+ elseif (nv == 5) then
vname = 'RTM_WT_'//trim(rtm_tracers(nt))
lname = 'water storage in tributary channels in cell'
uname = 'm3'
dfld => rtmCTL%wt(:,nt)
- elseif (nv == 7) then
+ elseif (nv == 6) then
vname = 'RTM_WR_'//trim(rtm_tracers(nt))
lname = 'water storage in main channel in cell'
uname = 'm3'
dfld => rtmCTL%wr(:,nt)
- elseif (nv == 8) then
+ elseif (nv == 7) then
vname = 'RTM_EROUT_'//trim(rtm_tracers(nt))
lname = 'instataneous flow out of main channel in cell'
uname = 'm3/s'
@@ -454,7 +449,6 @@ subroutine RtmRestart(ncid, flag)
if (abs(rtmCTL%volr(n,nt)) > 1.e30) rtmCTL%volr(n,nt) = 0.
if (abs(rtmCTL%runoff(n,nt)) > 1.e30) rtmCTL%runoff(n,nt) = 0.
if (abs(rtmCTL%dvolrdt(n,nt)) > 1.e30) rtmCTL%dvolrdt(n,nt) = 0.
- if (abs(rtmCTL%fluxout(n,nt)) > 1.e30) rtmCTL%fluxout(n,nt) = 0.
if (abs(rtmCTL%wh(n,nt)) > 1.e30) rtmCTL%wh(n,nt) = 0.
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.
@@ -464,7 +458,6 @@ subroutine RtmRestart(ncid, flag)
do nt = 1,nt_rtm
rtmCTL%runofflnd(n,nt) = rtmCTL%runoff(n,nt)
rtmCTL%dvolrdtlnd(n,nt)= rtmCTL%dvolrdt(n,nt)
- rtmCTL%volrlnd(n,nt) = rtmCTL%volr(n,nt)
end do
elseif (rtmCTL%mask(n) >= 2) then
do nt = 1,nt_rtm
diff --git a/models/rof/mosart/src/riverroute/RunoffMod.F90 b/models/rof/mosart/src/riverroute/RunoffMod.F90
index b103ab9d9ce9..5ddebd00e99f 100644
--- a/models/rof/mosart/src/riverroute/RunoffMod.F90
+++ b/models/rof/mosart/src/riverroute/RunoffMod.F90
@@ -42,36 +42,44 @@ module RunoffMod
integer , pointer :: dsig(:) ! downstream index, global index
integer , pointer :: outletg(:) ! outlet index, global index
- ! - local runtime
- real(r8), pointer :: runoff(:,:) ! RTM flow (m3 H2O/s)
- real(r8), pointer :: runoffall(:,:) ! RTM flow (m3 H2O/s)
+ ! - global
+ integer , pointer :: mask(:) ! general mask of cell 1=land, 2=ocean, 3=outlet
+ real(r8), pointer :: rlon(:) ! rtm longitude list, 1d
+ real(r8), pointer :: rlat(:) ! rtm latitude list, 1d
+ real(r8) :: totarea ! global area
+ integer :: numr ! rtm gdc global number of cells
+
+ ! - local
+ integer :: begr,endr ! local start/stop indices
+ integer :: lnumr ! local number of cells
+
+ ! - 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 :: 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 :: volrlnd(:,:) ! RTM storage masked for land (m3)
- real(r8), pointer :: fluxout(:,:) ! RTM cell tracer outlflux (m3/s)
real(r8), pointer :: fthresh(:) ! RTM water flood threshold
- real(r8), pointer :: flood(:) ! RTM water (flood) sent back to clm (mm/s)
+
+ ! - restarts
real(r8), pointer :: wh(:,:) ! MOSART hillslope surface water storage (m)
real(r8), pointer :: wt(:,:) ! MOSART sub-network water storage (m3)
real(r8), pointer :: wr(:,:) ! MOSART main channel water storage (m3)
real(r8), pointer :: erout(:,:) ! MOSART flow out of the main channel, instantaneous (m3/s)
- ! - global
- integer , pointer :: mask(:) ! general mask of cell 1=land, 2=ocean, 3=outlet
- real(r8), pointer :: rlon(:) ! rtm longitude list, 1d
- real(r8), pointer :: rlat(:) ! rtm latitude list, 1d
- real(r8) :: totarea ! global area
- integer :: numr ! rtm gdc global number of cells
+ ! inputs
+ 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 :: qdto(:,:) ! coupler diret-to-ocean forcing [m3/s]
- ! - local
- integer :: begr,endr ! local start/stop indices
- integer :: lnumr ! local number of cells
+ ! - outputs
+ real(r8), pointer :: flood(:) ! coupler return flood water sent back to clm [m3/s]
+ real(r8), pointer :: runoff(:,:) ! coupler return mosart basin derived flow [m3/s]
+ real(r8), pointer :: direct(:,:) ! coupler return direct flow [m3/s]
- ! - 1d field pointers for history files (currently needed)
+ ! - history (currently needed)
real(r8), pointer :: runofflnd_nt1(:)
real(r8), pointer :: runofflnd_nt2(:)
real(r8), pointer :: runoffocn_nt1(:)
@@ -82,11 +90,14 @@ module RunoffMod
real(r8), pointer :: dvolrdtocn_nt2(:)
real(r8), pointer :: volr_nt1(:)
real(r8), pointer :: volr_nt2(:)
-
- real(r8), pointer :: qsur(:)
- real(r8), pointer :: qsub(:)
- real(r8), pointer :: qgwl(:)
- real(r8), pointer :: qdto(:)
+ real(r8), pointer :: qsur_nt1(:)
+ real(r8), pointer :: qsur_nt2(:)
+ real(r8), pointer :: qsub_nt1(:)
+ real(r8), pointer :: qsub_nt2(:)
+ real(r8), pointer :: qgwl_nt1(:)
+ real(r8), pointer :: qgwl_nt2(:)
+ real(r8), pointer :: qdto_nt1(:)
+ real(r8), pointer :: qdto_nt2(:)
end type runoff_flow
@@ -177,8 +188,8 @@ module RunoffMod
type TstatusFlux
! hillsloope
!! states
- real(r8), pointer :: wh(:,:) ! storage of surface water, [m3]
- real(r8), pointer :: dwh(:,:) ! change of water storage, [m3]
+ real(r8), pointer :: wh(:,:) ! storage of surface water, [m]
+ real(r8), pointer :: dwh(:,:) ! change of water storage, [m/s]
real(r8), pointer :: yh(:,:) ! depth of surface water, [m]
real(r8), pointer :: wsat(:,:) ! storage of surface water within saturated area at hillslope [m]
real(r8), pointer :: wunsat(:,:) ! storage of surface water within unsaturated area at hillslope [m]
@@ -225,7 +236,10 @@ module RunoffMod
real(r8), pointer :: erlateral(:,:) ! 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]
real(r8), pointer :: eroutUp(:,:) ! outflow sum of upstream gridcells, instantaneous (m3/s)
+ real(r8), pointer :: eroutUp_avg(:,:) ! outflow sum of upstream gridcells, average [m3/s]
+ real(r8), pointer :: erlat_avg(:,:) ! erlateral average [m3/s]
real(r8), pointer :: flow(:,:) ! streamflow from the outlet of the reach, [m3/s]
real(r8), pointer :: erin1(:,:) ! inflow from upstream links during previous step, used for Muskingum method, [m3/s]
real(r8), pointer :: erin2(:,:) ! inflow from upstream links during current step, used for Muskingum method, [m3/s]
@@ -270,7 +284,6 @@ subroutine RunoffInit(begr, endr, numr)
integer :: ier
allocate(rtmCTL%runoff(begr:endr,nt_rtm), &
- rtmCTL%runoffall(begr:endr,nt_rtm), &
rtmCTL%dvolrdt(begr:endr,nt_rtm), &
rtmCTL%runofflnd(begr:endr,nt_rtm), &
rtmCTL%dvolrdtlnd(begr:endr,nt_rtm), &
@@ -278,8 +291,6 @@ subroutine RunoffInit(begr, endr, numr)
rtmCTL%dvolrdtocn(begr:endr,nt_rtm), &
rtmCTL%area(begr:endr), &
rtmCTL%volr(begr:endr,nt_rtm), &
- rtmCTL%volrlnd(begr:endr,nt_rtm), &
- rtmCTL%fluxout(begr:endr,nt_rtm), &
rtmCTL%lonc(begr:endr), &
rtmCTL%latc(begr:endr), &
rtmCTL%dsig(begr:endr), &
@@ -294,18 +305,27 @@ subroutine RunoffInit(begr, endr, numr)
rtmCTL%dvolrdtlnd_nt2(begr:endr), &
rtmCTL%dvolrdtocn_nt1(begr:endr), &
rtmCTL%dvolrdtocn_nt2(begr:endr), &
- rtmCTL%mask(numr), &
+ rtmCTL%qsur_nt1(begr:endr), &
+ rtmCTL%qsur_nt2(begr:endr), &
+ rtmCTL%qsub_nt1(begr:endr), &
+ rtmCTL%qsub_nt2(begr:endr), &
+ rtmCTL%qgwl_nt1(begr:endr), &
+ rtmCTL%qgwl_nt2(begr:endr), &
+ rtmCTL%qdto_nt1(begr:endr), &
+ rtmCTL%qdto_nt2(begr:endr), &
+ rtmCTL%mask(begr:endr), &
rtmCTL%gindex(begr:endr), &
rtmCTL%fthresh(begr:endr), &
rtmCTL%flood(begr:endr), &
+ rtmCTL%direct(begr:endr,nt_rtm), &
rtmCTL%wh(begr:endr,nt_rtm), &
rtmCTL%wt(begr:endr,nt_rtm), &
rtmCTL%wr(begr:endr,nt_rtm), &
rtmCTL%erout(begr:endr,nt_rtm), &
- rtmCTL%qsur(begr:endr), &
- rtmCTL%qsub(begr:endr), &
- rtmCTL%qgwl(begr:endr), &
- rtmCTL%qdto(begr:endr), &
+ rtmCTL%qsur(begr:endr,nt_rtm), &
+ rtmCTL%qsub(begr:endr,nt_rtm), &
+ rtmCTL%qgwl(begr:endr,nt_rtm), &
+ rtmCTL%qdto(begr:endr,nt_rtm), &
stat=ier)
if (ier /= 0) then
write(iulog,*)'Rtmini ERROR allocation of runoff local arrays'
@@ -313,20 +333,19 @@ subroutine RunoffInit(begr, endr, numr)
end if
rtmCTL%runoff(:,:) = 0._r8
- rtmCTL%runoffall(:,:) = 0._r8
rtmCTL%runofflnd(:,:) = spval
rtmCTL%runoffocn(:,:) = spval
rtmCTL%dvolrdt(:,:) = 0._r8
rtmCTL%dvolrdtlnd(:,:) = spval
rtmCTL%dvolrdtocn(:,:) = spval
rtmCTL%volr(:,:) = 0._r8
- rtmCTL%volrlnd(:,:) = spval
rtmCTL%flood(:) = 0._r8
+ rtmCTL%direct(:,:) = 0._r8
- rtmCTL%qsur(:) = 0._r8
- rtmCTL%qsub(:) = 0._r8
- rtmCTL%qgwl(:) = 0._r8
- rtmCTL%qdto(:) = 0._r8
+ rtmCTL%qsur(:,:) = 0._r8
+ rtmCTL%qsub(:,:) = 0._r8
+ rtmCTL%qgwl(:,:) = 0._r8
+ rtmCTL%qdto(:,:) = 0._r8
end subroutine RunoffInit