Adding level-pool reservoir simulation capability to v01 (#90)

* add folder to start reservoir work

* set debuglevel/verbosity to maximum for test

* Add reservoir demo script

* attribute module

* Add initial dummy reservoir calculation

* reorganize functional inputs

* created globals and reran df from utilities

* minor changes to flow/resevoir code

* moved waterbodies key values to supernetwork_data

* Sequential Network Execution

Works for Parallel, no waterbodies; serial no waterbodies; serial with waterbodies.

NOT YET WORKING for parallel with waterbodies (which may indicate deeper issue...)

* fix referencing in flowveldepth -- full execution now works

* remove never-used file Reservoir.f90

* move reservoir parameters to nested dict

* black nhd_reach_utilities

* split qlateral

Also moved flowveldepth into a purely functional call

* revert nhd_network_utilities_v02 (it should not have changed...)

* Update nhd_network_utilities_v02.py

* minor updates

* blackened

* blackened

* force single precision in reservoir calculation

* added route and replace function for testing of Mainstems_CONUS input DF

* Merge Jdhondia branch + many fixes to parallelization, etc.

More info:
Added qlateral-time-step subdivision capability
+        "--qlateral-time-step",
+        "--qts_subdivisions",

Detect and handle WSL compile (WSL solution appears to lose fidelity if compiled with optimizations)

fixing logic for network head segments downstream of other networks

add cumulative qlat and cumulative total volume in flowveldepth and modify write functions to handle

gather Qlat from all segments of the waterbody

Outstanding issues:
Head segments is never hit for quc
Tuple (False,) is unexplained

* remove tuple-wrapped False

* Use ts for flow backreference

using flowarr[-1] works within the loop but obtains the
wrong value for the arrays references upstream tributaries
-- the upstreams have completed calculation and will
have a complete time series to choose from.

* Add more usage examples

* Use empty dict for flowveldepth_connect

* Reduce volume of functional return to avoid memory error on process pickle

https://bugs.python.org/issue17560

Reviewing that link suggests that we may be following an anti-pattern in pusing
*pushing so much data into a multi-process. Hence, we are looking at the
v02 method using threads...

* Add TODO about network sorting

* Add RnR test

* fix rnr name

* fix return value to have key

* close pool only if parallel

also blackened file

* make networks global

* use static tuple inside temporary ordering list

* added sort to ordered_networks

* updated code for writing nc file for waterbodies

* added custom input arg, added init_waterbody_states and q_initial, changed json test file columns

* addes q_initial_states

* minor changes and added functions to utilities file

* separated functions to remove cascading effect

* sort waterbodies_df and correct search method

This change cut execution time in half or better

* renamed folder to file

* working custom json input and functions

* renamed inputs and removed for loop

* added waterbody_init and renamed a few things

* formatted with black

* added Waterbody_ID_crosswalks to def get_reservoir_restart_from_wrf_hydro(

* simplify references for short-time-step assumption

also fixed qlaterals.

* make placeholders to read warm/initial states

* update download link to latest version of route-link file

* symlink varPrecision

* blackend

* reorganize arguments

* add default output text and nc folders

* make channel initial states work for RnR

* set default initial state dataframe (with all zeros)

* cleanup and blacken

* cleanup and blacken

* Add Waterbody Initial State Parsing

* Blackend; add variable for xs in waterbody crosswalk

* use initial states in actual calculation

* Set custom file to read Pocono LakeParm

* add custom json example usage

* Add Pocono1 Test case

* finalize water body capability

* blackend

* Move All Inputs to Json

* add empty dict for empty custom input frames

* fix levelpool height assignment

Anyuse prior to this fix is completely wrong; the input variable was being dropped
without any assignment

* Cleanup reservoir test script and add docstring

* blackened

* remove old method

* write qlateral to flowveldepth to allow time-step adjustment

* Add sort networks argument

Also clean up and add more useful comments

* read waterbody parameters from separate utility function

* BUGFIX accumulate qup in initial states from all upstreams

The calculation previously was only re-assigning with "=" instead of "+="

* add timing statements for full program and main loop

* multiply by dt for correct qlat accumulation

* lowercase for uname.release

Co-authored-by: Jacob Hreha <jacob.hreha@noaa.gov>
Co-authored-by: jdhondia <juzer.dhondia@noaa.gov>
Former-commit-id: ee21efdbc669631a267a356f08769a8e8ce1c870

Author: 3 people

src/fortran_routing/mc_pylink_v00/Reservoir_singleTS/module_levelpool.f90

@@ -0,0 +1,237 @@
+! This module derived from the NWM repo found here: 
+! https://github.com/NCAR/wrf_hydro_nwm_public/blob
+! /master/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F
+
+! This module defines and instantiates objects
+! for a level pool type reservoir. The level
+! pool reservoir type inherits input and
+! output types from the reservoir base
+! module and calls instantiation of these into
+! sub-objects. The level pool reservoir type
+! also points to types for level pool properties
+! and state and calls instantiation of these into
+! sub-objects. This module also contains the
+! subroutine to run level pool reservoir that is
+! derived from the reservoir base type interface
+! to run reservoir. Running level pool will
+! then call the LEVELPOOL_PHYSICS subroutine, which
+! processes the given inputs, properties, and
+! state for a particular level pool reservoir and
+! returns the output/outflow.
+
+module module_levelpool
+
+    use precis
+    implicit none
+
+contains
+
+    ! ------------------------------------------------
+    !   SUBROUTINE LEVELPOOL
+    ! ------------------------------------------------
+
+    subroutine levelpool_physics(dt,qi0,qi1,ql,ar,we,maxh,wc,wl,dl,oe,oc,oa,H0,H1,qo1)
+
+        !! ----------------------------  argument variables
+        !! All elevations should be relative to a common base (often belev(k))
+
+        ! integer, intent(IN) :: ln      ! lake number
+        real(prec), intent(IN)    :: dt      ! routing period [s]
+        real(prec), intent(IN)    :: qi0     ! inflow at previous timestep (cms)
+        real(prec), intent(IN)    :: qi1     ! inflow at current timestep (cms)
+        real(prec), intent(IN)    :: ql      ! lateral inflow
+        real(prec), intent(IN)    :: ar      ! area of reservoir (km^2)
+        real(prec), intent(IN)    :: we      ! bottom of weir elevation
+        real(prec), intent(IN)    :: wc      ! weir coeff.
+        real(prec), intent(IN)    :: wl      ! weir length (m)
+        real(prec), intent(IN)    :: dl      ! dam length(m)
+        real(prec), intent(IN)    :: oe      ! orifice elevation
+        real(prec), intent(IN)    :: oc      ! orifice coeff.
+        real(prec), intent(IN)    :: oa      ! orifice area (m^2)
+        real(prec), intent(IN)    :: maxh    ! max depth of reservoir before overtop (m)
+        real(prec), intent(IN)    :: H0      ! water elevation height (m)
+        real(prec), intent(OUT)   :: H1      ! water elevation height (m)
+        real(prec), intent(OUT)   :: qo1     ! outflow at current timestep
+
+        !!DJG Add lake option switch here...move up to namelist in future versions...
+        integer :: LAKE_OPT            ! Lake model option (move to namelist later)
+        real(prec)    :: H                   ! Temporary assign of incoming lake el. (m)
+
+        !! ----------------------------  local variables
+        real(prec) :: sap                    ! local surface area values
+        real(prec) :: discharge              ! storage discharge m^3/s
+        real(prec) :: tmp1, tmp2
+        real(prec) :: dh, dh1, dh2, dh3      ! Depth in weir, and height function for 3 order RK
+        real(prec) :: It, Itdt_3, Itdt_2_3   ! inflow hydrographs
+        real(prec) :: maxWeirDepth           !maximum capacity of weir
+        real(prec) :: exp32, quot23, zero, gravity  ! internal variables
+        !real :: hdiff_vol, qdiff_vol   ! water balance check variables
+        !! ----------------------------  subroutine body: from chow, mad mays. pg. 252
+        !! -- determine from inflow hydrograph
+
+
+        !!DJG Set hardwire for LAKE_OPT...move specification of this to namelist in
+        !future versions...
+        LAKE_OPT = 2
+        exp32 = 3._prec/2._prec
+        quot23 = 0.667_prec
+        zero = 0.0_prec
+        gravity = 9.81_prec
+        !hdiff_vol = zero
+        !qdiff_vol = zero
+
+        !!DJG IF-block for lake model option  1 - outflow=inflow, 2 - Chow et al level
+        !pool, .....
+        if (LAKE_OPT == 1) then     ! If-block for simple pass through scheme....
+
+           qo1 = qi1                 ! Set outflow equal to inflow at current time
+           H = H0                    ! Set new lake water elevation to incoming lake el.
+
+        else if (LAKE_OPT == 2) then   ! If-block for Chow et al level pool scheme
+
+           H = H0
+           It = qi0
+           Itdt_3   = qi0 + ((qi1 + ql - qi0) * 0.33_prec)
+           Itdt_2_3 = qi0 + ((qi1 + ql - qi0) * 0.67_prec)
+           maxWeirDepth =  maxh - we
+
+           !assume vertically walled reservoir
+           !remove this when moving to a variable head area volume
+           sap = ar * 1.0E6_prec
+
+           !-- determine Q(dh) from elevation-discharge relationship
+           !-- and dh1
+           dh = H - we
+           if (dh > maxWeirDepth) then
+              dh = maxWeirDepth
+           endif
+
+           tmp1 = oc * oa * sqrt(2._prec * gravity * ( H - oe )) !orifice at capacity
+           tmp2 = wc * wl * (dh ** (exp32))  !weir flows at capacity
+
+           !determine the discharge based on current height
+           if(H > maxh) then
+             discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(exp32)) !overtop
+           else if (dh > zero) then              !! orifice and weir discharge
+             discharge = tmp1 + tmp2
+           else if ( H > oe ) then     !! only orifice flow
+             discharge = oc * oa * sqrt(2._prec * gravity * ( H - oe ) )
+           else
+             discharge = zero   !in the dead pool
+           endif
+
+           if (sap > zero) then
+              dh1 = ((It - discharge)/sap)*dt
+           else
+              dh1 = zero
+           endif
+
+           !-- determine Q(H + dh1/3) from elevation-discharge relationship
+           !-- dh2
+           dh = (H+dh1/3._prec) - we
+           if (dh > maxWeirDepth) then
+              dh = maxWeirDepth
+           endif
+
+           tmp1 = oc * oa * sqrt(2._prec * gravity * ( (H+dh1/3._prec) - oe ) )
+           tmp2 = wc * wl * (dh ** (exp32))
+
+           !determine the discharge based on current height
+           if(H > maxh) then
+             discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(exp32)) !overtop
+           else if (dh > zero) then              !! orifice and weir discharge
+             discharge = tmp1 + tmp2
+           else if ( (H+dh1/3._prec) > oe ) then     !! only orifice flow,not full
+             discharge = oc * oa * sqrt(2._prec * gravity * ( (H+dh1/3._prec) - oe ) )
+           else
+             discharge = zero
+            endif
+
+
+           if (sap > zero) then
+              dh2 = ((Itdt_3 - discharge)/sap)*dt
+           else
+              dh2 = zero
+           endif
+
+           !-- determine Q(H + 2/3 dh2) from elevation-discharge relationship
+           !-- dh3
+           dh = (H + (quot23*dh2)) - we
+           if (dh > maxWeirDepth) then
+              dh = maxWeirDepth
+           endif
+
+           tmp1 = oc * oa * sqrt(2._prec * gravity * ( (H+dh2*quot23) - oe ) )
+           tmp2 = wc * wl * (dh ** (exp32))
+
+           !determine the discharge based on current height
+           if(H > maxh) then  ! overtop condition, not good!
+              discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(exp32)) !overtop
+           else if (dh > zero) then              !! orifice and weir discharge
+              discharge = tmp1 + tmp2
+           else if ( (H+dh2*quot23) > oe ) then     !! only orifice flow,not full
+              discharge = oc * oa * sqrt(2._prec * gravity * ( (H+dh2*quot23) - oe ) )
+           else
+              discharge = zero
+           endif
+
+           if (sap > zero) then
+              dh3 = ((Itdt_2_3 - discharge)/sap)*dt
+           else
+              dh3 = zero
+           endif
+
+           !-- determine dh and H
+           dh = (dh1/4._prec) + (0.75_prec*dh3)
+           H = H + dh
+
+           !-- compute final discharge
+           dh = H - we
+           if (dh > maxWeirDepth) then
+              dh = maxWeirDepth
+           endif
+
+           tmp1 = oc * oa * sqrt(2._prec * gravity * ( H - oe ) )
+           tmp2 = wc * wl * (dh ** (exp32))
+
+           !determine the discharge based on current height
+           if(H > maxh) then  ! overtop condition, not good!
+              discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(exp32)) !overtop
+           else if (dh > zero) then              !! orifice and overtop discharge
+              discharge = tmp1 + tmp2
+           else if ( H > oe ) then     !! only orifice flow,not full
+              discharge = oc * oa * sqrt(2._prec * gravity * ( H - oe ) )
+           else
+              discharge = zero
+           endif
+
+           H1 = H
+           qo1  = discharge  ! return the flow rate from reservoir
+
+        !#ifdef HYDRO_D
+        !#ifndef NCEP_WCOSS
+        !   ! Water balance check
+        !   qdiff_vol = (qi1+ql-qo1)*dt !m3
+        !   hdiff_vol = (H-H0)*sap    !m3
+        !22 format(f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f6.0,2x,f20.1,2x,f20.1)
+        !   open (unit=67, &
+        !     file='lake_massbalance_out.txt', status='unknown',position='append')
+        !   write(67,22) H0, H, qi1, ql, qo1, dt, qdiff_vol, hdiff_vol
+        !   close(67)
+        !#endif
+        !#endif
+
+        !23 format('botof H dh orf wr Q',f8.4,2x,f8.4,2x,f8.3,2x,f8.3,2x,f8.2)
+        !24 format('ofonl H dh sap Q ',f8.4,2x,f8.4,2x,f8.0,2x,f8.2)
+
+
+        else   ! ELSE for LAKE_OPT....
+        endif  ! ENDIF for LAKE_OPT....
+
+        return
+
+    ! ----------------------------------------------------------------
+    end subroutine levelpool_physics
+    ! ----------------------------------------------------------------
+
+end module module_levelpool

src/fortran_routing/mc_pylink_v00/Reservoir_singleTS/reservoirs_nwm.py

@@ -0,0 +1,110 @@
+import traceback
+
+debuglevel = 0
+COMPILE = True
+if COMPILE:
+    try:
+        import subprocess
+
+        fortran_compile_call = []
+        fortran_compile_call.append(r"f2py3")
+        fortran_compile_call.append(r"-c")
+        fortran_compile_call.append(r"varPrecision.f90")
+        fortran_compile_call.append(r"module_levelpool.f90")
+        fortran_compile_call.append(r"-m")
+        fortran_compile_call.append(r"levelpools")
+        # fortran_compile_call.append(r"--opt='-fdefault-real-8'")
+
+        if debuglevel <= -1:
+            print(fortran_compile_call)
+        if debuglevel <= -2:
+            subprocess.run(fortran_compile_call)
+        else:
+            subprocess.run(
+                fortran_compile_call,
+                stdout=subprocess.DEVNULL,
+                stderr=subprocess.DEVNULL,
+            )
+        from levelpools import *
+
+        # import mc_sseg_stime as muskingcunge_module
+    except Exception as e:
+        print(e)
+        if debuglevel <= -1:
+            traceback.print_exc()
+else:
+    from levelpools import *
+
+    # import mc_sseg_stime as muskingcunge_module
+
+
+def reservoirs_calc(
+    dt=None,
+    qi0=None,
+    qi1=None,
+    ql=None,
+    ar=None,
+    we=None,
+    maxh=None,
+    wc=None,
+    wl=None,
+    dl=None,
+    oe=None,
+    oc=None,
+    oa=None,
+    h0=None,
+):
+    """This test function demonstrates a simple connection between
+       the module_levelpool routing function and python via f2py
+       inputs are described by the fortran header:
+        real(prec), intent(IN)    :: qi0     ! inflow at previous timestep (cms)
+        real(prec), intent(IN)    :: qi1     ! inflow at current timestep (cms)
+        real(prec), intent(IN)    :: ql      ! lateral inflow
+        real(prec), intent(IN)    :: ar      ! area of reservoir (km^2)
+        real(prec), intent(IN)    :: we      ! bottom of weir elevation
+        real(prec), intent(IN)    :: wc      ! weir coeff.
+        real(prec), intent(IN)    :: wl      ! weir length (m)
+        real(prec), intent(IN)    :: dl      ! dam length(m)
+        real(prec), intent(IN)    :: oe      ! orifice elevation
+        real(prec), intent(IN)    :: oc      ! orifice coeff.
+        real(prec), intent(IN)    :: oa      ! orifice area (m^2)
+        real(prec), intent(IN)    :: maxh    ! max depth of reservoir before overtop (m)
+        real(prec), intent(IN)    :: H0      ! water elevation height (m)
+
+
+       outputs are, in order the new water depth and new outflow. 
+        real(prec), intent(OUT)   :: H1      ! water elevation height (m)
+        real(prec), intent(OUT)   :: qo1     ! outflow at current timestep
+       """
+
+    # call Fortran routine
+    return module_levelpool.levelpool_physics(
+        dt, qi0, qi1, ql, ar, we, maxh, wc, wl, dl, oe, oc, oa, h0,
+    )
+    # return hc, qdc
+
+
+def main():
+
+    # dummy execution of reservoirs_calc()
+
+    dt = 300
+    qi0 = 500
+    qi1 = 600
+    ql = 400
+    ar = 22
+    we = 20
+    maxh = 30
+    wc = 1.6
+    wl = 20
+    dl = 200
+    oe = 1
+    oc = 0.6
+    oa = 0.2
+    h0 = 22
+
+    print(reservoirs_calc(dt, qi0, qi1, ql, ar, we, maxh, wc, wl, dl, oe, oc, oa, h0,))
+
+
+if __name__ == "__main__":
+    main()

src/fortran_routing/mc_pylink_v00/Reservoir_singleTS/varPrecision.f90

@@ -0,0 +1 @@
+../MC_singleSeg_singleTS/varPrecision.f90