As a followup to #80 we do want in place a working version of FLOOD for MOSART. This will take some work to put into place. And we should have some tests for it when it is installed.
Here is what @swensosc suggests we should do, first remove the contents of the RtmFloodInit subroutine and then do the following...
index c597256..fe42015 100644
--- a/src/riverroute/RtmMod.F90
+++ b/src/riverroute/RtmMod.F90
@@ -974,24 +974,9 @@ subroutine Rtmini
call RunoffInit(rtmCTL%begr, rtmCTL%endr, rtmCTL%numr)
- !-------------------------------------------------------
- ! Initialize mosart flood - rtmCTL%fthresh and evel
- !-------------------------------------------------------
-
- if (do_rtmflood) then
- write(iulog,*) subname,' Flood not validated in this version, abort'
- call shr_sys_abort(subname//' Flood feature unavailable')
- call RtmFloodInit (frivinp_rtm, rtmCTL%begr, rtmCTL%endr, rtmCTL%fthresh, evel)
- else
- effvel(:) = effvel0 ! downstream velocity (m/s)
- rtmCTL%fthresh(:) = abs(spval)
- do nt = 1,nt_rtm
- do nr = rtmCTL%begr,rtmCTL%endr
- evel(nr,nt) = effvel(nt)
- enddo
- enddo
- end if
+!scs: I think all of the code, even the evel/effvel stuff can be removed; it appears to be local
+
!-------------------------------------------------------
! Initialize runoff data type
!-------------------------------------------------------
@@ -1420,6 +1405,12 @@ subroutine Rtmrun(rstwr,nlend,rdate)
real(r8) :: qgwl_volume ! volume of runoff during time step [m3]
real(r8) :: irrig_volume ! volume of irrigation demand during time step [m3]
+ real(r8) :: scalar
+ real(r8) :: bheight ! height above bankfull (m)
+ real(r8) :: swidth ! river surface width (m)
+ real(r8) :: rarea_cross ! river cross section area (m2)
+ real(r8) :: rarea_surface ! river surface area (m2)
+
character(len=*),parameter :: subname = '(Rtmrun) '
!-----------------------------------------------------------------------
@@ -1542,33 +1533,55 @@ subroutine Rtmrun(rstwr,nlend,rdate)
! rtmCTL%flood is m3/s here
!-----------------------------------
- call t_startf('mosartr_flood')
- nt = 1
- rtmCTL%flood = 0._r8
- do nr = rtmCTL%begr,rtmCTL%endr
- ! initialize rtmCTL%flood to zero
- if (rtmCTL%mask(nr) == 1) then
- if (rtmCTL%volr(nr,nt) > rtmCTL%fthresh(nr)) then
- ! determine flux that is sent back to the land
- ! this is in m3/s
- rtmCTL%flood(nr) = &
- (rtmCTL%volr(nr,nt)-rtmCTL%fthresh(nr)) / (delt_coupling)
+ !scs
+ if (do_rtmflood) then
+ ! determine flux that is sent back to the land, in m3/s
+
+ call t_startf('mosartr_flood')
+ nt = 1
+ rtmCTL%flood = 0._r8
+ do nr = rtmCTL%begr,rtmCTL%endr
+ ! initialize rtmCTL%flood to zero
+ if (rtmCTL%mask(nr) == 1) then
+ if (TRunoff%yr(nr,1) > TUnit%rdepth(nr)) then
+ ! compute flooding flux when water depth is greater than bankfull depth
+ scalar = 5e-2_r8
+
+ ! height of water surface above bankfull height
+ bheight = TRunoff%yr(nr,1) - TUnit%rdepth(nr)
+
+ if (bheight > 0._r8) then
+ ! river surface width
+ swidth = TUnit%rwidth(nr) + 2._r8*bheight/0.1_r8
+
+ ! river cross section area
+ rarea_cross = (TUnit%rwidth(nr) * TUnit%rdepth(nr)) + 0.5*bheight*(swidth + TUnit%rwidth(nr))
+
+ ! river surface area
+ rarea_surface = swidth * TRunoff%wr(nr,1) / TRunoff%mr(nr,1)
+
+ ! surface area of river sets scale for scalar
+ rtmCTL%flood(nr) = scalar * bheight * rarea_surface / delt_coupling
+
+ ! limit flood water to above bankfull water amount
+ if (TRunoff%wr(nr,1) - (rtmCTL%flood(nr) * delt_coupling) < (TRunoff%wr(nr,1) * TUnit%rwidth(nr) * TUnit%rdepth(nr) / TRunoff%mr(nr,1))) then
+ rtmCTL%flood(nr) = (TRunoff%wr(nr,1) - (TRunoff%wr(nr,1) * TUnit%rwidth(nr) * TUnit%rdepth(nr) / TRunoff%mr(nr,1)))/delt_coupling
+ endif
+
+ ! remove flood water from wr (main channel)
+ TRunoff%wr(nr,1) = TRunoff%wr(nr,1) - rtmCTL%flood(nr) * delt_coupling
+
+ if(TRunoff%wr(nr,1) < 0._r8) then
+ write(iulog,*) ' negative_wr ',TRunoff%wr(nr,1), rtmCTL%flood(nr) * delt_coupling
+ endif
+ endif
+ endif
- ! rtmCTL%flood will be sent back to land - so must subtract this
- ! from the input runoff from land
- ! tcraig, comment - this seems like an odd approach, you
- ! might create negative forcing. why not take it out of
- ! the volr directly? it's also odd to compute this
- ! at the initial time of the time loop. why not do
- ! 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)
endif
- endif
- enddo
- call t_stopf('mosartr_flood')
-
+ enddo
+ call t_stopf('mosartr_flood')
+ endif
+
!-----------------------------------------------------
! DIRECT sMAT transfer to outlet point using sMat
! Remember to subtract water from TRunoff forcing
@@ -1892,6 +1905,11 @@ subroutine Rtmrun(rstwr,nlend,rdate)
rtmCTL%wr = TRunoff%wr
rtmCTL%erout = TRunoff%erout
+ !scs: set river water height for history (just tracer 1 i.e. liquid)
+ rtmCTL%rheight = TRunoff%yr(:,1)
+ rtmCTL%fheight = TRunoff%yr(:,1) - TUnit%rdepth
+ rtmCTL%flood_flux = 1e3*rtmCTL%flood/rtmCTL%area ! mm/s
+
do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
volr_init = rtmCTL%volr(nr,nt)
@ekluzek This should not precede the outstanding PR's that are already open and awaiting disposition. Why hasn't #74 been merged? I feel that you have broken a commitment.
As a followup to #80 we do want in place a working version of FLOOD for MOSART. This will take some work to put into place. And we should have some tests for it when it is installed.
Here is what @swensosc suggests we should do, first remove the contents of the RtmFloodInit subroutine and then do the following...