Add lake_option to hydro_namelist, reservoirs to their own new namelist

src/Routing/Reservoirs/Level_Pool/module_levelpool.F

@@ -56,7 +56,8 @@ module module_levelpool
     subroutine levelpool_init(this, water_elevation,  &
         lake_area, weir_elevation, weir_coeffecient, &
         weir_length, dam_length, orifice_elevation, orifice_coefficient, &
-        orifice_area, max_depth, lake_number)
+        orifice_area, max_depth, lake_number, lake_opt)
+
         implicit none
         class(levelpool), intent(inout) :: this ! object being initialized
         real, intent(inout) :: water_elevation           ! meters AMSL
@@ -69,7 +70,9 @@ subroutine levelpool_init(this, water_elevation,  &
         real, intent(in)    :: orifice_coefficient       ! orifice coefficient
         real, intent(in)    :: orifice_area              ! orifice area (meters^2)
         real, intent(in)    :: max_depth                 ! max depth of reservoir before overtop (meters)
-        integer(kind=int64), intent(in) :: lake_number               ! lake number
+        integer(kind=int64), intent(in) :: lake_number   ! lake number
+        integer, intent(in) :: lake_opt                  ! bypass lake physics (2 to use pass-through)
+
         character(len=15)   :: lake_number_string
 
 #ifdef RESERVOIR_D
@@ -114,7 +117,7 @@ subroutine levelpool_init(this, water_elevation,  &
                 call this%properties%init( lake_area,  &
                     weir_elevation, weir_coeffecient, weir_length, dam_length, &
                     orifice_elevation, orifice_coefficient, &
-                    orifice_area, max_depth, lake_number )
+                    orifice_area, max_depth, lake_number, lake_opt )
             end if
             this%pointer_allocation_guard = .true.
 
@@ -169,6 +172,7 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, &
         this%state%water_elevation = water_elevation
 
         call LEVELPOOL_PHYSICS(this%properties%lake_number,                  &
+                               this%properties%lake_opt,                   &
                                previous_timestep_inflow,                     &
                                this%input%inflow,                            &
                                this%output%outflow,                          &
@@ -217,7 +221,7 @@ end subroutine run_levelpool_reservoir
     !   SUBROUTINE LEVELPOOL
     ! ------------------------------------------------
 
-    subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa)
+    subroutine LEVELPOOL_PHYSICS(ln,lake_opt,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa)
 
         !! ----------------------------  argument variables
         !! All elevations should be relative to a common base (often belev(k))
@@ -238,9 +242,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa
         real, intent(IN)    :: oa      ! orifice area (m^2)
         real, intent(IN)    :: maxh    ! max depth of reservoir before overtop (m)
         integer(kind=int64), intent(IN) :: ln      ! lake number
+        integer, intent(in) :: lake_opt ! reservoir physics options (1: levelpool, 2: passthrough)
 
-        !!DJG Add lake option switch here...move up to namelist in future versions...
-        integer :: LAKE_OPT            ! Lake model option (move to namelist later)
         real    :: Htmp                ! Temporary assign of incoming lake el. (m)
 
         !! ----------------------------  local variables
@@ -254,22 +257,20 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa
         !! ----------------------------  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
         Htmp = H   !temporary set of incoming lake water elevation...
         !hdiff_vol = 0.0
         !qdiff_vol = 0.0
 
         !!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....
-
+        if (LAKE_OPT == 2) then     ! If-block for simple pass through scheme....
+#ifdef RESERVOIR_D
+           write(6,*) "LEVELPOOL LAKE_OPT=2, using reservoir passthrough"
+#endif
            qo1 = qi1                 ! Set outflow equal to inflow at current time
            H = Htmp                  ! Set new lake water elevation to incoming lake el.
 
-        else if (LAKE_OPT == 2) then   ! If-block for Chow et al level pool scheme
+        else if (LAKE_OPT == 1) then   ! If-block for Chow et al level pool scheme
 
            It = qi0
            Itdt_3   = qi0 + ((qi1 + ql - qi0) * 0.33)

src/Routing/Reservoirs/Level_Pool/module_levelpool_properties.F

@@ -23,6 +23,7 @@ module module_levelpool_properties
         real :: orifice_area             ! orifice area (meters^2)
         real :: max_depth                ! max depth of reservoir before overtop (meters)
         integer(kind=int64) :: lake_number           ! lake number
+        integer :: lake_opt              ! reservoir physics options (1: levelpool, 2: passthrough)
 
     contains
 
@@ -36,7 +37,7 @@ module module_levelpool_properties
     !Level Pool Properties Constructor
     subroutine levelpool_properties_init(this, lake_area, &
         weir_elevation, weir_coeffecient, weir_length, dam_length, orifice_elevation, &
-        orifice_coefficient, orifice_area, max_depth, lake_number)
+        orifice_coefficient, orifice_area, max_depth, lake_number, lake_opt)
         implicit none
         class(levelpool_properties_interface), intent(inout) :: this ! the type object being initialized
         real, intent(in)    :: lake_area      	        ! area of lake (km^2)
@@ -49,6 +50,7 @@ subroutine levelpool_properties_init(this, lake_area, &
         real, intent(in)    :: orifice_area             ! orifice area (meters^2)
         real, intent(in)    :: max_depth                ! max depth of reservoir before overtop (meters)
         integer(kind=int64), intent(in) :: lake_number              ! lake number
+        integer             :: lake_opt                 ! reservoir physics options (1: levelpool, 2: passthrough)
 
         ! Assign the values passed in to a particular level pool reservoir
         ! properties object's variables.
@@ -62,6 +64,7 @@ subroutine levelpool_properties_init(this, lake_area, &
         this%max_depth = max_depth
         this%lake_number = lake_number
         this%dam_length = dam_length
+        this%lake_opt = lake_opt
 
     end subroutine levelpool_properties_init
 

src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F

@@ -201,7 +201,7 @@ subroutine hybrid_init(this, water_elevation,  &
             ! Initialize level pool reservoir
             call this%state%levelpool_ptr%init(water_elevation, lake_area, &
             weir_elevation, weir_coeffecient, weir_length, dam_length, orifice_elevation, &
-            orifice_coefficient, orifice_area, lake_max_water_elevation, lake_number)
+            orifice_coefficient, orifice_area, lake_max_water_elevation, lake_number, 1)
 
         end if
     end subroutine hybrid_init

src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F

@@ -156,7 +156,7 @@ subroutine rfc_forecasts_init(this, water_elevation,  &
         ! Initialize level pool reservoir
         call this%state%levelpool_ptr%init(water_elevation, lake_area, &
         weir_elevation, weir_coeffecient, weir_length, dam_length, orifice_elevation, &
-        orifice_coefficient, orifice_area, lake_max_water_elevation, lake_number)
+        orifice_coefficient, orifice_area, lake_max_water_elevation, lake_number, 1)
 
         ! Call to initialize time series data object
         call time_series_data%init(start_date, time_series_path, forecast_lookback_hours, &

src/Routing/module_HYDRO_io.F

@@ -10358,8 +10358,8 @@ subroutine read_NSIMLAKES(NLAKES,route_lake_f)
           endif
        else
 !yw for IOC reach based routing, if netcdf lake file is not set from the hydro.namelist,
-!    we will assume that no lake will be assimulated.
-          write(6,*) "No lake nectdf file defined. NLAKES is set to be zero."
+!    we will assume that no lake will be assimilated.
+          write(6,*) "Lakes have been disabled -- NLAKES will be set to zero."
           NLAKES = 0
       endif
 #ifdef MPP_LAND

src/Routing/module_RT.F

@@ -853,7 +853,8 @@ subroutine LandRT_ini(did)
                      rt_domain(did)%ORIFICEC(lake_index),         &
                      rt_domain(did)%ORIFICEA(lake_index),         &
                      rt_domain(did)%LAKEMAXH(lake_index),         &
-                     rt_domain(did)%LAKEIDM(lake_index)            )
+                     rt_domain(did)%LAKEIDM(lake_index),          &
+                     nlst(did)%lake_option)
 
           type is (persistence_levelpool_hybrid)
               call reservoir%init(                                     &

src/hydro_namelists.json

@@ -22,6 +22,7 @@
             "igrid": 1,
             "io_config_outputs": 0,
             "io_form_outputs": 4,
+            "lake_option": 1,
             "land_spatial_meta_flnm": "NULL_specified_in_domain.json",
             "lsmout_domain": 1,
             "nsoil": 4,
@@ -31,7 +32,6 @@
             "output_channelbucket_influx": 0,
             "output_gw": 1,
             "ovrtswcrt": 1,
-            "reservoir_persistence_usgs": false,
             "restart_file": "NULL_specified_in_domain.json",
             "route_lake_f": "NULL_specified_in_domain.json",
             "route_link_f": "NULL_specified_in_domain.json",
@@ -50,6 +50,9 @@
             "udmp_opt": 1,
             "zsoil8": [-0.1, -0.4, -1.0, -2.0]
         },
+        "reservoir_nlist": {
+            "reservoir_persistence_usgs": false
+        },
         "nudging_nlist": {
             "biaswindowbeforet0": false,
             "invdisttimeweightbias": true,
@@ -109,7 +112,9 @@
         "hydro_nlist": {
             "lsmout_domain": 0,
             "output_gw": 0,
-            "output_channelbucket_influx": 2,
+            "output_channelbucket_influx": 2
+        },
+        "reservoir_nlist": {
             "reservoir_persistence_usgs": false
         },
         "nudging_nlist": {}
@@ -186,7 +191,8 @@
             "channel_option": 3,
             "dtrt_ch": 10,
             "rstrt_swc": 1,
-            "udmp_opt": 0
+            "udmp_opt": 0,
+            "lake_option": 0
         },
         "nudging_nlist": {
             "biaswindowbeforet0": "",
@@ -211,6 +217,7 @@
             "gwbasmskfil": "NULL_specified_in_domain.json",
             "dtrt_ch": 10,
             "outlake": 0,
+            "lake_option": 0,
             "rstrt_swc": 0,
             "t0outputflag": 0,
             "udmp_opt": 0

src/template/HYDRO/hydro.namelist

@@ -154,10 +154,40 @@ compound_channel = .FALSE.
 ! Switch to activate channel-loss option (0=no, 1=yes) [Requires Kchan in RouteLink]
 ! channel_loss_option = 0
 
+! Lake / Reservoir options (0=lakes off, 1=level pool (typical default),
+!                           2=passthrough, 3=reservoir DA [see &reservoir_nlist below])
+lake_option = 1
+
 ! Specify the lake parameter file (e.g.: "LAKEPARM.nc").
 ! Note REQUIRED if lakes are on.
 route_lake_f = "./DOMAIN/LAKEPARM.nc"
 
+! Switch to activate baseflow bucket model...(0=none, 1=exp. bucket, 2=pass-through,
+! 4=exp. bucket with area normalized parameters)
+! Option 4 is currently only supported if using reach-based routing with UDMP=1.
+GWBASESWCRT = 1
+
+! Switch to activate bucket model loss (0=no, 1=yes)
+! This option is currently only supported if using reach-based routing with UDMP=1.
+bucket_loss = 0
+
+! Groundwater/baseflow 2d mask specified on land surface model grid (e.g.: "GWBASINS.nc")
+! Note: Only required if baseflow  model is active (1 or 2) and UDMP_OPT=0.
+gwbasmskfil = "./DOMAIN/GWBASINS.nc"
+
+! Groundwater bucket parameter file (e.g.: "GWBUCKPARM.nc")
+GWBUCKPARM_file = "./DOMAIN/GWBUCKPARM.nc"
+
+! User defined mapping, such as NHDPlus: 0=no (default), 1=yes
+UDMP_OPT = 0
+
+! If on, specify the user-defined mapping file (e.g.: "spatialweights.nc")
+!udmap_file = "./DOMAIN/spatialweights.nc"
+
+/
+
+&reservoir_nlist
+
 ! Specify the reservoir parameter file
 reservoir_parameter_file = "./DOMAIN/persistence_parm.nc"
 
@@ -190,28 +220,6 @@ reservoir_rfc_forecasts_time_series_path = "./rfc_timeseries/"
 ! Specify lookback hours to read reservoir RFC forecasts
 reservoir_rfc_forecasts_lookback_hours = 28
 
-! Switch to activate baseflow bucket model...(0=none, 1=exp. bucket, 2=pass-through,
-! 4=exp. bucket with area normalized parameters)
-! Option 4 is currently only supported if using reach-based routing with UDMP=1.
-GWBASESWCRT = 1
-
-! Switch to activate bucket model loss (0=no, 1=yes)
-! This option is currently only supported if using reach-based routing with UDMP=1.
-bucket_loss = 0
-
-! Groundwater/baseflow 2d mask specified on land surface model grid (e.g.: "GWBASINS.nc")
-! Note: Only required if baseflow  model is active (1 or 2) and UDMP_OPT=0.
-gwbasmskfil = "./DOMAIN/GWBASINS.nc"
-
-! Groundwater bucket parameter file (e.g.: "GWBUCKPARM.nc")
-GWBUCKPARM_file = "./DOMAIN/GWBUCKPARM.nc"
-
-! User defined mapping, such as NHDPlus: 0=no (default), 1=yes
-UDMP_OPT = 0
-
-! If on, specify the user-defined mapping file (e.g.: "spatialweights.nc")
-!udmap_file = "./DOMAIN/spatialweights.nc"
-
 /
 
 &NUDGING_nlist