From c6b507085968933bb8e4ecfb8359fb510879351a Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Tue, 7 Jul 2020 13:49:03 -0600 Subject: [PATCH 01/12] new shape options were added to levelpool module. LAKEPARM.nc and hydro.namelist were modified accordingly. --- trunk/NDHMS/Data_Rec/module_namelist.F | 14 +- trunk/NDHMS/Data_Rec/namelist.inc | 3 +- trunk/NDHMS/Data_Rec/rt_include.inc | 7 +- trunk/NDHMS/OrchestratorLayer/config.f90 | 16 +- .../Reservoirs/Level_Pool/module_levelpool.F | 163 +++++++++++++++--- .../Level_Pool/module_levelpool_properties.F | 16 +- .../module_persistence_levelpool_hybrid.F | 9 +- .../RFC_Forecasts/module_rfc_forecasts.F | 8 +- trunk/NDHMS/Routing/module_HYDRO_io.F | 24 ++- trunk/NDHMS/Routing/module_RT.F | 17 ++ 10 files changed, 230 insertions(+), 47 deletions(-) diff --git a/trunk/NDHMS/Data_Rec/module_namelist.F b/trunk/NDHMS/Data_Rec/module_namelist.F index 758db5c85..9b7cb3167 100644 --- a/trunk/NDHMS/Data_Rec/module_namelist.F +++ b/trunk/NDHMS/Data_Rec/module_namelist.F @@ -46,7 +46,8 @@ subroutine read_rt_nlst(nlst) GWBASESWCRT, GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, & sys_cpl, rst_typ, rst_bi_in, rst_bi_out, & gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, gwsoilcpl, & - UDMP_OPT, io_form_outputs, bucket_loss + UDMP_OPT, io_form_outputs, bucket_loss, & + lake_shape_option real:: DTRT_TER,DTRT_CH,dxrt, gwChanCondConstIn, gwChanCondConstOut, gwIhShift character(len=256) :: route_topo_f="" character(len=256) :: route_chan_f="" @@ -138,8 +139,8 @@ subroutine read_rt_nlst(nlst) CHRTOUT_DOMAIN,CHANOBS_DOMAIN,CHRTOUT_GRID,LSMOUT_DOMAIN,& RTOUT_DOMAIN, output_gw, outlake, & frxst_pts_out, udmap_file, UDMP_OPT, GWBUCKPARM_file, bucket_loss, & - io_config_outputs, io_form_outputs, hydrotbl_f, t0OutputFlag, output_channelBucket_influx - + io_config_outputs, io_form_outputs, hydrotbl_f, t0OutputFlag, output_channelBucket_influx, & + lake_shape_option #ifdef WRF_HYDRO_NUDGING namelist /NUDGING_nlist/ nudgingParamFile, netwkReExFile, & readTimesliceParallel, temporalPersistence, & @@ -359,6 +360,7 @@ subroutine read_rt_nlst(nlst) call mpp_land_bcast_int1(RTOUT_DOMAIN) call mpp_land_bcast_int1(UDMP_OPT) call mpp_land_bcast_int1(reservoir_data_ingest) + !call mpp_land_bcast_int1(lake_shape_option) #ifdef WRF_HYDRO_NUDGING call mpp_land_bcast_char(256, nudgingParamFile ) call mpp_land_bcast_char(256, netwkReExFile ) @@ -401,6 +403,7 @@ subroutine read_rt_nlst(nlst) nlst%DTRT_TER = DTRT_TER nlst%DTRT_CH = DTRT_CH nlst%DTCT = DTRT_CH ! small time step for grid based channel routing + nlst%lake_shape_option = lake_shape_option #ifdef MPP_LAND if(my_id .eq. IO_id) then @@ -519,6 +522,7 @@ subroutine read_rt_nlst(nlst) write(6,*) " nlst%compound_channel ", compound_channel write(6,*) " nlst%route_lake_f ",route_lake_f write(6,*) " nlst%reservoir_parameter_file ", reservoir_parameter_file + write(6,*) " nlst%lake_shape_option ", lake_shape_option write(6,*) " nlst%reservoir_usgs_timeslice_path ", reservoir_usgs_timeslice_path write(6,*) " nlst%reservoir_usace_timeslice_path ", reservoir_usace_timeslice_path write(6,*) " nlst%reservoir_observation_lookback_hours ", nlst%reservoir_observation_lookback_hours @@ -803,7 +807,9 @@ subroutine rt_nlst_check(nlst) if (.not. fileExists) call hydro_stop('hydro.namelist ERROR: reservoir_parameter_file not found.') endif end if - + if( (nlst%lake_shape_option .lt. 0 ) .or. (nlst%lake_shape_option .gt. 5) ) then + call hydro_stop('hydro.namelist ERROR: Invalid lake_shape_option specified') + endif if(nlst%reservoir_persistence_usgs) then if(len(trim(nlst%reservoir_usgs_timeslice_path)) .eq. 0) then call hydro_stop('hydro.namelist ERROR: You MUST specify a reservoir_usgs timeslice_path for reservoir USGS persistence capability.') diff --git a/trunk/NDHMS/Data_Rec/namelist.inc b/trunk/NDHMS/Data_Rec/namelist.inc index 9598ef7bf..7dd440363 100644 --- a/trunk/NDHMS/Data_Rec/namelist.inc +++ b/trunk/NDHMS/Data_Rec/namelist.inc @@ -34,7 +34,8 @@ SUBRTSWCRT, OVRTSWCRT, AGGFACTRT, & GWBASESWCRT, GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, & sys_cpl, gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, & - gwsoilcpl, UDMP_OPT, bucket_loss + gwsoilcpl, UDMP_OPT, bucket_loss, & + lake_shape_option logical:: GwPreDiag, GwSpinUp real:: DTRT_TER,DTRT_CH, DTCT, dxrt0, gwChanCondConstIn, gwChanCondConstOut, gwIhShift character(len=256) :: route_topo_f="" diff --git a/trunk/NDHMS/Data_Rec/rt_include.inc b/trunk/NDHMS/Data_Rec/rt_include.inc index 420718f6a..d1812534d 100644 --- a/trunk/NDHMS/Data_Rec/rt_include.inc +++ b/trunk/NDHMS/Data_Rec/rt_include.inc @@ -229,9 +229,10 @@ REAL, allocatable, DIMENSION(:) :: LAKEMAXH !maximum depth (m) REAL, allocatable, DIMENSION(:) :: WEIRC !coeff of overtop weir REAL, allocatable, DIMENSION(:) :: WEIRH !depth of Lake coef - - - + REAL, allocatable, DIMENSION(:) :: LAKEMAXDEPTH !max depth (m) based on actual depth + REAL, allocatable, DIMENSION(:) :: LAKEVOL !volume of lake based on bathymetry (km^3) + REAL, allocatable, DIMENSION(:) :: LAKEPOLYA !coefficient a in polynomial function V=a*h^b + REAL, allocatable, DIMENSION(:) :: LAKEPOLYB !coefficient b in polynomial function V=a*h^b !DJG Modified namelist for routing and agg. variables real Z_tmp diff --git a/trunk/NDHMS/OrchestratorLayer/config.f90 b/trunk/NDHMS/OrchestratorLayer/config.f90 index 8e15f6374..7f6d66298 100644 --- a/trunk/NDHMS/OrchestratorLayer/config.f90 +++ b/trunk/NDHMS/OrchestratorLayer/config.f90 @@ -101,7 +101,8 @@ module config_base SUBRTSWCRT,OVRTSWCRT,AGGFACTRT, & GWBASESWCRT, GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, & sys_cpl, gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, & - gwsoilcpl, UDMP_OPT + gwsoilcpl, UDMP_OPT, & + lake_shape_option logical:: GwPreDiag, GwSpinUp real:: DTRT_TER,DTRT_CH, DTCT, dxrt0, gwChanCondConstIn, gwChanCondConstOut, gwIhShift character(len=256) :: route_topo_f="" @@ -413,7 +414,9 @@ subroutine rt_nlst_check(self) if (.not. fileExists) call hydro_stop('hydro.namelist ERROR: reservoir_parameter_file not found.') endif end if - + if( (self%lake_shape_option .lt. 0 ) .or. (self%lake_shape_option .gt. 5) ) then + call hydro_stop('hydro.namelist ERROR: Invalid lake_shape_option specified') + endif if(self%reservoir_persistence_usgs) then if(len(trim(self%reservoir_usgs_timeslice_path)) .eq. 0) then call hydro_stop('hydro.namelist ERROR: You MUST specify a reservoir_usgs_timeslice_path for & @@ -462,7 +465,8 @@ subroutine init_namelist_rt_field(did) GWBASESWCRT, GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, & sys_cpl, rst_typ, rst_bi_in, rst_bi_out, & gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, gwsoilcpl, & - UDMP_OPT, io_form_outputs, bucket_loss + UDMP_OPT, io_form_outputs, bucket_loss, & + lake_shape_option real:: DTRT_TER,DTRT_CH,dxrt, gwChanCondConstIn, gwChanCondConstOut, gwIhShift character(len=256) :: route_topo_f="" character(len=256) :: route_chan_f="" @@ -552,7 +556,8 @@ subroutine init_namelist_rt_field(did) CHRTOUT_DOMAIN,CHANOBS_DOMAIN,CHRTOUT_GRID,LSMOUT_DOMAIN,& RTOUT_DOMAIN, output_gw, outlake, & frxst_pts_out, udmap_file, UDMP_OPT, GWBUCKPARM_file, bucket_loss, & - io_config_outputs, io_form_outputs, hydrotbl_f, t0OutputFlag, output_channelBucket_influx + io_config_outputs, io_form_outputs, hydrotbl_f, t0OutputFlag, output_channelBucket_influx, & + lake_shape_option #ifdef WRF_HYDRO_NUDGING namelist /NUDGING_nlist/ nudgingParamFile, netwkReExFile, & @@ -713,6 +718,7 @@ subroutine init_namelist_rt_field(did) nlst(did)%DTRT_TER = DTRT_TER nlst(did)%DTRT_CH = DTRT_CH nlst(did)%DTCT = DTRT_CH ! small time step for grid based channel routing + nlst(did)%lake_shape_option = lake_shape_option ! Some fields haven't been initialized yet (e.g. DT) @@ -735,6 +741,8 @@ subroutine init_namelist_rt_field(did) call hydro_stop("module_namelist: DT not a multiple of DTRT_CH") endif + print *, "Lake Shape Option is:", nlst(did)%lake_shape_option + nlst(did)%SUBRTSWCRT = SUBRTSWCRT nlst(did)%OVRTSWCRT = OVRTSWCRT nlst(did)%dxrt0 = dxrt diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index d064b29e2..6e10da382 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -22,7 +22,7 @@ module module_levelpool use module_reservoir, only: reservoir, reservoir_input, & reservoir_output use module_hydro_stop, only: HYDRO_stop - + use config_base, only: nlst #ifdef RESERVOIR_D use module_reservoir_utilities, only: create_levelpool_diagnostic_log_file, & log_levelpool_diagnostic_data @@ -55,7 +55,9 @@ 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, & + max_depth_full_lake, lake_vol, lake_polya, lake_polyb, & + lake_number) implicit none class(levelpool), intent(inout) :: this ! object being initialized real, intent(inout) :: water_elevation ! meters AMSL @@ -67,7 +69,11 @@ subroutine levelpool_init(this, water_elevation, & real, intent(in) :: orifice_elevation ! orifice elevation (meters AMSL) 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) + real, intent(in) :: max_depth ! max depth of reservoir before overtop (meters) based on active reservoir depth + real, intent(in) :: max_depth_full_lake ! max depth of reservoir before overtop (meters) based on full lake depth + real, intent(in) :: lake_vol ! volume of lake based on bathymetry (km^3) + real, intent(in) :: lake_polya ! coefficient a in polynomial function V=a*h^b + real, intent(in) :: lake_polyb ! coefficient b in polynomial function V=a*h^b integer, intent(in) :: lake_number ! lake number character(len=15) :: lake_number_string @@ -113,7 +119,9 @@ 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, & + max_depth_full_lake, lake_vol, lake_polya, lake_polyb, & + lake_number ) end if this%pointer_allocation_guard = .true. @@ -182,8 +190,11 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, & this%properties%dam_length, & this%properties%orifice_elevation, & this%properties%orifice_coefficient, & - this%properties%orifice_area & - ) + this%properties%orifice_area, & + this%properties%max_depth_full_lake, & + this%properties%lake_vol, & + this%properties%lake_polya, & + this%properties%lake_polyb) ! Update output variable returned from this subroutine outflow = this%output%outflow @@ -216,7 +227,8 @@ 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,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa, & + maxdepth,vol,polya,polyb) !! ---------------------------- argument variables !! All elevations should be relative to a common base (often belev(k)) @@ -235,15 +247,28 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa real, intent(IN) :: oe ! orifice elevation real, intent(IN) :: oc ! orifice coeff. real, intent(IN) :: oa ! orifice area (m^2) - real, intent(IN) :: maxh ! max depth of reservoir before overtop (m) + real, intent(IN) :: maxh ! max depth of reservoir before overtop (m) based on active reservoir depth + real, intent(IN) :: maxdepth ! max depth of reservoir before overtop (m) based on full lake depth + real, intent(IN) :: vol ! volume of lake based on bathymetry (km^3) + real, intent(IN) :: polya ! coefficient a in polynomial function V=a*h^b + real, intent(IN) :: polyb ! coefficient b in polynomial function V=a*h^b integer, intent(IN) :: ln ! lake number - !!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) + integer, parameter :: did = 1 + integer :: lake_shape_option ! Lake model option + ! lake_shape_option options: + ! -1 ---> Simple pass through scheme + ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) + ! 1 ---> BOX shape with depth based on full lake depth + ! 2 ---> CONIC shape with depth based on full lake depth + ! 3 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 4 ---> ELLIPSOID shape with depth based on full lake depth + ! 5 ---> H-A-V RELATIONSHIP depth based on full lake depth !! ---------------------------- local variables - real :: sap ! local surface area values + real :: Htmp ! Temporary assign of incoming lake el. (m) + real :: sap ! local surface area values in m2 + real :: vtp ! local total volume values in m3 real :: discharge ! storage discharge m^3/s real :: tmp1, tmp2 real :: dh, dh1, dh2, dh3 ! Depth in weir, and height function for 3 order RK @@ -253,31 +278,38 @@ 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 + ! read lake shape function from hydro.namelist: + lake_shape_option = nlst(did)%lake_shape_option + !print *, "lake_shape_option is:", lake_shape_option - !!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-block for lake shape option: + if (lake_shape_option == -1) then ! If-block for simple pass through scheme.... 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-block for Chow et al level pool scheme + + ! Define reservoir max depth and rescale orifice and weir elevation + ! if using lake_shape_option >= 1: + if (lake_shape_option == 0) then ! define max depth based on active reservoir depth (NWM v1.x - v2.1) + ! ADD ADJUSTENTS FOR LAKE_SHAPE_OPT = 0 HERE! + else ! define max depth based on full lake depth + ! ADD ADJUSTENTS FOR LAKE_SHAPE_OPT > 0 HERE! + endif It = qi0 Itdt_3 = qi0 + ((qi1 + ql - qi0) * 0.33) Itdt_2_3 = qi0 + ((qi1 + ql - qi0) * 0.67) maxWeirDepth = maxh - we - !assume vertically walled reservoir - !remove this when moving to a variable head area volume + !convert area and volume to m2 and m3 sap = ar * 1.0E6 + vtp = vol * 1.0E9 !-- determine Q(dh) from elevation-discharge relationship !-- and dh1 @@ -301,7 +333,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0) then - dh1 = ((It - discharge)/sap)*dt + !dh1 = ((It - discharge)/sap)*dt + call STAGE_STORAGE(dh1,maxdepth,H,It,discharge,sap,vtp,dt,lake_shape_option,polya,polya) else dh1 = 0.0 endif @@ -329,7 +362,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0.0) then - dh2 = ((Itdt_3 - discharge)/sap)*dt + !dh2 = ((Itdt_3 - discharge)/sap)*dt + call STAGE_STORAGE(dh2,maxdepth,H,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polya) else dh2 = 0.0 endif @@ -356,7 +390,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0.0) then - dh3 = ((Itdt_2_3 - discharge)/sap)*dt + !dh3 = ((Itdt_2_3 - discharge)/sap)*dt + call STAGE_STORAGE(dh3,maxdepth,H,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polya) else dh3 = 0.0 endif @@ -403,14 +438,86 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa !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.... - + endif ! ENDIF for lake_shape_option.... + !print *, "TEST STAGE_STORAGE..." + !print *, dh, dh1, dh2, dh3 + !print *, dh, dh1_test, dh2_test, dh3_test return ! ---------------------------------------------------------------- end subroutine LEVELPOOL_PHYSICS ! ---------------------------------------------------------------- + ! ------------------------------------------------ + ! SUBROUTINE STAGE-STORAGE + ! ------------------------------------------------ + ! This subroutine calculates change in sotrage water level (Deltah) given the initial water level (h1), + ! reservoir shape option (lake_shape_option), inflow (Inflow_dt), and discharge (res_discharge), and surface area (surf_area) + subroutine STAGE_STORAGE(Deltah,Dmax,h1,Inflow_dt,res_discharge,surf_area,total_vol, & + Deltat,LAKE_SHAPE_OPT, polyhV_a, polyhV_b) + + !! ---------------------------- argument variables + !! All elevations should be relative to a common base (often belev(k)) + + real, intent(OUT) :: Deltah ! change in sotrage water level (m) + real, intent(IN) :: Dmax ! Maximum depth of reservoir (m) + real, intent(IN) :: h1 ! initial water level (m) + real, intent(IN) :: Inflow_dt ! inflow (m^3/s) + real, intent(IN) :: res_discharge ! storage discharge (m^3/s) + real, intent(IN) :: surf_area ! surface area (m2) + real, intent(IN) :: total_vol ! total volume (m3) + real, intent(IN) :: Deltat ! routing period (s) + real, intent(IN) :: polyhV_a ! coefficient a in polynomial function V=a*h^b + real, intent(IN) :: polyhV_b ! coefficient b in polynomial function V=a*h^b + integer, intent(IN) :: LAKE_SHAPE_OPT ! LAKE_SHAPE_OPT (lake_shape_option) options from namelist: + ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) + ! 1 ---> BOX shape with depth based on full lake depth + ! 2 ---> CONIC shape with depth based on full lake depth + ! 3 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 4 ---> ELLIPSOID shape with depth based on full lake depth + ! 5 ---> H-A-V RELATIONSHIP depth based on full lake depth + + !! ---------------------------- local variables + real, parameter :: pi = 3.14159265 + real :: R_cone ! cone base radius (m) + real :: S1 ! initial reservoir storage (m3) + real :: S2 ! reservoir storage after change (m3) + real :: dS ! change in reservoir storage (m3) + real :: h2 ! cone base radius (m) + !! ---------------------------- + + ! change in storage (m3) + dS = (Inflow_dt - res_discharge)*Deltat + + if ( (LAKE_SHAPE_OPT == 0) .OR. (LAKE_SHAPE_OPT == 1) ) then ! BOX shape + h2 = (dS / surf_area) + h1 + + elseif (LAKE_SHAPE_OPT == 2) then ! CONIC shape + R_cone = (surf_area/pi)**0.5 ! cone base radius (m) + S1 = (pi/3.) * (R_cone/Dmax)**2. * (h1**3.) ! storage before change in volume (m3) + S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) + h2 = (3.*S2/pi*(Dmax/R_cone)**2.)**(1./3.) ! elevation after change in storage in (m) + + elseif (LAKE_SHAPE_OPT == 3) then ! RECTANGLUAR PRISM shape + S1 = (surf_area/2./Dmax)*(h1**2.) ! storage before change in volume (m3) + S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) + h2 = (2.*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) + + elseif (LAKE_SHAPE_OPT == 4) then ! ELLIPSOID shape + !S1 = (surf_area/2/Dmax)*(h1**2) ! storage before change in volume (m3) + !S2 = S1 + dS ! storage after change in volume (m3) + !h2 = (2*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) + + elseif (LAKE_SHAPE_OPT == 5) then ! H-A-V RELATIONSHIP + S1 = polyhV_a*(h1**polyhV_b) * 1.0E9 ! storage before change in volume (m3) + S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) + h2 = (S2/1.0E9/polyhV_a)**(1./polyhV_a) ! elevation after change in storage in (m) + endif + + Deltah = h2-h1 ! change in sotrage water level (m) + return + + ! ---------------------------------------------------------------- + end subroutine STAGE_STORAGE + ! ---------------------------------------------------------------- end module module_levelpool diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool_properties.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool_properties.F index a285da9a1..b406456a0 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool_properties.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool_properties.F @@ -21,6 +21,10 @@ module module_levelpool_properties real :: orifice_coefficient ! orifice coefficient real :: orifice_area ! orifice area (meters^2) real :: max_depth ! max depth of reservoir before overtop (meters) + real :: max_depth_full_lake ! max depth of reservoir before overtop (meters) based on full lake depth + real :: lake_vol ! volume of lake based on bathymetry (km^3) + real :: lake_polya ! coefficient a in polynomial function V=a*h^b + real :: lake_polyb ! coefficient b in polynomial function V=a*h^b integer :: lake_number ! lake number contains @@ -35,7 +39,9 @@ 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, & + max_depth_full_lake, lake_vol, lake_polya, lake_polyb, & + lake_number) implicit none class(levelpool_properties_interface), intent(inout) :: this ! the type object being initialized real, intent(in) :: lake_area ! area of lake (km^2) @@ -47,6 +53,10 @@ subroutine levelpool_properties_init(this, lake_area, & 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) + real, intent(in) :: max_depth_full_lake ! max depth of reservoir before overtop (meters) based on full lake depth + real, intent(in) :: lake_vol ! volume of lake based on bathymetry (km^3) + real, intent(in) :: lake_polya ! coefficient a in polynomial function V=a*h^b + real, intent(in) :: lake_polyb ! coefficient b in polynomial function V=a*h^b integer, intent(in) :: lake_number ! lake number ! Assign the values passed in to a particular level pool reservoir @@ -59,6 +69,10 @@ subroutine levelpool_properties_init(this, lake_area, & this%orifice_coefficient = orifice_coefficient this%orifice_area = orifice_area this%max_depth = max_depth + this%max_depth_full_lake = max_depth_full_lake + this%lake_vol = lake_vol + this%lake_polya = lake_polya + this%lake_polyb = lake_polyb this%lake_number = lake_number this%dam_length = dam_length diff --git a/trunk/NDHMS/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F b/trunk/NDHMS/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F index 945e761eb..c408ecf1d 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F +++ b/trunk/NDHMS/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F @@ -86,7 +86,10 @@ subroutine hybrid_init(this, water_elevation, & integer, intent(in) :: observation_update_time_interval_seconds character(len=15) :: lake_number_string character(len=15) :: reservoir_type_string - + real :: max_depth_full_lake ! max depth of reservoir before overtop (meters) based on full lake depth + real :: lake_vol ! volume of lake based on bathymetry (km^3) + real :: lake_polya ! coefficient a in polynomial function V=a*h^b + real :: lake_polyb ! coefficient b in polynomial function V=a*h^b #ifdef RESERVOIR_D ! Create diagnostic log file only for development/debugging purposes @@ -200,7 +203,9 @@ 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, & + max_depth_full_lake, lake_vol, lake_polya, lake_polyb, & + lake_number) end if end subroutine hybrid_init diff --git a/trunk/NDHMS/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F b/trunk/NDHMS/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F index ba03a5605..7151c2193 100644 --- a/trunk/NDHMS/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F +++ b/trunk/NDHMS/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F @@ -80,6 +80,10 @@ subroutine rfc_forecasts_init(this, water_elevation, & integer :: update_offset_seconds character(len=15) :: lake_number_string character(len=15) :: reservoir_type_string + real :: max_depth_full_lake ! max depth of reservoir before overtop (meters) based on full lake depth + real :: lake_vol ! volume of lake based on bathymetry (km^3) + real :: lake_polya ! coefficient a in polynomial function V=a*h^b + real :: lake_polyb ! coefficient b in polynomial function V=a*h^b #ifdef RESERVOIR_D ! Create diagnostic log file only for development/debugging purposes @@ -155,7 +159,9 @@ 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, & + max_depth_full_lake, lake_vol, lake_polya, lake_polyb, & + lake_number) ! Call to initialize time series data object call time_series_data%init(start_date, time_series_path, forecast_lookback_hours, & diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index 4e60407dc..b5f318d92 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -7354,6 +7354,7 @@ subroutine READ_CHROUTING1( & integer, intent(INOUT), dimension(:) :: reservoir_type character(len=*), intent(in) :: reservoir_parameter_file real, intent(INOUT), dimension(:) :: LATLAKE,LONLAKE,ELEVLAKE +real, allocatable, dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB real, intent(INOUT), dimension(:) :: ChSSlp, Bw, Tw real, intent(INOUT), dimension(:) :: Tw_CC, n_CC ! channel properties of compund @@ -7368,6 +7369,7 @@ subroutine READ_CHROUTING1( & integer :: did logical :: fexist + did = 1 !--------------------------------------------------------- @@ -7470,7 +7472,9 @@ subroutine READ_CHROUTING1( & call read_route_lake_netcdf(trim(route_lake_f),HRZAREA, & LAKEMAXH, WEIRH, WEIRC,WEIRL, DAML, ORIFICEC, & ORIFICEA, ORIFICEE, reservoir_type_specified, reservoir_type, & - reservoir_parameter_file, LAKEIDM, latlake, lonlake, ELEVLAKE, NLAKES) + reservoir_parameter_file, LAKEIDM, latlake, lonlake, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & + ELEVLAKE, NLAKES) else open(unit=79,file=trim(route_lake_f), form='formatted',status='old') write(6,*) "Before read LAKEPARM from text ", trim(route_lake_f) @@ -7508,6 +7512,10 @@ subroutine READ_CHROUTING1( & call mpp_land_bcast_real(NLAKES,ORIFICEE) call mpp_land_bcast_real(NLAKES,LATLAKE ) call mpp_land_bcast_real(NLAKES,LONLAKE ) + call mpp_land_bcast_real(NLAKES,LAKEMAXDEPTH ) + call mpp_land_bcast_real(NLAKES,LAKEVOL ) + call mpp_land_bcast_real(NLAKES,LAKEPOLYA ) + call mpp_land_bcast_real(NLAKES,LAKEPOLYB ) call mpp_land_bcast_real(NLAKES,ELEVLAKE) call mpp_land_bcast_int(NLAKES, reservoir_type) endif @@ -8285,6 +8293,7 @@ subroutine readLinkSL( GNLINKSL,NLINKSL,route_link_f, route_lake_f, maxorder, & character(len=*), intent(in) :: reservoir_parameter_file REAL, intent(out), dimension(:) :: HRZAREA,LAKEMAXH, WEIRC, WEIRL, DAML, ORIFICEC, WEIRH, & ORIFICEA, ORIFICEE, ELEVLAKE + real, allocatable, dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB !end NLAKES INTEGER, dimension(GNLINKSL) :: tmpLAKEIDA, tmpLINKID, tmpTO_NODE, tmpTYPEL, tmpORDER @@ -8370,7 +8379,9 @@ subroutine readLinkSL( GNLINKSL,NLINKSL,route_link_f, route_lake_f, maxorder, & call read_route_lake_netcdf(route_lake_f,HRZAREA, & LAKEMAXH, WEIRH, WEIRC, WEIRL, DAML, ORIFICEC, & ORIFICEA, ORIFICEE, reservoir_type_specified, reservoir_type, reservoir_parameter_file, & - LAKEIDM, latlake, lonlake, ELEVLAKE, NLAKES) + LAKEIDM, latlake, lonlake, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & + ELEVLAKE, NLAKES) endif !!- initialize channel if missing in input @@ -8907,7 +8918,9 @@ subroutine read_route_lake_netcdf(route_lake_file, & HRZAREA, LAKEMAXH, WEIRH, WEIRC, WEIRL, DAML, & ORIFICEC, ORIFICEA, ORIFICEE, reservoir_type_specified, & reservoir_type, reservoir_parameter_file, & - LAKEIDM, lakelat, lakelon, ELEVLAKE, NLAKES) + LAKEIDM, lakelat, lakelon, & + LAKEMAXDEPTH, LAKEVOL, LAKE_POLYA, LAKE_POLYB, & + ELEVLAKE, NLAKES) implicit none character(len=*), intent(in) :: route_lake_file @@ -8917,6 +8930,7 @@ subroutine read_route_lake_netcdf(route_lake_file, & integer, dimension(:), intent(out) :: LAKEIDM real, dimension(:), intent(out) :: HRZAREA, LAKEMAXH, WEIRC, WEIRL, WEIRH, DAML real, dimension(:), intent(out) :: ORIFICEC, ORIFICEA, ORIFICEE, lakelat, lakelon + real, dimension(:), intent(out) :: LAKEVOL, LAKEMAXDEPTH, LAKE_POLYA, LAKE_POLYB real, dimension(:), intent(out) :: ELEVLAKE integer, dimension(:), intent(out) :: reservoir_type @@ -8948,6 +8962,10 @@ subroutine read_route_lake_netcdf(route_lake_file, & call get_1d_netcdf_real(ncid, 'OrificeE', ORIFICEE, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'lat', lakelat, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'lon', lakelon, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_a', LAKE_POLYA, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_b', LAKE_POLYB, 'read_route_lake_netcdf', .TRUE.) !remove the alt var. and add initial fractional depth var. LKR/DY call get_1d_netcdf_real(ncid, 'ifd', ELEVLAKE, 'read_route_lake_netcdf', .FALSE.) diff --git a/trunk/NDHMS/Routing/module_RT.F b/trunk/NDHMS/Routing/module_RT.F index e71627b42..5bb765e3b 100644 --- a/trunk/NDHMS/Routing/module_RT.F +++ b/trunk/NDHMS/Routing/module_RT.F @@ -383,6 +383,11 @@ subroutine rt_allocate(did,ix,jx,ixrt,jxrt,nsoil,CHANRTSWCRT) allocate( rt_domain(did)%ORIFICEC(NLAKES) ) allocate( rt_domain(did)%ORIFICEA(NLAKES) ) allocate( rt_domain(did)%ORIFICEE(NLAKES) ) + allocate( rt_domain(did)%LAKEMAXDEPTH(NLAKES) ) + allocate( rt_domain(did)%LAKEVOL(NLAKES) ) + allocate( rt_domain(did)%LAKEPOLYA(NLAKES) ) + allocate( rt_domain(did)%LAKEPOLYB(NLAKES) ) + rt_domain(did)%HRZAREA = 0.0 rt_domain(did)%WEIRH = 0.0 @@ -394,6 +399,10 @@ subroutine rt_allocate(did,ix,jx,ixrt,jxrt,nsoil,CHANRTSWCRT) rt_domain(did)%ORIFICEC = 0.0 rt_domain(did)%ORIFICEA = 0.0 rt_domain(did)%ORIFICEE = 0.0 + rt_domain(did)%LAKEMAXDEPTH = 0.0 + rt_domain(did)%LAKEVOL = 0.0 + rt_domain(did)%LAKEPOLYA = 0.0 + rt_domain(did)%LAKEPOLYB = 0.0 rt_domain(did)%reservoir_type = 1 rt_domain(did)%final_reservoir_type = 1 rt_domain(did)%reservoir_assimilated_value = -9999.0 @@ -834,6 +843,10 @@ 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)%LAKEMAXDEPTH(lake_index), & + rt_domain(did)%LAKEVOL(lake_index), & + rt_domain(did)%LAKEPOLYA(lake_index), & + rt_domain(did)%LAKEPOLYB(lake_index), & rt_domain(did)%LAKEIDM(lake_index) ) type is (persistence_levelpool_hybrid) @@ -1656,6 +1669,10 @@ subroutine LandRT_ini(did) if(allocated(rt_domain(did)%ORIFICEC)) deallocate(rt_domain(did)%ORIFICEC) if(allocated(rt_domain(did)%ORIFICEA)) deallocate(rt_domain(did)%ORIFICEA) if(allocated(rt_domain(did)%LAKEMAXH)) deallocate(rt_domain(did)%LAKEMAXH) + if(allocated(rt_domain(did)%LAKEMAXDEPTH)) deallocate(rt_domain(did)%LAKEMAXDEPTH) + if(allocated(rt_domain(did)%LAKEVOL)) deallocate(rt_domain(did)%LAKEVOL) + if(allocated(rt_domain(did)%LAKEPOLYA)) deallocate(rt_domain(did)%LAKEPOLYA) + if(allocated(rt_domain(did)%LAKEPOLYB)) deallocate(rt_domain(did)%LAKEPOLYB) end subroutine LandRT_ini From 5a708d2d1665738cacb5d0ea04974cfef63c3bee Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Fri, 24 Jul 2020 10:18:05 -0600 Subject: [PATCH 02/12] adjustments made to levelpool physics and also reading the new variables from LAKEPARM file --- trunk/NDHMS/Data_Rec/module_namelist.F | 2 +- trunk/NDHMS/OrchestratorLayer/config.f90 | 2 +- .../Reservoirs/Level_Pool/module_levelpool.F | 85 +++--- trunk/NDHMS/Routing/module_HYDRO_io.F | 28 +- trunk/NDHMS/Routing/module_RT.F | 2 + trunk/NDHMS/Run/hydro.namelist | 265 ++++++++++++++++++ 6 files changed, 328 insertions(+), 56 deletions(-) create mode 100644 trunk/NDHMS/Run/hydro.namelist diff --git a/trunk/NDHMS/Data_Rec/module_namelist.F b/trunk/NDHMS/Data_Rec/module_namelist.F index 9b7cb3167..c97daefbd 100644 --- a/trunk/NDHMS/Data_Rec/module_namelist.F +++ b/trunk/NDHMS/Data_Rec/module_namelist.F @@ -807,7 +807,7 @@ subroutine rt_nlst_check(nlst) if (.not. fileExists) call hydro_stop('hydro.namelist ERROR: reservoir_parameter_file not found.') endif end if - if( (nlst%lake_shape_option .lt. 0 ) .or. (nlst%lake_shape_option .gt. 5) ) then + if( (nlst%lake_shape_option .lt. 0 ) .or. (nlst%lake_shape_option .gt. 3) ) then call hydro_stop('hydro.namelist ERROR: Invalid lake_shape_option specified') endif if(nlst%reservoir_persistence_usgs) then diff --git a/trunk/NDHMS/OrchestratorLayer/config.f90 b/trunk/NDHMS/OrchestratorLayer/config.f90 index 7f6d66298..324cd2485 100644 --- a/trunk/NDHMS/OrchestratorLayer/config.f90 +++ b/trunk/NDHMS/OrchestratorLayer/config.f90 @@ -414,7 +414,7 @@ subroutine rt_nlst_check(self) if (.not. fileExists) call hydro_stop('hydro.namelist ERROR: reservoir_parameter_file not found.') endif end if - if( (self%lake_shape_option .lt. 0 ) .or. (self%lake_shape_option .gt. 5) ) then + if( (self%lake_shape_option .lt. 0 ) .or. (self%lake_shape_option .gt. 3) ) then call hydro_stop('hydro.namelist ERROR: Invalid lake_shape_option specified') endif if(self%reservoir_persistence_usgs) then diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index 6e10da382..52e96d2ba 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -259,11 +259,9 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa ! lake_shape_option options: ! -1 ---> Simple pass through scheme ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) - ! 1 ---> BOX shape with depth based on full lake depth - ! 2 ---> CONIC shape with depth based on full lake depth - ! 3 ---> RECTANGLUAR PRISM shape with depth based on full lake depth - ! 4 ---> ELLIPSOID shape with depth based on full lake depth - ! 5 ---> H-A-V RELATIONSHIP depth based on full lake depth + ! 1 ---> CONIC shape with depth based on full lake depth + ! 2 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 3 ---> H-A-V RELATIONSHIP based on full lake depth !! ---------------------------- local variables real :: Htmp ! Temporary assign of incoming lake el. (m) @@ -273,8 +271,11 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa real :: tmp1, tmp2 real :: dh, dh1, dh2, dh3 ! Depth in weir, and height function for 3 order RK real :: It, Itdt_3, Itdt_2_3 ! inflow hydrographs - real :: maxWeirDepth !maximum capacity of weir - !real :: hdiff_vol, qdiff_vol ! water balance check variables + real :: maxWeirDepth ! maximum capacity of weir + real :: hin ! initial depth in storage from lake bottom (m) + real :: Dactivepool ! active depth in storage (m) + real :: LkBottom ! lake bottom elevation ASL (m) + !real :: hdiff_vol, qdiff_vol ! water balance check variables !! ---------------------------- subroutine body: from chow, mad mays. pg. 252 !! -- determine from inflow hydrograph @@ -294,14 +295,12 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa else ! If-block for Chow et al level pool scheme - ! Define reservoir max depth and rescale orifice and weir elevation - ! if using lake_shape_option >= 1: - if (lake_shape_option == 0) then ! define max depth based on active reservoir depth (NWM v1.x - v2.1) - ! ADD ADJUSTENTS FOR LAKE_SHAPE_OPT = 0 HERE! - else ! define max depth based on full lake depth - ! ADD ADJUSTENTS FOR LAKE_SHAPE_OPT > 0 HERE! - endif + ! Water depth in reservoir from lake bottom: + Dactivepool = (maxh - oe) * 3. + LkBottom = maxh - Dactivepool + hin = maxdepth - Dactivepool + (H - LkBottom) + ! Inflow calculations: It = qi0 Itdt_3 = qi0 + ((qi1 + ql - qi0) * 0.33) Itdt_2_3 = qi0 + ((qi1 + ql - qi0) * 0.67) @@ -334,7 +333,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0) then !dh1 = ((It - discharge)/sap)*dt - call STAGE_STORAGE(dh1,maxdepth,H,It,discharge,sap,vtp,dt,lake_shape_option,polya,polya) + call STAGE_STORAGE(dh1,maxdepth,hin,It,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh1 = 0.0 endif @@ -363,7 +362,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0.0) then !dh2 = ((Itdt_3 - discharge)/sap)*dt - call STAGE_STORAGE(dh2,maxdepth,H,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polya) + call STAGE_STORAGE(dh2,maxdepth,hin,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh2 = 0.0 endif @@ -391,7 +390,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0.0) then !dh3 = ((Itdt_2_3 - discharge)/sap)*dt - call STAGE_STORAGE(dh3,maxdepth,H,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polya) + call STAGE_STORAGE(dh3,maxdepth,hin,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh3 = 0.0 endif @@ -439,9 +438,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa !24 format('ofonl H dh sap Q ',f8.4,2x,f8.4,2x,f8.0,2x,f8.2) endif ! ENDIF for lake_shape_option.... - !print *, "TEST STAGE_STORAGE..." - !print *, dh, dh1, dh2, dh3 - !print *, dh, dh1_test, dh2_test, dh3_test + return ! ---------------------------------------------------------------- @@ -453,7 +450,7 @@ end subroutine LEVELPOOL_PHYSICS ! ------------------------------------------------ ! This subroutine calculates change in sotrage water level (Deltah) given the initial water level (h1), ! reservoir shape option (lake_shape_option), inflow (Inflow_dt), and discharge (res_discharge), and surface area (surf_area) - subroutine STAGE_STORAGE(Deltah,Dmax,h1,Inflow_dt,res_discharge,surf_area,total_vol, & + subroutine STAGE_STORAGE(Deltah,Dmax, h1,Inflow_dt,res_discharge,surf_area,total_vol, & Deltat,LAKE_SHAPE_OPT, polyhV_a, polyhV_b) !! ---------------------------- argument variables @@ -461,7 +458,6 @@ subroutine STAGE_STORAGE(Deltah,Dmax,h1,Inflow_dt,res_discharge,surf_area,total_ real, intent(OUT) :: Deltah ! change in sotrage water level (m) real, intent(IN) :: Dmax ! Maximum depth of reservoir (m) - real, intent(IN) :: h1 ! initial water level (m) real, intent(IN) :: Inflow_dt ! inflow (m^3/s) real, intent(IN) :: res_discharge ! storage discharge (m^3/s) real, intent(IN) :: surf_area ! surface area (m2) @@ -471,11 +467,9 @@ subroutine STAGE_STORAGE(Deltah,Dmax,h1,Inflow_dt,res_discharge,surf_area,total_ real, intent(IN) :: polyhV_b ! coefficient b in polynomial function V=a*h^b integer, intent(IN) :: LAKE_SHAPE_OPT ! LAKE_SHAPE_OPT (lake_shape_option) options from namelist: ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) - ! 1 ---> BOX shape with depth based on full lake depth - ! 2 ---> CONIC shape with depth based on full lake depth - ! 3 ---> RECTANGLUAR PRISM shape with depth based on full lake depth - ! 4 ---> ELLIPSOID shape with depth based on full lake depth - ! 5 ---> H-A-V RELATIONSHIP depth based on full lake depth + ! 1 ---> CONIC shape with depth based on full lake depth + ! 2 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 3 ---> H-A-V RELATIONSHIP based on full lake depth !! ---------------------------- local variables real, parameter :: pi = 3.14159265 @@ -483,38 +477,35 @@ subroutine STAGE_STORAGE(Deltah,Dmax,h1,Inflow_dt,res_discharge,surf_area,total_ real :: S1 ! initial reservoir storage (m3) real :: S2 ! reservoir storage after change (m3) real :: dS ! change in reservoir storage (m3) - real :: h2 ! cone base radius (m) + real :: h1 ! initial depth in storage (m) + real :: h2 ! final depth in storage (m) !! ---------------------------- ! change in storage (m3) dS = (Inflow_dt - res_discharge)*Deltat - if ( (LAKE_SHAPE_OPT == 0) .OR. (LAKE_SHAPE_OPT == 1) ) then ! BOX shape + if (LAKE_SHAPE_OPT == 0) then ! BOX shape h2 = (dS / surf_area) + h1 - elseif (LAKE_SHAPE_OPT == 2) then ! CONIC shape - R_cone = (surf_area/pi)**0.5 ! cone base radius (m) - S1 = (pi/3.) * (R_cone/Dmax)**2. * (h1**3.) ! storage before change in volume (m3) - S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) - h2 = (3.*S2/pi*(Dmax/R_cone)**2.)**(1./3.) ! elevation after change in storage in (m) - - elseif (LAKE_SHAPE_OPT == 3) then ! RECTANGLUAR PRISM shape - S1 = (surf_area/2./Dmax)*(h1**2.) ! storage before change in volume (m3) + elseif (LAKE_SHAPE_OPT == 1) then ! CONIC shape + !R_cone = (surf_area/pi)**0.5 ! cone base radius (m) + R_cone = (3.*total_vol/(pi*Dmax))**0.5 ! cone base radius (m) + S1 = (pi/3.) * (h1**3.) * (R_cone/Dmax)**2. ! storage before change in volume (m3) S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) - h2 = (2.*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) + h2 = (3.*S2/pi*(Dmax/R_cone)**2.)**(1./3.) ! elevation after change in storage in (m) - elseif (LAKE_SHAPE_OPT == 4) then ! ELLIPSOID shape - !S1 = (surf_area/2/Dmax)*(h1**2) ! storage before change in volume (m3) - !S2 = S1 + dS ! storage after change in volume (m3) - !h2 = (2*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) + elseif (LAKE_SHAPE_OPT == 2) then ! RECTANGLUAR PRISM shape + S1 = (surf_area/2./Dmax)*(h1**2.) ! storage before change in volume (m3) + S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) + h2 = (2.*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) - elseif (LAKE_SHAPE_OPT == 5) then ! H-A-V RELATIONSHIP - S1 = polyhV_a*(h1**polyhV_b) * 1.0E9 ! storage before change in volume (m3) - S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) - h2 = (S2/1.0E9/polyhV_a)**(1./polyhV_a) ! elevation after change in storage in (m) + elseif (LAKE_SHAPE_OPT == 3) then ! H-A-V RELATIONSHIP + !S1 = 1000000000.0*(1.5175592E-08)*(h1**3.516952) ! storage before change in volume (m3) + S1 = 1.0E9*polyhV_a*(h1**polyhV_b) + S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) + h2 = (S2/1.0E9/polyhV_a)**(1./polyhV_b) ! elevation after change in storage in (m) endif - Deltah = h2-h1 ! change in sotrage water level (m) return ! ---------------------------------------------------------------- diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index b5f318d92..b3ec07845 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -8183,6 +8183,7 @@ subroutine read_routelink(& n_CC, LAKEIDA, HRZAREA, & LAKEMAXH, WEIRH, WEIRC, WEIRL, DAML, & ORIFICEC, ORIFICEA, ORIFICEE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & reservoir_type_specified, reservoir_type, & reservoir_parameter_file, LATLAKE, & LONLAKE, ELEVLAKE, LAKEIDM, LAKEIDX, & @@ -8202,6 +8203,7 @@ subroutine read_routelink(& real, intent(INOUT), dimension(:) :: HRZAREA real, intent(INOUT), dimension(:) :: LAKEMAXH, WEIRH, WEIRC, WEIRL, DAML real, intent(INOUT), dimension(:) :: ORIFICEC, ORIFICEA, ORIFICEE +real, intent(INOUT), dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB logical, intent(IN) :: reservoir_type_specified integer, intent(INOUT), dimension(:) :: reservoir_type character(len=*), intent(in) :: reservoir_parameter_file @@ -8229,7 +8231,8 @@ subroutine read_routelink(& MannN, So, ChSSlp, Bw, Tw, Tw_CC, n_CC, LAKEIDA, HRZAREA, & LAKEMAXH, WEIRH, WEIRC, WEIRL, DAML, ORIFICEC, & ORIFICEA, ORIFICEE, reservoir_type_specified, reservoir_type, reservoir_parameter_file, & - gages, gageMiss, LAKEIDM, NLAKES, latlake, lonlake,ELEVLAKE) + gages, gageMiss, LAKEIDM, NLAKES, latlake, lonlake,ELEVLAKE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB) !--- get the lake configuration here. #ifdef MPP_LAND @@ -8254,6 +8257,10 @@ subroutine read_routelink(& call mpp_land_bcast_real(NLAKES,ORIFICEC) call mpp_land_bcast_real(NLAKES,ORIFICEA) call mpp_land_bcast_real(NLAKES,ORIFICEE) + call mpp_land_bcast_real(NLAKES,LAKEMAXDEPTH) + call mpp_land_bcast_real(NLAKES,LAKEVOL) + call mpp_land_bcast_real(NLAKES,LAKEPOLYA) + call mpp_land_bcast_real(NLAKES,LAKEPOLYB) call mpp_land_bcast_real(NLAKES,LATLAKE ) call mpp_land_bcast_real(NLAKES,LONLAKE ) call mpp_land_bcast_real(NLAKES,ELEVLAKE) @@ -8271,7 +8278,8 @@ subroutine readLinkSL( GNLINKSL,NLINKSL,route_link_f, route_lake_f, maxorder, & MannN, So, ChSSlp, Bw, Tw, Tw_CC, n_CC, LAKEIDA, HRZAREA, & LAKEMAXH,WEIRH, WEIRC, WEIRL, DAML, ORIFICEC, & ORIFICEA, ORIFICEE, reservoir_type_specified, reservoir_type, reservoir_parameter_file, & - gages, gageMiss, LAKEIDM,NLAKES, latlake, lonlake, ELEVLAKE) + gages, gageMiss, LAKEIDM,NLAKES, latlake, lonlake, ELEVLAKE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB) implicit none character(len=*) :: route_link_f,route_lake_f @@ -8293,7 +8301,7 @@ subroutine readLinkSL( GNLINKSL,NLINKSL,route_link_f, route_lake_f, maxorder, & character(len=*), intent(in) :: reservoir_parameter_file REAL, intent(out), dimension(:) :: HRZAREA,LAKEMAXH, WEIRC, WEIRL, DAML, ORIFICEC, WEIRH, & ORIFICEA, ORIFICEE, ELEVLAKE - real, allocatable, dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB + real, intent(out), dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB !end NLAKES INTEGER, dimension(GNLINKSL) :: tmpLAKEIDA, tmpLINKID, tmpTO_NODE, tmpTYPEL, tmpORDER @@ -8438,6 +8446,10 @@ subroutine readLinkSL( GNLINKSL,NLINKSL,route_link_f, route_lake_f, maxorder, & call mpp_land_bcast_real(NLAKES, ORIFICEC) call mpp_land_bcast_real(NLAKES, ORIFICEA) call mpp_land_bcast_real(NLAKES, ORIFICEE) + call mpp_land_bcast_real(NLAKES, LAKEMAXDEPTH) + call mpp_land_bcast_real(NLAKES, LAKEVOL) + call mpp_land_bcast_real(NLAKES, LAKEPOLYA) + call mpp_land_bcast_real(NLAKES, LAKEPOLYB) call mpp_land_bcast_int(NLAKES, LAKEIDM) call mpp_land_bcast_real(NLAKES, ELEVLAKE) call mpp_land_bcast_int(NLAKES, reservoir_type) @@ -8919,7 +8931,7 @@ subroutine read_route_lake_netcdf(route_lake_file, & ORIFICEC, ORIFICEA, ORIFICEE, reservoir_type_specified, & reservoir_type, reservoir_parameter_file, & LAKEIDM, lakelat, lakelon, & - LAKEMAXDEPTH, LAKEVOL, LAKE_POLYA, LAKE_POLYB, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & ELEVLAKE, NLAKES) implicit none @@ -8930,7 +8942,7 @@ subroutine read_route_lake_netcdf(route_lake_file, & integer, dimension(:), intent(out) :: LAKEIDM real, dimension(:), intent(out) :: HRZAREA, LAKEMAXH, WEIRC, WEIRL, WEIRH, DAML real, dimension(:), intent(out) :: ORIFICEC, ORIFICEA, ORIFICEE, lakelat, lakelon - real, dimension(:), intent(out) :: LAKEVOL, LAKEMAXDEPTH, LAKE_POLYA, LAKE_POLYB + real, dimension(:), intent(out) :: LAKEVOL, LAKEMAXDEPTH, LAKEPOLYA, LAKEPOLYB real, dimension(:), intent(out) :: ELEVLAKE integer, dimension(:), intent(out) :: reservoir_type @@ -8964,8 +8976,8 @@ subroutine read_route_lake_netcdf(route_lake_file, & call get_1d_netcdf_real(ncid, 'lon', lakelon, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'poly_a', LAKE_POLYA, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'poly_b', LAKE_POLYB, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_b', LAKEPOLYB, 'read_route_lake_netcdf', .TRUE.) !remove the alt var. and add initial fractional depth var. LKR/DY call get_1d_netcdf_real(ncid, 'ifd', ELEVLAKE, 'read_route_lake_netcdf', .FALSE.) @@ -8987,6 +8999,8 @@ subroutine read_route_lake_netcdf(route_lake_file, & print*,'HRZAREA', HRZAREA(ii) print*,'LAKEMAXH', LAKEMAXH(ii), 'WEIRC', WEIRC(ii), 'WEIRL', WEIRL(ii), 'DAML', DAML(ii) print*,'ORIFICEC', ORIFICEC(ii), 'ORIFICEA', ORIFICEA(ii), 'ORIFICEE', ORIFICEE(ii) + print*,'LAKEMAXDEPTH', LAKEMAXDEPTH(ii), 'LAKEVOL', LAKEVOL(ii) + print*,'LAKEPOLYA', LAKEPOLYA(ii), 'LAKEPOLYB', LAKEPOLYB(ii) print*,"finish read_route_lake_netcdf" #endif diff --git a/trunk/NDHMS/Routing/module_RT.F b/trunk/NDHMS/Routing/module_RT.F index 5bb765e3b..caa9848f8 100644 --- a/trunk/NDHMS/Routing/module_RT.F +++ b/trunk/NDHMS/Routing/module_RT.F @@ -787,6 +787,8 @@ subroutine LandRT_ini(did) rt_domain(did)%DAML, & rt_domain(did)%ORIFICEC, rt_domain(did)%ORIFICEA, & rt_domain(did)%ORIFICEE, & + rt_domain(did)%LAKEMAXDEPTH, rt_domain(did)%LAKEVOL, & + rt_domain(did)%LAKEPOLYA, rt_domain(did)%LAKEPOLYB, & nlst(did)%reservoir_type_specified, & rt_domain(did)%reservoir_type, & nlst(did)%reservoir_parameter_file, & diff --git a/trunk/NDHMS/Run/hydro.namelist b/trunk/NDHMS/Run/hydro.namelist new file mode 100644 index 000000000..f1f2c0b30 --- /dev/null +++ b/trunk/NDHMS/Run/hydro.namelist @@ -0,0 +1,265 @@ +&HYDRO_nlist +!!!! ---------------------- SYSTEM COUPLING ----------------------- !!!! + +! Specify what is being coupled: 1=HRLDAS (offline Noah-LSM), 2=WRF, 3=NASA/LIS, 4=CLM +sys_cpl = 1 + +!!!! ------------------- MODEL INPUT DATA FILES ------------------- !!!! + +! Specify land surface model gridded input data file (e.g.: "geo_em.d01.nc") +GEO_STATIC_FLNM = "./DOMAIN/geo_em.d01.nc" + +! Specify the high-resolution routing terrain input data file (e.g.: "Fulldom_hires.nc") +GEO_FINEGRID_FLNM = "./DOMAIN/Fulldom_hires.nc" + +! Specify the spatial hydro parameters file (e.g.: "hydro2dtbl.nc") +! If you specify a filename and the file does not exist, it will be created for you. +HYDROTBL_F = "./DOMAIN/hydro2dtbl.nc" + +! Specify spatial metadata file for land surface grid. (e.g.: "GEOGRID_LDASOUT_Spatial_Metadata.nc") +LAND_SPATIAL_META_FLNM = "./DOMAIN/GEOGRID_LDASOUT_Spatial_Metadata.nc" + +! Specify the name of the restart file if starting from restart...comment out with '!' if not... +RESTART_FILE = 'RESTART/HYDRO_RST.2011-08-26_00:00_DOMAIN1' + +!!!! --------------------- MODEL SETUP OPTIONS -------------------- !!!! + +! Specify the domain or nest number identifier...(integer) +IGRID = 1 + +! Specify the restart file write frequency...(minutes) +! A value of -99999 will output restarts on the first day of the month only. +rst_dt = 120 + +! Reset the LSM soil states from the high-res routing restart file (1=overwrite, 0=no overwrite) +! NOTE: Only turn this option on if overland or subsurface rotuing is active! +rst_typ = 1 + +! Restart file format control +rst_bi_in = 0 !0: use netcdf input restart file (default) + !1: use parallel io for reading multiple restart files, 1 per core +rst_bi_out = 0 !0: use netcdf output restart file (default) + !1: use parallel io for outputting multiple restart files, 1 per core + +! Restart switch to set restart accumulation variables to 0 (0=no reset, 1=yes reset to 0.0) +RSTRT_SWC = 0 + +! Specify baseflow/bucket model initialization...(0=cold start from table, 1=restart file) +GW_RESTART = 1 + +!!!! -------------------- MODEL OUTPUT CONTROL -------------------- !!!! + +! Specify the output file write frequency...(minutes) +out_dt = 60 + +! Specify the number of output times to be contained within each output history file...(integer) +! SET = 1 WHEN RUNNING CHANNEL ROUTING ONLY/CALIBRATION SIMS!!! +! SET = 1 WHEN RUNNING COUPLED TO WRF!!! +SPLIT_OUTPUT_COUNT = 1 + +! Specify the minimum stream order to output to netcdf point file...(integer) +! Note: lower value of stream order produces more output. +order_to_write = 1 + +! Flag to turn on/off new I/O routines: 0 = deprecated output routines (use when running with Noah LSM), +! 1 = with scale/offset/compression, ! 2 = with scale/offset/NO compression, +! 3 = compression only, 4 = no scale/offset/compression (default) +io_form_outputs = 4 + +! Realtime run configuration option: +! 0=all (default), 1=analysis, 2=short-range, 3=medium-range, 4=long-range, 5=retrospective, +! 6=diagnostic (includes all of 1-4 outputs combined) +io_config_outputs = 0 + +! Option to write output files at time 0 (restart cold start time): 0=no, 1=yes (default) +t0OutputFlag = 1 + +! Options to output channel & bucket influxes. Only active for UDMP_OPT=1. +! Nonzero choice requires that out_dt above matches NOAH_TIMESTEP in namelist.hrldas. +! 0=None (default), 1=channel influxes (qSfcLatRunoff, qBucket) +! 2=channel+bucket fluxes (qSfcLatRunoff, qBucket, qBtmVertRunoff_toBucket) +! 3=channel accumulations (accSfcLatRunoff, accBucket) *** NOT TESTED *** +output_channelBucket_influx = 0 + +! Output netcdf file control +CHRTOUT_DOMAIN = 1 ! Netcdf point timeseries output at all channel points (1d) + ! 0 = no output, 1 = output +CHANOBS_DOMAIN = 0 ! Netcdf point timeseries at forecast points or gage points (defined in Routelink) + ! 0 = no output, 1 = output at forecast points or gage points. +CHRTOUT_GRID = 0 ! Netcdf grid of channel streamflow values (2d) + ! 0 = no output, 1 = output + ! NOTE: Not available with reach-based routing +LSMOUT_DOMAIN = 0 ! Netcdf grid of variables passed between LSM and routing components (2d) + ! 0 = no output, 1 = output + ! NOTE: No scale_factor/add_offset available +RTOUT_DOMAIN = 1 ! Netcdf grid of terrain routing variables on routing grid (2d) + ! 0 = no output, 1 = output +output_gw = 1 ! Netcdf GW output + ! 0 = no output, 1 = output +outlake = 1 ! Netcdf grid of lake values (1d) + ! 0 = no output, 1 = output +frxst_pts_out = 0 ! ASCII text file of forecast points or gage points (defined in Routelink) + ! 0 = no output, 1 = output + +!!!! ------------ PHYSICS OPTIONS AND RELATED SETTINGS ------------ !!!! + +! Specify the number of soil layers (integer) and the depth of the bottom of each layer... (meters) +! Notes: In Version 1 of WRF-Hydro these must be the same as in the namelist.input file. +! Future versions will permit this to be different. +NSOIL=4 +ZSOIL8(1) = -0.10 +ZSOIL8(2) = -0.40 +ZSOIL8(3) = -1.00 +ZSOIL8(4) = -2.00 + +! Specify the grid spacing of the terrain routing grid...(meters) +DXRT = 250.0 + +! Specify the integer multiple between the land model grid and the terrain routing grid...(integer) +AGGFACTRT = 4 + +! Specify the channel routing model timestep...(seconds) +DTRT_CH = 10 + +! Specify the terrain routing model timestep...(seconds) +DTRT_TER = 10 + +! Switch to activate subsurface routing...(0=no, 1=yes) +SUBRTSWCRT = 1 + +! Switch to activate surface overland flow routing...(0=no, 1=yes) +OVRTSWCRT = 1 + +! Specify overland flow routing option: 1=Seepest Descent (D8) 2=CASC2D (not active) +! NOTE: Currently subsurface flow is only steepest descent +rt_option = 1 + +! Switch to activate channel routing...(0=no, 1=yes) +CHANRTSWCRT = 1 + +! Specify channel routing option: 1=Muskingam-reach, 2=Musk.-Cunge-reach, 3=Diff.Wave-gridded +channel_option = 3 + +! Specify the reach file for reach-based routing options (e.g.: "Route_Link.nc") +!route_link_f = "./DOMAIN/Route_Link.nc" + +! If using channel_option=2, activate the compound channel formulation? (Default=.FALSE.) +! This option is currently only supported if using reach-based routing with UDMP=1. +compound_channel = .FALSE. + +! Specify the lake parameter file (e.g.: "LAKEPARM.nc"). +! Note REQUIRED if lakes are on. +route_lake_f = "./DOMAIN/LAKEPARM.nc" + +! Specify the reservoir parameter file +reservoir_parameter_file = "./DOMAIN/persistence_parm.nc" + +! If using USGS persistence reservoirs, set to True. (default=.FALSE.) +reservoir_persistence_usgs = .FALSE. + +! Specify the path to the timeslice files to be used by USGS reservoirs +reservoir_usgs_timeslice_path = "./usgs_timeslices/" + +! If using USACE persistence reservoirs, set to True. (default=.FALSE.) +reservoir_persistence_usace = .FALSE. + +! Specify the path to the timeslice files to be used by USACE reservoirs +reservoir_usace_timeslice_path = "./usace_timeslices/" + +! Specify lookback hours to read reservoir observation data +reservoir_observation_lookback_hours = 24 + +! Specify update time interval in seconds to read new reservoir observation data +! The default is 86400 (seconds per day). Set to 3600 for standard and extended AnA simulations. +! Set to 1000000000 for short range and medium range forecasts. +reservoir_observation_update_time_interval_seconds = 3600 + +! If using RFC forecast reservoirs, set to True. (default=.FALSE.) +reservoir_rfc_forecasts = .FALSE. + +! Specify the path to the RFC time series files to be used by reservoirs +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 + +! Path to the "timeslice" observation files. +timeSlicePath = "./nudgingTimeSliceObs/" + +nudgingParamFile = "DOMAIN/nudgingParams.nc" + +! Nudging restart file = "nudgingLastObsFile" +! nudgingLastObsFile defaults to '', which will look for nudgingLastObs.YYYY-mm-dd_HH:MM:SS.nc +! **AT THE INITALIZATION TIME OF THE RUN**. Set to a missing file to use no restart. +!nudgingLastObsFile = '/a/nonexistent/file/gives/nudging/cold/start' + +!! Parallel input of nudging timeslice observation files? +readTimesliceParallel = .TRUE. + +! temporalPersistence defaults to true, only runs if necessary params present. +temporalPersistence = .FALSE. + +! The total number of last (obs, modeled) pairs to save in nudgingLastObs for +! removal of bias. This is the maximum array length. (This option is active when persistBias=FALSE) +! (Default=960=10days @15min obs resolution, if all the obs are present and longer if not.) +nLastObs = 960 + +! If using temporalPersistence the last observation persists by default. +! This option instead persists the bias after the last observation. +persistBias = .FALSE. + +! AnA (FALSE) vs Forecast (TRUE) bias persistence. +! If persistBias: Does the window for calculating the bias end at +! model init time (=t0)? +! FALSE = window ends at model time (moving), +! TRUE = window ends at init=t0(fcst) time. +! (If commented out, Default=FALSE) +! Note: Perfect restart tests require this option to be .FALSE. +biasWindowBeforeT0 = .FALSE. + +! If persistBias: Only use this many last (obs, modeled) pairs. (If Commented out, Default=-1*nLastObs) +! > 0: apply an age-based filter, units=hours. +! = 0: apply no additional filter, use all available/usable obs. +! < 0: apply an count-based filter, units=count +maxAgePairsBiasPersist = -960 + +! If persistBias: The minimum number of last (obs, modeled) pairs, with age less than +! maxAgePairsBiasPersist, required to apply a bias correction. (default=8) +minNumPairsBiasPersist = 8 + +! If persistBias: give more weight to observations closer in time? (default=FALSE) +invDistTimeWeightBias = .TRUE. + +! If persistBias: "No constructive interference in bias correction?", Reduce the bias adjustment +! when the model and the bias adjustment have the same sign relative to the modeled flow at t0? +! (default=FALSE) +! Note: Perfect restart tests require this option to be .FALSE. +noConstInterfBias = .FALSE. + +/ From 0a21b45b8b95d9798cc46ae1e9b5b6af93a6d47d Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Fri, 24 Jul 2020 13:20:17 -0600 Subject: [PATCH 03/12] fix minor issue in levelpool code --- trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F | 2 ++ 1 file changed, 2 insertions(+) diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index 52e96d2ba..a1b6bb8f2 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -506,6 +506,8 @@ subroutine STAGE_STORAGE(Deltah,Dmax, h1,Inflow_dt,res_discharge,surf_area,total h2 = (S2/1.0E9/polyhV_a)**(1./polyhV_b) ! elevation after change in storage in (m) endif + Deltah = h2-h1 ! change in sotrage water level (m) + return ! ---------------------------------------------------------------- From 314369b99596aefe4d94c61d707f8a8d36470df0 Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Mon, 27 Jul 2020 16:46:49 -0600 Subject: [PATCH 04/12] Add default value of 0 for lake_shape_option --- trunk/NDHMS/OrchestratorLayer/config.f90 | 5 +- trunk/NDHMS/Run/hydro.namelist | 265 ----------------------- 2 files changed, 3 insertions(+), 267 deletions(-) delete mode 100644 trunk/NDHMS/Run/hydro.namelist diff --git a/trunk/NDHMS/OrchestratorLayer/config.f90 b/trunk/NDHMS/OrchestratorLayer/config.f90 index 324cd2485..1a5ba7e7e 100644 --- a/trunk/NDHMS/OrchestratorLayer/config.f90 +++ b/trunk/NDHMS/OrchestratorLayer/config.f90 @@ -465,8 +465,8 @@ subroutine init_namelist_rt_field(did) GWBASESWCRT, GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, & sys_cpl, rst_typ, rst_bi_in, rst_bi_out, & gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, gwsoilcpl, & - UDMP_OPT, io_form_outputs, bucket_loss, & - lake_shape_option + UDMP_OPT, io_form_outputs, bucket_loss + integer :: lake_shape_option = 0 real:: DTRT_TER,DTRT_CH,dxrt, gwChanCondConstIn, gwChanCondConstOut, gwIhShift character(len=256) :: route_topo_f="" character(len=256) :: route_chan_f="" @@ -592,6 +592,7 @@ subroutine init_namelist_rt_field(did) reservoir_rfc_forecasts = .FALSE. reservoir_rfc_forecasts_lookback_hours = 24 reservoir_type_specified = .FALSE. + lake_shape_option = 0 #ifdef WRF_HYDRO_NUDGING ! Default values for NUDGING_nlist diff --git a/trunk/NDHMS/Run/hydro.namelist b/trunk/NDHMS/Run/hydro.namelist deleted file mode 100644 index f1f2c0b30..000000000 --- a/trunk/NDHMS/Run/hydro.namelist +++ /dev/null @@ -1,265 +0,0 @@ -&HYDRO_nlist -!!!! ---------------------- SYSTEM COUPLING ----------------------- !!!! - -! Specify what is being coupled: 1=HRLDAS (offline Noah-LSM), 2=WRF, 3=NASA/LIS, 4=CLM -sys_cpl = 1 - -!!!! ------------------- MODEL INPUT DATA FILES ------------------- !!!! - -! Specify land surface model gridded input data file (e.g.: "geo_em.d01.nc") -GEO_STATIC_FLNM = "./DOMAIN/geo_em.d01.nc" - -! Specify the high-resolution routing terrain input data file (e.g.: "Fulldom_hires.nc") -GEO_FINEGRID_FLNM = "./DOMAIN/Fulldom_hires.nc" - -! Specify the spatial hydro parameters file (e.g.: "hydro2dtbl.nc") -! If you specify a filename and the file does not exist, it will be created for you. -HYDROTBL_F = "./DOMAIN/hydro2dtbl.nc" - -! Specify spatial metadata file for land surface grid. (e.g.: "GEOGRID_LDASOUT_Spatial_Metadata.nc") -LAND_SPATIAL_META_FLNM = "./DOMAIN/GEOGRID_LDASOUT_Spatial_Metadata.nc" - -! Specify the name of the restart file if starting from restart...comment out with '!' if not... -RESTART_FILE = 'RESTART/HYDRO_RST.2011-08-26_00:00_DOMAIN1' - -!!!! --------------------- MODEL SETUP OPTIONS -------------------- !!!! - -! Specify the domain or nest number identifier...(integer) -IGRID = 1 - -! Specify the restart file write frequency...(minutes) -! A value of -99999 will output restarts on the first day of the month only. -rst_dt = 120 - -! Reset the LSM soil states from the high-res routing restart file (1=overwrite, 0=no overwrite) -! NOTE: Only turn this option on if overland or subsurface rotuing is active! -rst_typ = 1 - -! Restart file format control -rst_bi_in = 0 !0: use netcdf input restart file (default) - !1: use parallel io for reading multiple restart files, 1 per core -rst_bi_out = 0 !0: use netcdf output restart file (default) - !1: use parallel io for outputting multiple restart files, 1 per core - -! Restart switch to set restart accumulation variables to 0 (0=no reset, 1=yes reset to 0.0) -RSTRT_SWC = 0 - -! Specify baseflow/bucket model initialization...(0=cold start from table, 1=restart file) -GW_RESTART = 1 - -!!!! -------------------- MODEL OUTPUT CONTROL -------------------- !!!! - -! Specify the output file write frequency...(minutes) -out_dt = 60 - -! Specify the number of output times to be contained within each output history file...(integer) -! SET = 1 WHEN RUNNING CHANNEL ROUTING ONLY/CALIBRATION SIMS!!! -! SET = 1 WHEN RUNNING COUPLED TO WRF!!! -SPLIT_OUTPUT_COUNT = 1 - -! Specify the minimum stream order to output to netcdf point file...(integer) -! Note: lower value of stream order produces more output. -order_to_write = 1 - -! Flag to turn on/off new I/O routines: 0 = deprecated output routines (use when running with Noah LSM), -! 1 = with scale/offset/compression, ! 2 = with scale/offset/NO compression, -! 3 = compression only, 4 = no scale/offset/compression (default) -io_form_outputs = 4 - -! Realtime run configuration option: -! 0=all (default), 1=analysis, 2=short-range, 3=medium-range, 4=long-range, 5=retrospective, -! 6=diagnostic (includes all of 1-4 outputs combined) -io_config_outputs = 0 - -! Option to write output files at time 0 (restart cold start time): 0=no, 1=yes (default) -t0OutputFlag = 1 - -! Options to output channel & bucket influxes. Only active for UDMP_OPT=1. -! Nonzero choice requires that out_dt above matches NOAH_TIMESTEP in namelist.hrldas. -! 0=None (default), 1=channel influxes (qSfcLatRunoff, qBucket) -! 2=channel+bucket fluxes (qSfcLatRunoff, qBucket, qBtmVertRunoff_toBucket) -! 3=channel accumulations (accSfcLatRunoff, accBucket) *** NOT TESTED *** -output_channelBucket_influx = 0 - -! Output netcdf file control -CHRTOUT_DOMAIN = 1 ! Netcdf point timeseries output at all channel points (1d) - ! 0 = no output, 1 = output -CHANOBS_DOMAIN = 0 ! Netcdf point timeseries at forecast points or gage points (defined in Routelink) - ! 0 = no output, 1 = output at forecast points or gage points. -CHRTOUT_GRID = 0 ! Netcdf grid of channel streamflow values (2d) - ! 0 = no output, 1 = output - ! NOTE: Not available with reach-based routing -LSMOUT_DOMAIN = 0 ! Netcdf grid of variables passed between LSM and routing components (2d) - ! 0 = no output, 1 = output - ! NOTE: No scale_factor/add_offset available -RTOUT_DOMAIN = 1 ! Netcdf grid of terrain routing variables on routing grid (2d) - ! 0 = no output, 1 = output -output_gw = 1 ! Netcdf GW output - ! 0 = no output, 1 = output -outlake = 1 ! Netcdf grid of lake values (1d) - ! 0 = no output, 1 = output -frxst_pts_out = 0 ! ASCII text file of forecast points or gage points (defined in Routelink) - ! 0 = no output, 1 = output - -!!!! ------------ PHYSICS OPTIONS AND RELATED SETTINGS ------------ !!!! - -! Specify the number of soil layers (integer) and the depth of the bottom of each layer... (meters) -! Notes: In Version 1 of WRF-Hydro these must be the same as in the namelist.input file. -! Future versions will permit this to be different. -NSOIL=4 -ZSOIL8(1) = -0.10 -ZSOIL8(2) = -0.40 -ZSOIL8(3) = -1.00 -ZSOIL8(4) = -2.00 - -! Specify the grid spacing of the terrain routing grid...(meters) -DXRT = 250.0 - -! Specify the integer multiple between the land model grid and the terrain routing grid...(integer) -AGGFACTRT = 4 - -! Specify the channel routing model timestep...(seconds) -DTRT_CH = 10 - -! Specify the terrain routing model timestep...(seconds) -DTRT_TER = 10 - -! Switch to activate subsurface routing...(0=no, 1=yes) -SUBRTSWCRT = 1 - -! Switch to activate surface overland flow routing...(0=no, 1=yes) -OVRTSWCRT = 1 - -! Specify overland flow routing option: 1=Seepest Descent (D8) 2=CASC2D (not active) -! NOTE: Currently subsurface flow is only steepest descent -rt_option = 1 - -! Switch to activate channel routing...(0=no, 1=yes) -CHANRTSWCRT = 1 - -! Specify channel routing option: 1=Muskingam-reach, 2=Musk.-Cunge-reach, 3=Diff.Wave-gridded -channel_option = 3 - -! Specify the reach file for reach-based routing options (e.g.: "Route_Link.nc") -!route_link_f = "./DOMAIN/Route_Link.nc" - -! If using channel_option=2, activate the compound channel formulation? (Default=.FALSE.) -! This option is currently only supported if using reach-based routing with UDMP=1. -compound_channel = .FALSE. - -! Specify the lake parameter file (e.g.: "LAKEPARM.nc"). -! Note REQUIRED if lakes are on. -route_lake_f = "./DOMAIN/LAKEPARM.nc" - -! Specify the reservoir parameter file -reservoir_parameter_file = "./DOMAIN/persistence_parm.nc" - -! If using USGS persistence reservoirs, set to True. (default=.FALSE.) -reservoir_persistence_usgs = .FALSE. - -! Specify the path to the timeslice files to be used by USGS reservoirs -reservoir_usgs_timeslice_path = "./usgs_timeslices/" - -! If using USACE persistence reservoirs, set to True. (default=.FALSE.) -reservoir_persistence_usace = .FALSE. - -! Specify the path to the timeslice files to be used by USACE reservoirs -reservoir_usace_timeslice_path = "./usace_timeslices/" - -! Specify lookback hours to read reservoir observation data -reservoir_observation_lookback_hours = 24 - -! Specify update time interval in seconds to read new reservoir observation data -! The default is 86400 (seconds per day). Set to 3600 for standard and extended AnA simulations. -! Set to 1000000000 for short range and medium range forecasts. -reservoir_observation_update_time_interval_seconds = 3600 - -! If using RFC forecast reservoirs, set to True. (default=.FALSE.) -reservoir_rfc_forecasts = .FALSE. - -! Specify the path to the RFC time series files to be used by reservoirs -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 - -! Path to the "timeslice" observation files. -timeSlicePath = "./nudgingTimeSliceObs/" - -nudgingParamFile = "DOMAIN/nudgingParams.nc" - -! Nudging restart file = "nudgingLastObsFile" -! nudgingLastObsFile defaults to '', which will look for nudgingLastObs.YYYY-mm-dd_HH:MM:SS.nc -! **AT THE INITALIZATION TIME OF THE RUN**. Set to a missing file to use no restart. -!nudgingLastObsFile = '/a/nonexistent/file/gives/nudging/cold/start' - -!! Parallel input of nudging timeslice observation files? -readTimesliceParallel = .TRUE. - -! temporalPersistence defaults to true, only runs if necessary params present. -temporalPersistence = .FALSE. - -! The total number of last (obs, modeled) pairs to save in nudgingLastObs for -! removal of bias. This is the maximum array length. (This option is active when persistBias=FALSE) -! (Default=960=10days @15min obs resolution, if all the obs are present and longer if not.) -nLastObs = 960 - -! If using temporalPersistence the last observation persists by default. -! This option instead persists the bias after the last observation. -persistBias = .FALSE. - -! AnA (FALSE) vs Forecast (TRUE) bias persistence. -! If persistBias: Does the window for calculating the bias end at -! model init time (=t0)? -! FALSE = window ends at model time (moving), -! TRUE = window ends at init=t0(fcst) time. -! (If commented out, Default=FALSE) -! Note: Perfect restart tests require this option to be .FALSE. -biasWindowBeforeT0 = .FALSE. - -! If persistBias: Only use this many last (obs, modeled) pairs. (If Commented out, Default=-1*nLastObs) -! > 0: apply an age-based filter, units=hours. -! = 0: apply no additional filter, use all available/usable obs. -! < 0: apply an count-based filter, units=count -maxAgePairsBiasPersist = -960 - -! If persistBias: The minimum number of last (obs, modeled) pairs, with age less than -! maxAgePairsBiasPersist, required to apply a bias correction. (default=8) -minNumPairsBiasPersist = 8 - -! If persistBias: give more weight to observations closer in time? (default=FALSE) -invDistTimeWeightBias = .TRUE. - -! If persistBias: "No constructive interference in bias correction?", Reduce the bias adjustment -! when the model and the bias adjustment have the same sign relative to the modeled flow at t0? -! (default=FALSE) -! Note: Perfect restart tests require this option to be .FALSE. -noConstInterfBias = .FALSE. - -/ From 13f0f5849dec08cf047ea9288b56a6279d13dd7c Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Tue, 28 Jul 2020 15:41:49 -0600 Subject: [PATCH 05/12] Update the code to read from old LAKEPARM file when lake_shape_option is not defined or defined 0 --- trunk/NDHMS/Routing/module_HYDRO_io.F | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index b3ec07845..e392a7f81 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -8974,10 +8974,23 @@ subroutine read_route_lake_netcdf(route_lake_file, & call get_1d_netcdf_real(ncid, 'OrificeE', ORIFICEE, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'lat', lakelat, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'lon', lakelon, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.) - call get_1d_netcdf_real(ncid, 'poly_b', LAKEPOLYB, 'read_route_lake_netcdf', .TRUE.) + + ! read variables for lake_shape_option not equal to 0 (default Box shape) + if(nlst(did)%lake_shape_option /= 0) then + print*, "Different Lake Shape Option (rather than Default Box Shape) is defined in hydro.namelist." + print*, "Variables are all required in LAKEPARM file: LkVol, LkMxDepth, poly_a, poly_b!" + call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.) + call get_1d_netcdf_real(ncid, 'poly_b', LAKEPOLYB, 'read_route_lake_netcdf', .TRUE.) + else + print*, "Default Lake Shape Option (Box Shape) is defined in hydro.namelist" + LAKEVOL = 0.0 + LAKEMAXDEPTH = 0.0 + LAKEPOLYA = 0.0 + LAKEPOLYB = 0.0 + end if + !remove the alt var. and add initial fractional depth var. LKR/DY call get_1d_netcdf_real(ncid, 'ifd', ELEVLAKE, 'read_route_lake_netcdf', .FALSE.) From b2fdc1f38c7372606f7636a59e7d7ab76cfe03d1 Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Wed, 29 Jul 2020 13:37:28 -0600 Subject: [PATCH 06/12] Switched back dh calculations in LP to old scheme to check for diffs in engineering tests --- .../Routing/Reservoirs/Level_Pool/module_levelpool.F | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index a1b6bb8f2..be6dbf5de 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -332,8 +332,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0) then - !dh1 = ((It - discharge)/sap)*dt - call STAGE_STORAGE(dh1,maxdepth,hin,It,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + dh1 = ((It - discharge)/sap)*dt + !call STAGE_STORAGE(dh1,maxdepth,hin,It,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh1 = 0.0 endif @@ -361,8 +361,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0.0) then - !dh2 = ((Itdt_3 - discharge)/sap)*dt - call STAGE_STORAGE(dh2,maxdepth,hin,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + dh2 = ((Itdt_3 - discharge)/sap)*dt + !call STAGE_STORAGE(dh2,maxdepth,hin,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh2 = 0.0 endif @@ -389,8 +389,8 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0.0) then - !dh3 = ((Itdt_2_3 - discharge)/sap)*dt - call STAGE_STORAGE(dh3,maxdepth,hin,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + dh3 = ((Itdt_2_3 - discharge)/sap)*dt + !call STAGE_STORAGE(dh3,maxdepth,hin,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) else dh3 = 0.0 endif From d3815c453a24dae31e05813864b515d3f54a4de6 Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Wed, 29 Jul 2020 15:57:25 -0600 Subject: [PATCH 07/12] Update LP to debug issues related to calculation of dh in shape option 0 --- .../Reservoirs/Level_Pool/module_levelpool.F | 29 ++++++++++++------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index be6dbf5de..e3f13a84d 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -332,9 +332,12 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0) then - dh1 = ((It - discharge)/sap)*dt - !call STAGE_STORAGE(dh1,maxdepth,hin,It,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) - else + if (lake_shape_option == 0) then + dh1 = ((It - discharge)/sap)*dt + else + call STAGE_STORAGE(dh1,maxdepth,hin,It,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + endif + else dh1 = 0.0 endif @@ -361,8 +364,11 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa if (sap > 0.0) then - dh2 = ((Itdt_3 - discharge)/sap)*dt - !call STAGE_STORAGE(dh2,maxdepth,hin,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + if (lake_shape_option == 0) then + dh2 = ((Itdt_3 - discharge)/sap)*dt + else + call STAGE_STORAGE(dh2,maxdepth,hin,Itdt_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + endif else dh2 = 0.0 endif @@ -389,8 +395,11 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa endif if (sap > 0.0) then - dh3 = ((Itdt_2_3 - discharge)/sap)*dt - !call STAGE_STORAGE(dh3,maxdepth,hin,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + if (lake_shape_option == 0) then + dh3 = ((Itdt_2_3 - discharge)/sap)*dt + else + call STAGE_STORAGE(dh3,maxdepth,hin,Itdt_2_3,discharge,sap,vtp,dt,lake_shape_option,polya,polyb) + endif else dh3 = 0.0 endif @@ -484,10 +493,10 @@ subroutine STAGE_STORAGE(Deltah,Dmax, h1,Inflow_dt,res_discharge,surf_area,total ! change in storage (m3) dS = (Inflow_dt - res_discharge)*Deltat - if (LAKE_SHAPE_OPT == 0) then ! BOX shape - h2 = (dS / surf_area) + h1 + !if (LAKE_SHAPE_OPT == 0) then ! BOX shape + ! h2 = (dS / surf_area) + h1 - elseif (LAKE_SHAPE_OPT == 1) then ! CONIC shape + if (LAKE_SHAPE_OPT == 1) then ! CONIC shape !R_cone = (surf_area/pi)**0.5 ! cone base radius (m) R_cone = (3.*total_vol/(pi*Dmax))**0.5 ! cone base radius (m) S1 = (pi/3.) * (h1**3.) * (R_cone/Dmax)**2. ! storage before change in volume (m3) From a2e0636c53e7429b23106fb07c356e09b249282c Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Wed, 29 Jul 2020 20:15:50 -0600 Subject: [PATCH 08/12] Add Lake Shape arguments to READ_CHROUTING1 et al --- trunk/NDHMS/Routing/module_HYDRO_io.F | 6 +++++- trunk/NDHMS/Routing/module_RT.F | 2 ++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index e392a7f81..e0b22fe0b 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -7282,6 +7282,7 @@ subroutine READ_CHROUTING1( & Tw_CC, n_CC, HRZAREA, LAKEMAXH, & WEIRH, WEIRC, WEIRL, DAML, & ORIFICEC, ORIFICEA, ORIFICEE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & reservoir_type_specified, reservoir_type, & reservoir_parameter_file, LATLAKE, LONLAKE, & ELEVLAKE, dist, ZELEV, LAKENODE, & @@ -7354,7 +7355,7 @@ subroutine READ_CHROUTING1( & integer, intent(INOUT), dimension(:) :: reservoir_type character(len=*), intent(in) :: reservoir_parameter_file real, intent(INOUT), dimension(:) :: LATLAKE,LONLAKE,ELEVLAKE -real, allocatable, dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB +real, intent(INOUT), dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB real, intent(INOUT), dimension(:) :: ChSSlp, Bw, Tw real, intent(INOUT), dimension(:) :: Tw_CC, n_CC ! channel properties of compund @@ -8551,6 +8552,7 @@ subroutine MPP_READ_CHROUTING_new(& n_CC, HRZAREA, LAKEMAXH, & WEIRH, WEIRC, WEIRL, DAML, & ORIFICEC, ORIFICEA, ORIFICEE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & reservoir_type_specified, reservoir_type, & reservoir_parameter_file, LATLAKE, LONLAKE, & ELEVLAKE, dist, ZELEV, LAKENODE, & @@ -8600,6 +8602,7 @@ subroutine MPP_READ_CHROUTING_new(& real, intent(INOUT), dimension(NLAKES) :: LATLAKE,LONLAKE,ELEVLAKE real, intent(INOUT), dimension(NLINKS) :: ChSSlp, Bw, Tw real, intent(INOUT), dimension(NLINKS) :: Tw_CC, n_CC +real, intent(INOUT), dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB character(len=* ) :: geo_finegrid_flnm, route_lake_f character(len=256) :: var_name @@ -8627,6 +8630,7 @@ subroutine MPP_READ_CHROUTING_new(& n_CC, HRZAREA, LAKEMAXH, & WEIRH, WEIRC, WEIRL, DAML, & ORIFICEC, ORIFICEA, ORIFICEE, & + LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB, & reservoir_type_specified, reservoir_type, & reservoir_parameter_file, LATLAKE, LONLAKE, & ELEVLAKE, dist, ZELEV, LAKENODE,& diff --git a/trunk/NDHMS/Routing/module_RT.F b/trunk/NDHMS/Routing/module_RT.F index caa9848f8..0d453c5f8 100644 --- a/trunk/NDHMS/Routing/module_RT.F +++ b/trunk/NDHMS/Routing/module_RT.F @@ -743,6 +743,8 @@ subroutine LandRT_ini(did) rt_domain(did)%WEIRL, rt_domain(did)%DAML, & rt_domain(did)%ORIFICEC, rt_domain(did)%ORIFICEA, & rt_domain(did)%ORIFICEE, & + rt_domain(did)%LAKEMAXDEPTH, rt_domain(did)%LAKEVOL, & + rt_domain(did)%LAKEPOLYA, rt_domain(did)%LAKEPOLYB, & nlst(did)%reservoir_type_specified, rt_domain(did)%reservoir_type, & nlst(did)%reservoir_parameter_file, & rt_domain(did)%LATLAKE, rt_domain(did)%LONLAKE, & From 79260edf1448920df03c16809cb51aa606a84f6f Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Tue, 18 Aug 2020 14:58:59 -0600 Subject: [PATCH 09/12] removed print statments used for testing lake shape option --- trunk/NDHMS/OrchestratorLayer/config.f90 | 1 - trunk/NDHMS/Routing/module_HYDRO_io.F | 3 --- 2 files changed, 4 deletions(-) diff --git a/trunk/NDHMS/OrchestratorLayer/config.f90 b/trunk/NDHMS/OrchestratorLayer/config.f90 index 1a5ba7e7e..61bd22172 100644 --- a/trunk/NDHMS/OrchestratorLayer/config.f90 +++ b/trunk/NDHMS/OrchestratorLayer/config.f90 @@ -742,7 +742,6 @@ subroutine init_namelist_rt_field(did) call hydro_stop("module_namelist: DT not a multiple of DTRT_CH") endif - print *, "Lake Shape Option is:", nlst(did)%lake_shape_option nlst(did)%SUBRTSWCRT = SUBRTSWCRT nlst(did)%OVRTSWCRT = OVRTSWCRT diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index e0b22fe0b..8cb8a6f79 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -8981,14 +8981,11 @@ subroutine read_route_lake_netcdf(route_lake_file, & ! read variables for lake_shape_option not equal to 0 (default Box shape) if(nlst(did)%lake_shape_option /= 0) then - print*, "Different Lake Shape Option (rather than Default Box Shape) is defined in hydro.namelist." - print*, "Variables are all required in LAKEPARM file: LkVol, LkMxDepth, poly_a, poly_b!" call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'poly_b', LAKEPOLYB, 'read_route_lake_netcdf', .TRUE.) else - print*, "Default Lake Shape Option (Box Shape) is defined in hydro.namelist" LAKEVOL = 0.0 LAKEMAXDEPTH = 0.0 LAKEPOLYA = 0.0 From 2de7efff40c614be51e9d98e6397c965557b1c6a Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Wed, 16 Sep 2020 10:15:05 -0600 Subject: [PATCH 10/12] fixed typo in lake shape options rectangular --> triangular --- .../NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F index e3f13a84d..9865d1dc3 100644 --- a/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/trunk/NDHMS/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -260,7 +260,7 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa ! -1 ---> Simple pass through scheme ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) ! 1 ---> CONIC shape with depth based on full lake depth - ! 2 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 2 ---> TRIANGULAR PRISM shape with depth based on full lake depth ! 3 ---> H-A-V RELATIONSHIP based on full lake depth !! ---------------------------- local variables @@ -477,7 +477,7 @@ subroutine STAGE_STORAGE(Deltah,Dmax, h1,Inflow_dt,res_discharge,surf_area,total integer, intent(IN) :: LAKE_SHAPE_OPT ! LAKE_SHAPE_OPT (lake_shape_option) options from namelist: ! 0 ---> BOX shape with depth based on active reservoir depth (NWM v1.x - v2.1) ! 1 ---> CONIC shape with depth based on full lake depth - ! 2 ---> RECTANGLUAR PRISM shape with depth based on full lake depth + ! 2 ---> TRIANGULAR PRISM shape with depth based on full lake depth ! 3 ---> H-A-V RELATIONSHIP based on full lake depth !! ---------------------------- local variables @@ -503,7 +503,7 @@ subroutine STAGE_STORAGE(Deltah,Dmax, h1,Inflow_dt,res_discharge,surf_area,total S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) h2 = (3.*S2/pi*(Dmax/R_cone)**2.)**(1./3.) ! elevation after change in storage in (m) - elseif (LAKE_SHAPE_OPT == 2) then ! RECTANGLUAR PRISM shape + elseif (LAKE_SHAPE_OPT == 2) then ! TRIANGULAR PRISM shape S1 = (surf_area/2./Dmax)*(h1**2.) ! storage before change in volume (m3) S2 = MAX(S1 + dS, 0.0) ! storage after change in volume (m3) h2 = (2.*S2*Dmax/surf_area)**0.5 ! elevation after change in storage in (m) From aef24cb650937c1e43b37317135c67a4be46642f Mon Sep 17 00:00:00 2001 From: bkhazaei <60622555+bkhazaei@users.noreply.github.com> Date: Thu, 17 Jun 2021 08:55:40 -0600 Subject: [PATCH 11/12] Update module_HYDRO_io.F Changed if(nlst(did)%lake_shape_option /= 0) to if(nlst(did)%lake_shape_option >= 0) in Routing/module_HYDRO_io.F to make pass-through scheme work if new LAKEPARM variables are not provided (backward compatible with LAKEPARM v2.1 and before). --- trunk/NDHMS/Routing/module_HYDRO_io.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index 8cb8a6f79..be94f5066 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -8980,7 +8980,7 @@ subroutine read_route_lake_netcdf(route_lake_file, & call get_1d_netcdf_real(ncid, 'lon', lakelon, 'read_route_lake_netcdf', .TRUE.) ! read variables for lake_shape_option not equal to 0 (default Box shape) - if(nlst(did)%lake_shape_option /= 0) then + if(nlst(did)%lake_shape_option >= 0) then call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.) From bc679791ed43ee80b61bef58fb289bf374998fe1 Mon Sep 17 00:00:00 2001 From: bkhazaei <60622555+bkhazaei@users.noreply.github.com> Date: Thu, 17 Jun 2021 14:38:18 -0600 Subject: [PATCH 12/12] Update module_HYDRO_io.F Changed >=0 to > 0 --- trunk/NDHMS/Routing/module_HYDRO_io.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index be94f5066..42583f0d3 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -8980,7 +8980,7 @@ subroutine read_route_lake_netcdf(route_lake_file, & call get_1d_netcdf_real(ncid, 'lon', lakelon, 'read_route_lake_netcdf', .TRUE.) ! read variables for lake_shape_option not equal to 0 (default Box shape) - if(nlst(did)%lake_shape_option >= 0) then + if(nlst(did)%lake_shape_option > 0) then call get_1d_netcdf_real(ncid, 'LkVol', LAKEVOL, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'LkMxDepth', LAKEMAXDEPTH, 'read_route_lake_netcdf', .TRUE.) call get_1d_netcdf_real(ncid, 'poly_a', LAKEPOLYA, 'read_route_lake_netcdf', .TRUE.)