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)
@@ -406,6 +407,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa
 
 
  else ! ELSE for LAKE_OPT....
+ call hydro_stop("Invalid lake option supplied to LEVELPOOL_PHYSICS()")
  endif ! ENDIF for LAKE_OPT....
 
  return

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