Merge remote-tracking branch 'liho745fork/liho745/river/bug-fix-4f7ae…

…11'(PR #6568)

MOSART can be forced to crash if negative channel storage values are produced mainly due to round errors,
when the channel storage is completely depleted in a single time step.
The crashing was reported in several cases of E3SM fully coupled runs,
where the kinematic waving routing method was invoked by default.
Jon Wolfe was able to reproduce the crashing error in
a low-res run and print out detailed information to identify the cause.
In this fix, a numerical treatment has been implemented
in the kinematic wave routing method, subroutine Routing_KW(),
to ensure that any channel storage does not get completely
depleted in a single time step.
This treatment is physically reasonable as long as MOSART runs
at a sub-daily time step.
Jon Wolfe has tested this treatment in his low-res coupled run,
and it helped avoid crashing.

[non-BFB] for MOSART

components/elm/cime_config/testdefs/testmods_dirs/elm/lnd_rof_2way/shell_commands

@@ -0,0 +1,2 @@
+./xmlchange LND_NCPL=48
+./xmlchange ROF_NCPL=48

components/elm/cime_config/testdefs/testmods_dirs/elm/lnd_rof_2way/user_nl_mosart

@@ -1 +1,2 @@
-inundflag = .true.
+inundflag = .true.
+delt_mosart = 1800

components/mosart/src/riverroute/MOSART_physics_mod.F90

@@ -680,10 +680,12 @@ subroutine Euler
     end do  ! DLevelH2R
     ! subcycling within MOSART ends
 
-! check for negative channel storage
-    if (negchan < -1.e-10) then
-       write(iulog,*) 'Warning: Negative channel storage found! ',negchan
-!       call shr_sys_abort('mosart: negative channel storage')
+   ! check for negative channel storage
+    if (negchan < -1.e-10 .and. negchan >= -1.e-8) then
+       write(iulog,*) 'Warning: Small negative channel storage found! ',negchan
+    elseif(negchan < -1.e-8) then
+       write(iulog,*) 'Error: Negative channel storage found! ',negchan
+       call shr_sys_abort('mosart: negative channel storage')
     endif
     TRunoff%flow = TRunoff%flow / Tctl%DLevelH2R
     TRunoff%erowm_regi(:,nt_nmud:nt_nsan) = TRunoff%erowm_regi(:,nt_nmud:nt_nsan) / Tctl%DLevelH2R
@@ -876,11 +878,7 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
               TRunoff%erout(iunit,nt) = -TRunoff%vr(iunit,nt) * TRunoff%mr(iunit,nt)
               if(-TRunoff%erout(iunit,nt) > TINYVALUE .and. TRunoff%wr(iunit,nt) + &
                  (TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt) + TRunoff%erout(iunit,nt)) * theDeltaT < TINYVALUE) then
-                 if (sediflag) then
-                  TRunoff%erout(iunit,nt) = -(TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt) + TRunoff%wr(iunit,nt)*MaxStorageDepleted/ theDeltaT)
-                 else
-                  TRunoff%erout(iunit,nt) = -(TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt) + TRunoff%wr(iunit,nt)/ theDeltaT)
-                 end if
+                 TRunoff%erout(iunit,nt) = -(TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt) + TRunoff%wr(iunit,nt)*MaxStorageDepleted/ theDeltaT)
                  if(TRunoff%mr(iunit,nt) > 0._r8) then
                     TRunoff%vr(iunit,nt) = -TRunoff%erout(iunit,nt) / TRunoff%mr(iunit,nt)
                  end if
@@ -918,16 +916,6 @@ subroutine Routing_KW(iunit, nt, theDeltaT)
 
     TRunoff%dwr(iunit,nt) = TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt) + TRunoff%erout(iunit,nt) + temp_gwl
 
-    !if(TRunoff%wr(iunit,nt) < TINYVALUE .and. abs(TRunoff%erout(iunit,nt))> TINYVALUE) then
-    !    write(unit=1111,fmt="(i10, 4(e20.11))") iunit, TRunoff%wr(iunit,nt), TRunoff%erout(iunit,nt), TRunoff%erlateral(iunit,nt) + TRunoff%erin(iunit,nt), TRunoff%dwr(iunit,nt) 
-    !    write(unit=1112,fmt="(2(i10), 4(e20.11))") iunit, TUnit%mask(iunit), TRunoff%vr(iunit,nt), TUnit%rlen(iunit), TUnit%rwidth(iunit), TUnit%areaTotal2(iunit)/TUnit%rwidth(iunit)/TUnit%rlen(iunit)
-    !end if
-
-!    if(iunit==490 .and. nt==1) then
-!        write(unit=1111,fmt="(3(e20.11), 5(f12.4))") TUnit%areaTotal2(iunit),TUnit%areaTotal(iunit),TUnit%area(iunit),TUnit%rdepth(iunit), TUnit%rwidth(iunit), TUnit%rslp(iunit), TUnit%nr(iunit), TUnit%nt(iunit)
-!        write(unit=1112,fmt="(6(e20.11))") TRunoff%wr(iunit,nt), TRunoff%dwr(iunit,nt), TRunoff%erlateral(iunit,nt), TRunoff%erin(iunit,nt), TRunoff%erout(iunit,nt), temp_gwl
-!    end if
-
 ! 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)