A working FLOOD version

Bắt đầu trong src/riverroute/RtmMod.F90, đọc Rtmini và Rtmrun, bao gồm đường dẫn RtmFloodInit hiện có và các phép tính lũ được đề xuất. Theo dõi cách TRunoff, TUnit và rtmCTL tương tác trước khi triển khai hành vi FLOOD; hoàn thành có nghĩa là có một phiên bản FLOOD được cài đặt và hoạt động, cùng với các bài kiểm thử bao phủ phiên bản đó.

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)