A working FLOOD version

[ekluzek]

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)

[jedwards4b]

@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.