diff --git a/trunk/NDHMS/Data_Rec/module_namelist.F b/trunk/NDHMS/Data_Rec/module_namelist.F index 758db5c85..c97daefbd 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. 3) ) 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..61bd22172 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. 3) ) 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 & @@ -463,6 +466,7 @@ subroutine init_namelist_rt_field(did) sys_cpl, rst_typ, rst_bi_in, rst_bi_out, & gwChanCondSw, GwPreCycles, GwSpinCycles, GwPreDiagInterval, gwsoilcpl, & 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="" @@ -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, & @@ -587,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 @@ -713,6 +719,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 +742,7 @@ subroutine init_namelist_rt_field(did) call hydro_stop("module_namelist: DT not a multiple of DTRT_CH") endif + 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..9865d1dc3 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,49 +247,68 @@ 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 ---> CONIC 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 - 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 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 + ! 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 + + ! 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) 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,8 +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 - 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 @@ -329,7 +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 + 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 @@ -356,7 +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 + 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 @@ -403,9 +446,7 @@ 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.... return @@ -413,4 +454,72 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa 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) :: 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 ---> CONIC 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 + 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 :: 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) then ! BOX shape + ! h2 = (dS / surf_area) + h1 + + 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) + 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 ! 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) + + 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 + + ! ---------------------------------------------------------------- + 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..42583f0d3 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,6 +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, 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 @@ -7368,6 +7370,7 @@ subroutine READ_CHROUTING1( & integer :: did logical :: fexist + did = 1 !--------------------------------------------------------- @@ -7470,7 +7473,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 +7513,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 @@ -8175,6 +8184,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, & @@ -8194,6 +8204,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 @@ -8221,7 +8232,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 @@ -8246,6 +8258,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) @@ -8263,7 +8279,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 @@ -8285,6 +8302,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, intent(out), dimension(:) :: LAKEMAXDEPTH, LAKEVOL, LAKEPOLYA, LAKEPOLYB !end NLAKES INTEGER, dimension(GNLINKSL) :: tmpLAKEIDA, tmpLINKID, tmpTO_NODE, tmpTYPEL, tmpORDER @@ -8370,7 +8388,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 @@ -8427,6 +8447,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) @@ -8528,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, & @@ -8577,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 @@ -8604,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,& @@ -8907,7 +8934,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, LAKEPOLYA, LAKEPOLYB, & + ELEVLAKE, NLAKES) implicit none character(len=*), intent(in) :: route_lake_file @@ -8917,6 +8946,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, LAKEPOLYA, LAKEPOLYB real, dimension(:), intent(out) :: ELEVLAKE integer, dimension(:), intent(out) :: reservoir_type @@ -8948,6 +8978,20 @@ 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.) + + ! read variables for lake_shape_option not equal to 0 (default Box shape) + 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.) + call get_1d_netcdf_real(ncid, 'poly_b', LAKEPOLYB, 'read_route_lake_netcdf', .TRUE.) + else + 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.) @@ -8969,6 +9013,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 e71627b42..0d453c5f8 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 @@ -734,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, & @@ -778,6 +789,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, & @@ -834,6 +847,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 +1673,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