From 5f404ff1d33f69b6a14b2afea53f5c065ad23a4e Mon Sep 17 00:00:00 2001 From: Anthony Craig Date: Sun, 29 Nov 2015 22:13:55 -0800 Subject: [PATCH] mosart code cleanup, add budget diagnostics, review history files *********1*********2*********3*********4*********5*********6*********7** Code cleanup to remove dead code and make the code clearer. More could be done. Update the direct output computation. Add budget diagnostics to verify conservation. In short tests, the model conserves except for in the frozen term in the euler solver. There is still a small issue there that needs further diagnosing. Review history and restart files and checked many fields. Updated the MOSART input file to adjust the rotated Antarctica. LG-20 [Non-BFB] --- .../namelist_defaults_mosart.xml | 2 +- models/rof/mosart/doc/ChangeLog | 23 ++ models/rof/mosart/src/cpl/rof_comp_esmf.F90 | 169 ++++---- models/rof/mosart/src/cpl/rof_comp_mct.F90 | 105 ++--- .../src/riverroute/MOSART_physics_mod.F90 | 61 +-- .../rof/mosart/src/riverroute/RtmHistFile.F90 | 35 +- .../rof/mosart/src/riverroute/RtmHistFlds.F90 | 104 +++-- models/rof/mosart/src/riverroute/RtmMod.F90 | 382 ++++++++++++------ .../rof/mosart/src/riverroute/RtmRestFile.F90 | 23 +- .../rof/mosart/src/riverroute/RunoffMod.F90 | 93 +++-- 10 files changed, 593 insertions(+), 404 deletions(-) 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