SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
waterdist_module Module Reference

Functions/Subroutines

subroutine drainage (is, state_is, StorCap, DrainEq, DrainCoef1, DrainCoef2, nsh_real, drain_is)
 
subroutine cal_water_storage (is, sfr_surf, PipeCapacity, RunoffToWater, pin, WU_surf, drain_surf, AddWater, addImpervious, nsh_real, state_in, frac_water2runoff, PervFraction, addVeg, SoilStoreCap, addWaterBody, FlowChange, StateLimit, runoffAGimpervious, runoffAGveg, runoffPipes, ev, soilstore_id, surplusWaterBody, SurplusEvap, runoffWaterBody, runoff, state_out)
 
subroutine cal_water_storage_surf (pin, nsh_real, PipeCapacity, RunoffToWater, addImpervious, addVeg, addWaterBody, FlowChange, SoilStoreCap_surf, StateLimit_surf, PervFraction, sfr_surf, drain_surf, AddWater_surf, frac_water2runoff_surf, WU_surf, ev_surf_in, state_surf_in, soilstore_surf_in, ev_surf_out, state_surf_out, soilstore_surf_out, runoff_surf, runoffAGimpervious_grid, runoffAGveg_grid, runoffPipes_grid, runoffWaterBody_grid)
 
subroutine cal_water_storage_building (pin, nsh_real, nlayer, sfr_roof, StateLimit_roof, SoilStoreCap_roof, WetThresh_roof, ev_roof_in, state_roof_in, soilstore_roof_in, sfr_wall, StateLimit_wall, SoilStoreCap_wall, WetThresh_wall, ev_wall_in, state_wall_in, soilstore_wall_in, ev_roof_out, state_roof_out, soilstore_roof_out, runoff_roof, ev_wall_out, state_wall_out, soilstore_wall_out, runoff_wall, state_building, soilstore_building, runoff_building, SoilStoreCap_building)
 
subroutine updateflood (is, runoff, sfr_surf, PipeCapacity, RunoffToWater, runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes)
 
subroutine redistributewater (SnowUse, WaterDist, sfr_surf, Drain, AddWaterRunoff, AddWater)
 
subroutine suews_update_soilmoist (NonWaterFraction, SoilStoreCap, sfr_surf, soilstore_id, SoilMoistCap, SoilState, vsmd, smd)
 
subroutine suews_cal_soilstate (SMDMethod, xsmd, NonWaterFraction, SoilMoistCap, SoilStoreCap, surf_chang_per_tstep, soilstore_id, soilstoreOld, sfr_surf, smd, smd_nsurf, tot_chang_per_tstep, SoilState)
 
subroutine suews_cal_horizontalsoilwater (sfr_surf, SoilStoreCap, SoilDepth, SatHydraulicConduct, SurfaceArea, NonWaterFraction, tstep_real, soilstore_id, runoffSoil, runoffSoil_per_tstep)
 
subroutine suews_cal_wateruse (nsh_real, wu_m3, SurfaceArea, sfr_surf, IrrFracPaved, IrrFracBldgs, IrrFracEveTr, IrrFracDecTr, IrrFracGrass, IrrFracBSoil, IrrFracWater, DayofWeek_id, WUProfA_24hr, WUProfM_24hr, InternalWaterUse_h, HDD_id, WUDay_id, WaterUseMethod, NSH, it, imin, DLS, wu_surf, wu_int, wu_ext)
 

Function/Subroutine Documentation

◆ cal_water_storage()

subroutine waterdist_module::cal_water_storage ( integer, intent(in)  is,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(in)  PipeCapacity,
real(kind(1d0)), intent(in)  RunoffToWater,
real(kind(1d0)), intent(in)  pin,
real(kind(1d0)), dimension(nsurf), intent(in)  WU_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  drain_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  AddWater,
real(kind(1d0)), intent(in)  addImpervious,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), dimension(nsurf), intent(in)  state_in,
real(kind(1d0)), dimension(nsurf), intent(in)  frac_water2runoff,
real(kind(1d0)), intent(in)  PervFraction,
real(kind(1d0)), intent(in)  addVeg,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap,
real(kind(1d0)), intent(in)  addWaterBody,
real(kind(1d0)), intent(in)  FlowChange,
real(kind(1d0)), dimension(nsurf), intent(in)  StateLimit,
real(kind(1d0)), intent(inout)  runoffAGimpervious,
real(kind(1d0)), intent(inout)  runoffAGveg,
real(kind(1d0)), intent(inout)  runoffPipes,
real(kind(1d0)), intent(inout)  ev,
real(kind(1d0)), dimension(nsurf), intent(inout)  soilstore_id,
real(kind(1d0)), intent(inout)  surplusWaterBody,
real(kind(1d0)), dimension(2), intent(inout)  SurplusEvap,
real(kind(1d0)), intent(inout)  runoffWaterBody,
real(kind(1d0)), dimension(nsurf), intent(out)  runoff,
real(kind(1d0)), dimension(nsurf), intent(out)  state_out 
)

Definition at line 92 of file suews_phys_waterdist.f95.

100 !------------------------------------------------------------------------------
101 !Calculation of storage change
102 ! TS 30 Nov 2019
103 ! - Allow irrigation on all surfaces (previously only on vegetated surfaces)
104 ! LJ 27 Jan 2016
105 ! - Removed tabs and cleaned the code
106 ! HCW 08 Dec 2015
107 ! -Added if-loop check for no Paved surfaces
108 ! LJ 6 May 2015
109 ! - Calculations of the piperunoff exceedings moved to separate subroutine updateFlood.
110 ! - Now also called from snow subroutine
111 ! - Evaporation is modified using EvapPart
112 ! - when no water on impervious surfaces, evap occurs above pervious surfaces instead
113 ! Rewritten by HCW 12 Feb 2015
114 ! - Old variable 'p' for water input to the surface renamed to 'p_mm'
115 ! - All water now added to p_mm first, before threshold checks or other calculations
116 ! - Water from other grids now added to p_mm (instead of state_id for impervious surfaces)
117 ! - Removed division of runoff by nsh, as whole model now runs at the same timestep
118 ! - Adjusted transfer of ev between surfaces to conserve mass (not depth)
119 ! - Volumes used for water transport between grids to account for SurfaceArea changing between grids
120 ! - Added threshold check for state_id(WaterSurf) - was going negative
121 ! Last modified HCW 09 Feb 2015
122 ! - Removed StorCap input because it is provided by module allocateArray
123 ! - Tidied and commented code
124 ! Modified by LJ in November 2012:
125 ! - P>10 was not taken into account for impervious surfaces - Was fixed.
126 ! - Above impervious surfaces possibility of the state_id to exceed max capacity was limited
127 ! although this should be possible - was fixed
128 ! Modified by LJ 10/2010
129 ! Rewritten mostly by LJ in 2010
130 ! To do:
131 ! - Finish area normalisation for RG2G & finish coding GridConnections
132 ! - What is the 10 mm hr-1 threshold for?
133 ! - Decide upon and correct storage capacities here & in evap subroutine
134 ! - FlowChange units should be mm hr-1 - need to update everywhere
135 ! - Add SurfaceFlood(is)?
136 ! - What happens if sfr_surf(is) = 0 or 1?
137 ! - Consider how irrigated trees actually works...
138 !------------------------------------------------------------------------------
139
140 IMPLICIT NONE
141
142 !Stores flood water when surface state_id exceeds storage capacity [mm]
143 INTEGER, INTENT(in) :: is ! surface type
144
145 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! surface fractions
146 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: AddWater !Water from other surfaces (WGWaterDist in SUEWS_ReDistributeWater.f95) [mm]
147 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: state_in !Wetness status of each surface type from previous timestep [mm]
148 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: frac_water2runoff !Fraction of water going to runoff/sub-surface soil (WGWaterDist) [-]
149 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap !Capacity of soil store for each surface [mm]
150 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: StateLimit !Limit for state_id of each surface type [mm] (specified in input files)
151 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: drain_surf !Drainage of each surface type [mm]
152 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: WU_surf !external water use of each surface type [mm]
153
154 REAL(KIND(1D0)), INTENT(in) :: PipeCapacity !Capacity of pipes to transfer water
155 REAL(KIND(1D0)), INTENT(in) :: RunoffToWater !Fraction of surface runoff going to water body
156 REAL(KIND(1D0)), INTENT(in) :: pin !Rain per time interval
157 REAL(KIND(1D0)), INTENT(in) :: addImpervious !Water from impervious surfaces of other grids [mm] for whole surface area
158 REAL(KIND(1D0)), INTENT(in) :: nsh_real !nsh cast as a real for use in calculations
159 REAL(KIND(1D0)), INTENT(in) :: PervFraction ! sum of surface cover fractions for impervious surfaces
160 REAL(KIND(1D0)), INTENT(in) :: addVeg !Water from vegetated surfaces of other grids [mm] for whole surface area
161 REAL(KIND(1D0)), INTENT(in) :: addWaterBody !Water from water surface of other grids [mm] for whole surface area
162 REAL(KIND(1D0)), INTENT(in) :: FlowChange !Difference between the input and output flow in the water body
163
164 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: soilstore_id !Soil moisture of each surface type [mm]
165 REAL(KIND(1D0)), DIMENSION(2), INTENT(inout) :: SurplusEvap !Surplus for evaporation in 5 min timestep
166
167 REAL(KIND(1D0)), DIMENSION(nsurf) :: chang !Change in state_id [mm]
168 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: runoff !Runoff from each surface type [mm]
169 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: state_out !Wetness status of each surface type [mm]
170
171 ! =============================
172 ! TS 01 Apr 2022:
173 ! these variables were used as inout variables in the original code for water transfer between grids
174 ! but they are not used in the new code, so they are removed here
175 REAL(KIND(1D0)), INTENT(inout) :: surplusWaterBody !Extra runoff that goes to water body [mm] as specified by RunoffToWater
176 REAL(KIND(1D0)), INTENT(inout) :: runoffAGimpervious !Above ground runoff from impervious surface [mm] for whole surface area
177 REAL(KIND(1D0)), INTENT(inout) :: runoffAGveg !Above ground runoff from vegetated surfaces [mm] for whole surface area
178 REAL(KIND(1D0)), INTENT(inout) :: runoffPipes !Runoff in pipes [mm] for whole surface area
179 REAL(KIND(1D0)), INTENT(inout) :: ev !Evaporation
180 REAL(KIND(1D0)), INTENT(inout) :: runoffWaterBody !Above ground runoff from water surface [mm] for whole surface area
181
182 REAL(KIND(1D0)) :: p_mm !Inputs to surface water balance
183
184 !Extra evaporation [mm] from impervious surfaces which cannot happen due to lack of water
185 REAL(KIND(1D0)) :: EvPart
186 REAL(KIND(1D0)), PARAMETER :: NotUsed = -55.5
187
188 !Threshold for intense precipitation [mm hr-1]
189 REAL(KIND(1D0)), PARAMETER :: IPThreshold_mmhr = 10 ! NB:this should be an input and can be specified. SG 25 Apr 2018
190
191 !Initialise extra evaporation to zero
192 evpart = 0
193
194 !Initialise runoff to zero
195 runoff(is) = 0
196
197 !SurfaceFlood(is) = 0 !!This probably needs to be carried over between timesteps, but reset for now
198
199 !==================================================================
200 ! Combine water inputs to the current surface
201 ! Add external water use for each surface type
202 p_mm = pin + wu_surf(is)
203
204 ! Add water from other surfaces within the same grid (RS2S) ----
205 ! AddWater is the water supplied to the current surface from other surfaces
206 ! i.e. drain*WaterDist (see SUEWS_ReDistributeWater)
207 p_mm = p_mm + addwater(is)
208 !==================================================================
209
210 !========surface-specific calculation=========================
211 SELECT CASE (is)
212 CASE (pavsurf, bldgsurf)
213 !==== Impervious surfaces (Paved, Buildings) ======================
214 ! Add water from neighbouring grids (RG2G)
215 ! Add to PavSurf only, as water cannot flow onto buildings
216 IF (is == pavsurf) THEN
217 IF (sfr_surf(pavsurf) /= 0) THEN ! If loop added HCW 08 Dec 2015
218 p_mm = p_mm + addimpervious/sfr_surf(pavsurf)
219 END IF
220 END IF
221
222 ! Calculate change in surface state_id (inputs - outputs)
223 chang(is) = p_mm - (drain_surf(is) + ev)
224
225 ! If p_mm is too large, excess goes to runoff (i.e. the rate of water supply is too fast)
226 ! and does not affect state_id
227 IF (p_mm > ipthreshold_mmhr/nsh_real) THEN
228 runoff(is) = runoff(is) + (p_mm - ipthreshold_mmhr/nsh_real)
229 chang(is) = ipthreshold_mmhr/nsh_real - (drain_surf(is) + ev)
230 END IF
231
232 ! Calculate updated state_id using chang
233 state_out(is) = state_in(is) + chang(is)
234
235 ! Check state_id is within physical limits between zero (dry) and max. storage capacity
236 IF (state_out(is) < 0.0) THEN ! Cannot have a negative surface state_id
237 ! If there is not sufficient water on the surface, then don't allow this evaporation to happen
238 ! Allow evaporation only until surface is dry (state_id(is)=0); additional evaporation -> evaporation surplus
239 surplusevap(is) = abs(state_out(is)) !Surplus evaporation is that which tries to evaporate non-existent water
240 ev = ev - surplusevap(is) !Limit evaporation according to water availability
241 state_out(is) = 0.0 !Now surface is dry
242 ! elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
243 ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
244 ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
245 ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
246 ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
247 END IF
248
249 ! Recalculate change in surface state_id from difference with previous timestep
250 chang(is) = state_out(is) - state_in(is)
251
252 ! Runoff -------------------------------------------------------
253 ! For impervious surfaces, some of drain(is) becomes runoff
254 runoff(is) = runoff(is) + drain_surf(is)*frac_water2runoff(is) !Drainage (that is not flowing to other surfaces) goes to runoff
255
256 !So, up to this point, runoff(is) can have contributions if
257 ! p_mm > ipthreshold (water input too fast)
258 ! state_id > StoreDrainPrm(6,is) (net water exceeds storage capacity)
259 ! WaterDist specifies some fraction of drain(is) -> runoff
260
261 CASE (conifsurf:bsoilsurf)
262 !==== For Conif, Decid, Grass, BSoil surfaces ==================
263 ! Transfer evaporation surplus from impervious surfaces to pervious surfaces
264 evpart = merge( &
265 dot_product(surplusevap(pavsurf:bldgsurf), sfr_surf(pavsurf:bldgsurf)/pervfraction), &
266 0d0, &
267 pervfraction /= 0)
268
269 ! Add surplus evaporation to ev for pervious surfaces
270 ev = ev + evpart
271
272 ! ---- Add water from neighbouring grids (RG2G) ----
273 ! Add to Grass and BSoil only, as water cannot flow onto trees
274 IF (is == grasssurf .OR. is == bsoilsurf) THEN
275 IF ((sfr_surf(grasssurf) + sfr_surf(bsoilsurf)) /= 0) THEN
276 p_mm = p_mm + addveg/(sfr_surf(grasssurf) + sfr_surf(bsoilsurf))
277 END IF
278 END IF
279
280 ! Calculate change in surface state_id (inputs - outputs)
281 chang(is) = p_mm - (drain_surf(is) + ev)
282
283 ! If p_mm is too large, excess goes to runoff (i.e. the rate of water supply is too fast)
284 ! and does not affect state_id
285 IF (p_mm > ipthreshold_mmhr/nsh_real) THEN
286 runoff(is) = runoff(is) + (p_mm - ipthreshold_mmhr/nsh_real)
287 chang(is) = ipthreshold_mmhr/nsh_real - (drain_surf(is) + ev)
288 END IF
289
290 ! Calculate updated state_id using chang
291 state_out(is) = state_in(is) + chang(is)
292
293 ! Check state_id is within physical limits between zero (dry) and max. storage capacity
294 IF (state_out(is) < 0.0) THEN ! Cannot have a negative surface state_id
295 ! If there is not sufficient water on the surface, then remove water from soilstore
296 ! Allow evaporation until soilstore_id is depleted and surface is dry
297 IF ((soilstore_id(is) + state_out(is)) >= 0) THEN
298 soilstore_id(is) = soilstore_id(is) + state_out(is)
299 state_out(is) = 0.0
300 ! If there is not sufficient water on the surface or soilstore, then don't allow this evaporation to happen
301 ELSE
302 ev = ev - abs(state_out(is)) !Limit evaporation according to water availability
303 state_out(is) = 0.0 !Now surface is dry
304 END IF
305
306 !elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
307 ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
308 ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
309 ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
310 ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
311 END IF
312
313 ! Recalculate change in surface state_id from difference with previous timestep
314 chang(is) = state_out(is) - state_in(is)
315
316 !Where should this go? Used to be before previous part!!
317 ! soilstore_id -------------------------------------------------
318 ! For pervious surfaces (not water), some of drain(is) goes to soil storage
319 ! Drainage (that is not flowing to other surfaces) goes to soil storages
320 soilstore_id(is) = soilstore_id(is) + drain_surf(is)*frac_water2runoff(is)
321
322 ! If soilstore is full, the excess will go to runoff
323 IF (soilstore_id(is) > soilstorecap(is)) THEN ! TODO: this should also go to flooding of some sort
324 runoff(is) = runoff(is) + (soilstore_id(is) - soilstorecap(is))
325 soilstore_id(is) = soilstorecap(is)
326 ELSEIF (soilstore_id(is) < 0) THEN !! QUESTION: But where does this lack of water go? !!Can this really happen here?
327 CALL errorhint(62, 'SUEWS_store: soilstore_id(is) < 0 ', soilstore_id(is), notused, is)
328 ! Code this properly - soilstore_id(is) < 0 shouldn't happen given the above loops
329 !soilstore_id(is)=0 !Groundwater / deeper soil should kick in
330 END IF
331
332 CASE (watersurf)
333 IF (sfr_surf(watersurf) /= 0) THEN
334
335 ! ---- Add water from neighbouring grids (RG2G) ----
336 p_mm = p_mm + addwaterbody/sfr_surf(watersurf)
337
338 ! Calculate change in surface state_id (inputs - outputs)
339 ! No drainage for water surface
340 ! FlowChange is the difference in input and output flows [mm hr-1]
341 chang(is) = p_mm + flowchange/nsh_real - ev
342
343 ! Calculate updated state_id using chang
344 state_out(is) = state_in(is) + chang(is)
345
346 ! Check state_id is within physical limits between zero (dry) and max. storage capacity
347 IF (state_out(is) < 0.0) THEN ! Cannot have a negative surface state_id
348 ! If there is not sufficient water on the surface, then don't allow this evaporation to happen
349 ev = ev - abs(state_out(is)) !Limit evaporation according to water availability
350 state_out(is) = 0.0 !Now surface is dry
351 !elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
352 ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
353 ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
354 ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
355 ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
356 END IF
357
358 ! Recalculate change in surface state_id from difference with previous timestep
359 chang(is) = state_out(is) - state_in(is)
360
361 ! If state_id exceeds limit, then excess goes to runoff (currently applies to water StoreDrainPrm only)
362 IF (state_out(watersurf) > statelimit(watersurf)) THEN
363 runoff(watersurf) = runoff(watersurf) + (state_out(watersurf) - statelimit(watersurf))
364 state_out(watersurf) = statelimit(watersurf)
365 runoffwaterbody = runoffwaterbody + runoff(watersurf)*sfr_surf(watersurf)
366 ELSE
367 state_out(watersurf) = state_out(watersurf) + surpluswaterbody
368 IF (state_out(watersurf) > statelimit(watersurf)) THEN
369 runoffwaterbody = runoffwaterbody + (state_out(watersurf) - statelimit(watersurf))*sfr_surf(watersurf)
370 state_out(watersurf) = statelimit(watersurf)
371 END IF
372 END IF
373
374 ! Recalculate change in surface state_id from difference with previous timestep
375 chang(is) = state_out(is) - state_in(is)
376 END IF
377 END SELECT
378 !==================================================================
379
380 !==== RUNOFF ======================================================
381 ! TODO: to consider areas here - SurfaceArea may vary between grids too
382 ! - also implement where water for next surface is calculated (RunoffFromGrid subroutine)
383 ! Calculations of the piperunoff exceedensances moved to separate subroutine so that from snow same
384 ! calculations can be made. LJ in May 2015
385
386 IF (is < watersurf) THEN !Not for water body
387 ! CALL updateFlood
388 CALL updateflood( &
389 is, runoff, & ! input:
390 sfr_surf, pipecapacity, runofftowater, &
391 runoffagimpervious, surpluswaterbody, runoffagveg, runoffpipes) ! inout:
392 END IF
393
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)

References allocatearray::bldgsurf, allocatearray::bsoilsurf, allocatearray::conifsurf, errorhint(), allocatearray::grasssurf, allocatearray::pavsurf, updateflood(), and allocatearray::watersurf.

Referenced by cal_water_storage_surf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_water_storage_building()

subroutine waterdist_module::cal_water_storage_building ( real(kind(1d0)), intent(in)  pin,
real(kind(1d0)), intent(in)  nsh_real,
integer, intent(in)  nlayer,
real(kind(1d0)), dimension(nlayer), intent(in)  sfr_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  StateLimit_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  SoilStoreCap_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  WetThresh_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  ev_roof_in,
real(kind(1d0)), dimension(nlayer), intent(in)  state_roof_in,
real(kind(1d0)), dimension(nlayer), intent(in)  soilstore_roof_in,
real(kind(1d0)), dimension(nlayer), intent(in)  sfr_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  StateLimit_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  SoilStoreCap_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  WetThresh_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  ev_wall_in,
real(kind(1d0)), dimension(nlayer), intent(in)  state_wall_in,
real(kind(1d0)), dimension(nlayer), intent(in)  soilstore_wall_in,
real(kind(1d0)), dimension(nlayer), intent(out)  ev_roof_out,
real(kind(1d0)), dimension(nlayer), intent(out)  state_roof_out,
real(kind(1d0)), dimension(nlayer), intent(out)  soilstore_roof_out,
real(kind(1d0)), dimension(nlayer), intent(out)  runoff_roof,
real(kind(1d0)), dimension(nlayer), intent(out)  ev_wall_out,
real(kind(1d0)), dimension(nlayer), intent(out)  state_wall_out,
real(kind(1d0)), dimension(nlayer), intent(out)  soilstore_wall_out,
real(kind(1d0)), dimension(nlayer), intent(out)  runoff_wall,
real(kind(1d0)), intent(out)  state_building,
real(kind(1d0)), intent(out)  soilstore_building,
real(kind(1d0)), intent(out)  runoff_building,
real(kind(1d0)), intent(out)  SoilStoreCap_building 
)

Definition at line 516 of file suews_phys_waterdist.f95.

525
526 IMPLICIT NONE
527
528 INTEGER, INTENT(in) :: nlayer !number of layers
529 REAL(KIND(1D0)), INTENT(in) :: pin !Rain per time interval
530 REAL(KIND(1D0)), INTENT(in) :: nsh_real !timesteps per hour
531
532 ! input for generic roof facets
533 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_roof
534 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: StateLimit_roof !Limit for state_id of each surface type [mm] (specified in input files)
535 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: WetThresh_roof
536 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: SoilStoreCap_roof
537 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: state_roof_in
538 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: soilstore_roof_in
539 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: ev_roof_in
540
541 ! input for generic wall facets
542 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_wall
543 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: StateLimit_wall
544 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: WetThresh_wall
545 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: SoilStoreCap_wall
546 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: state_wall_in
547 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: soilstore_wall_in
548 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: ev_wall_in
549
550 ! output for generic roof facets
551 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: ev_roof_out
552 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: state_roof_out
553 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: soilstore_roof_out
554 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: runoff_roof
555
556 ! output for generic wall facets
557 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: ev_wall_out
558 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: state_wall_out
559 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: soilstore_wall_out
560 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: runoff_wall
561
562 REAL(KIND(1D0)), INTENT(out) :: state_building !aggregated surface water of building facets [mm]
563 REAL(KIND(1D0)), INTENT(out) :: soilstore_building ! aggregated soilstore of building facets[mm]
564 REAL(KIND(1D0)), INTENT(out) :: runoff_building !aggregated Runoff of building facets [mm]
565 REAL(KIND(1D0)), INTENT(out) :: SoilStoreCap_building
566
567 REAL(KIND(1D0)) :: precip_excess_roof !precipitation excess above IPThreshold_mmhr [mm]
568 REAL(KIND(1D0)) :: precip_excess_wall !precipitation excess above limit [mm]
569 REAL(KIND(1D0)) :: pin_wall !input water to wall facet [mm]
570
571 REAL(KIND(1D0)), DIMENSION(nlayer) :: chang_roof !Runoff from each surface type [mm]
572 REAL(KIND(1D0)), DIMENSION(nlayer) :: chang_wall !Runoff from each surface type [mm]
573 REAL(KIND(1D0)), DIMENSION(nlayer) :: drain_roof !drainage to sewer [mm]
574 REAL(KIND(1D0)), DIMENSION(nlayer) :: drain_wall !drainage to sewer [mm]
575 REAL(KIND(1D0)), DIMENSION(nlayer) :: infil_roof !infiltration to replenish soil water [mm]
576 REAL(KIND(1D0)), DIMENSION(nlayer) :: infil_wall !infiltration to replenish soil water [mm]
577 REAL(KIND(1D0)), DIMENSION(nlayer) :: evap_roof !evapotranpiration from each surface type [mm]
578 REAL(KIND(1D0)), DIMENSION(nlayer) :: evap_wall !evapotranpiration from each surface type [mm]
579 ! REAL(KIND(1D0)), DIMENSION(nsurf) :: soilstore !Soil moisture of each surface type [mm]
580 ! REAL(KIND(1D0)), DIMENSION(nsurf) :: chang !Change in state_id [mm]
581 ! REAL(KIND(1D0)), DIMENSION(2) :: SurplusEvap !Surplus for evaporation in 5 min timestep
582
583 !Threshold for intense precipitation [mm hr-1]
584 REAL(KIND(1D0)), PARAMETER :: IPThreshold_mmhr = 10 ! NB:this should be an input and can be specified. SG 25 Apr 2018
585 REAL(KIND(1D0)) :: IPThreshold ! NB:this should be an input and can be specified. SG 25 Apr 2018
586
587 INTEGER :: n_layer, i_layer
588
589 ! Threshold for intense precipitation in each timestep
590 ipthreshold = ipthreshold_mmhr/nsh_real
591
592 ! precipitation excess above IPThreshold_mmhr
593 precip_excess_roof = pin - ipthreshold
594
595 ! n_layer = SIZE(sfr_roof)
596 DO i_layer = 1, nlayer
597
598 ! ===== roof part =====
599 ! initialize temporary variables
600 runoff_roof(i_layer) = 0.0
601 chang_roof(i_layer) = 0.0
602
603 ! let's assume so for now
604 ! TODO: introduce a more sophisticated drinage function for green roofs
605 drain_roof(i_layer) = state_roof_in(i_layer)*.0
606
607 ! infiltration to replenish soil water
608 infil_roof(i_layer) = state_roof_in(i_layer)*.3
609
610 IF (precip_excess_roof > 0) THEN
611 ! runoff generated from roof
612 runoff_roof(i_layer) = precip_excess_roof
613 chang_roof(i_layer) = ipthreshold - ev_roof_in(i_layer) - drain_roof(i_layer) - infil_roof(i_layer)
614 ELSE
615 chang_roof(i_layer) = pin - ev_roof_in(i_layer) - drain_roof(i_layer) - infil_roof(i_layer)
616 END IF
617
618 ! change in surface water
619 state_roof_out(i_layer) = state_roof_in(i_layer) + chang_roof(i_layer)
620
621 ! Check state_id is within physical limits between zero (dry) and max. storage capacity
622 IF (state_roof_out(i_layer) < 0.0) THEN ! Cannot have a negative surface state_id
623 ! If there is not sufficient water on the surface, then remove water from soilstore
624 ! Allow evaporation until soilstore_id is depleted and surface is dry
625 IF ((soilstore_roof_in(i_layer) + state_roof_out(i_layer)) >= 0) THEN
626 ! update soilstore with deficit in surface state
627 soilstore_roof_out(i_layer) = soilstore_roof_in(i_layer) + state_roof_out(i_layer)
628 ! If there is not sufficient water on the surface or soilstore, then don't allow this evaporation to happen
629 ELSE
630 ev_roof_out(i_layer) = ev_roof_in(i_layer) - abs(state_roof_out(i_layer)) !Limit evaporation according to water availability
631 END IF
632 ! force surface to dry
633 state_roof_out(i_layer) = 0.0
634 ELSE
635 ! If there is sufficient water on the surface, then allow surface water to replenish soilstore
636 IF (soilstore_roof_in(i_layer) + infil_roof(i_layer) > soilstorecap_roof(i_layer)) THEN
637 ! cap soilstore_id at max. storage capacity if saturated
638 soilstore_roof_out(i_layer) = soilstorecap_roof(i_layer)
639
640 ! excessive infiltration goes into runoff
641 runoff_roof = runoff_roof + (soilstore_roof_in(i_layer) + infil_roof(i_layer) - soilstorecap_roof(i_layer))
642 ELSE
643 soilstore_roof_out(i_layer) = soilstore_roof_in(i_layer) + infil_roof(i_layer)
644 END IF
645
646 END IF
647
648 ! ===== wall part =====
649 ! initialize temporary variables
650 runoff_wall(i_layer) = 0.0
651 chang_wall(i_layer) = 0.0
652 pin_wall = 0.0
653
654 ! runoff from roof goes into wall as water supply
655 ! NB: only a fraction of precipitation is diverted to the wall
656 pin_wall = runoff_roof(i_layer) + pin*.2
657 precip_excess_wall = pin_wall - statelimit_wall(i_layer)
658
659 ! let's assume so for now
660 ! TODO: introduce a more sophisticated drinage function for green roofs
661 drain_wall(i_layer) = state_wall_in(i_layer)*.0
662
663 ! infiltration to replenish soil water
664 infil_wall(i_layer) = state_wall_in(i_layer)*.1
665
666 IF (precip_excess_wall > 0) THEN
667 ! runoff generated from roof
668 runoff_wall(i_layer) = precip_excess_wall
669 chang_wall(i_layer) = statelimit_wall(i_layer) - ev_wall_in(i_layer) - drain_wall(i_layer) - infil_wall(i_layer)
670 ELSE
671 chang_wall(i_layer) = pin_wall - ev_wall_in(i_layer) - drain_wall(i_layer) - infil_wall(i_layer)
672 END IF
673
674 ! change in surface water
675 state_wall_out(i_layer) = state_wall_in(i_layer) + chang_wall(i_layer)
676
677 ! Check state_id is within physical limits between zero (dry) and max. storage capacity
678 IF (state_wall_out(i_layer) < 0.0) THEN ! Cannot have a negative surface state_id
679 ! If there is not sufficient water on the surface, then remove water from soilstore
680 ! Allow evaporation until soilstore_id is depleted and surface is dry
681 IF ((soilstore_wall_in(i_layer) + state_wall_out(i_layer)) >= 0) THEN
682 ! update soilstore with deficit in surface state
683 soilstore_wall_out(i_layer) = soilstore_wall_in(i_layer) + state_wall_out(i_layer)
684 ! If there is not sufficient water on the surface or soilstore, then don't allow this evaporation to happen
685 ELSE
686 ev_wall_out(i_layer) = ev_wall_in(i_layer) - abs(state_wall_out(i_layer)) !Limit evaporation according to water availability
687 END IF
688 ! force surface to dry
689 state_wall_out(i_layer) = 0.0
690 ELSE
691 ! If there is sufficient water on the surface, then allow surface water to replenish soilstore
692 IF (soilstore_wall_in(i_layer) + infil_wall(i_layer) > soilstorecap_wall(i_layer)) THEN
693 ! cap soilstore_id at max. storage capacity if saturated
694 soilstore_wall_out(i_layer) = soilstorecap_wall(i_layer)
695
696 ! excessive infiltration goes into runoff
697 runoff_wall(i_layer) = &
698 runoff_wall(i_layer) + &
699 (soilstore_wall_in(i_layer) + infil_wall(i_layer) - soilstorecap_wall(i_layer))
700 ELSE
701 soilstore_wall_out(i_layer) = soilstore_wall_in(i_layer) + infil_wall(i_layer)
702 END IF
703
704 END IF
705
706 END DO
707
708 ! diagnostics info:
709 ! call r8vec_print(SIZE(soilstore_roof_out), soilstore_roof_in, 'soilstore_roof_in in layer')
710 ! call r8vec_print(SIZE(soilstore_roof_out), SoilStoreCap_roof, 'SoilStoreCap_roof in layer')
711 ! call r8vec_print(SIZE(soilstore_roof_out), infil_roof, 'infil_roof in layer')
712 ! call r8vec_print(SIZE(soilstore_wall_out), soilstore_roof_out, 'soilstore_roof_out in layer')
713
714 ! aggregated values
715 state_building = dot_product(state_roof_out, sfr_roof) + dot_product(state_wall_out, sfr_wall)
716 soilstore_building = dot_product(soilstore_roof_out, sfr_roof) + dot_product(soilstore_wall_out, sfr_wall)
717 soilstorecap_building = dot_product(soilstorecap_roof, sfr_roof) + dot_product(soilstorecap_wall, sfr_wall)
718 ! only allow runoff from walls
719 runoff_building = dot_product(runoff_wall, sfr_wall)
720

Referenced by suews_driver::suews_cal_qe().

Here is the caller graph for this function:

◆ cal_water_storage_surf()

subroutine waterdist_module::cal_water_storage_surf ( real(kind(1d0)), intent(in)  pin,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), intent(in)  PipeCapacity,
real(kind(1d0)), intent(in)  RunoffToWater,
real(kind(1d0)), intent(in)  addImpervious,
real(kind(1d0)), intent(in)  addVeg,
real(kind(1d0)), intent(in)  addWaterBody,
real(kind(1d0)), intent(in)  FlowChange,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  StateLimit_surf,
real(kind(1d0)), intent(in)  PervFraction,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  drain_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  AddWater_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  frac_water2runoff_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  WU_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  ev_surf_in,
real(kind(1d0)), dimension(nsurf), intent(in)  state_surf_in,
real(kind(1d0)), dimension(nsurf), intent(in)  soilstore_surf_in,
real(kind(1d0)), dimension(nsurf), intent(out)  ev_surf_out,
real(kind(1d0)), dimension(nsurf), intent(out)  state_surf_out,
real(kind(1d0)), dimension(nsurf), intent(out)  soilstore_surf_out,
real(kind(1d0)), dimension(nsurf), intent(out)  runoff_surf,
real(kind(1d0)), intent(out)  runoffAGimpervious_grid,
real(kind(1d0)), intent(out)  runoffAGveg_grid,
real(kind(1d0)), intent(out)  runoffPipes_grid,
real(kind(1d0)), intent(out)  runoffWaterBody_grid 
)

Definition at line 398 of file suews_phys_waterdist.f95.

410 IMPLICIT NONE
411
412 !Stores flood water when surface state_id exceeds storage capacity [mm]
413 INTEGER :: is ! surface type
414
415 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! surface fractions
416 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: AddWater_surf !Water from other surfaces (WGWaterDist in SUEWS_ReDistributeWater.f95) [mm]
417 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: state_surf_in !Wetness status of each surface type from previous timestep [mm]
418 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: frac_water2runoff_surf !Fraction of water going to runoff/sub-surface soil (WGWaterDist) [-]
419 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap_surf !Capacity of soil store for each surface [mm]
420 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: StateLimit_surf !Limit for state_id of each surface type [mm] (specified in input files)
421 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: drain_surf !Drainage of each surface type [mm]
422 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: WU_surf !external water use of each surface type [mm]
423 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: ev_surf_in !Evaporation
424 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: soilstore_surf_in !Evaporation
425
426 ! REAL(KIND(1D0)), INTENT(in) :: NonWaterFraction !Capacity of pipes to transfer water
427 REAL(KIND(1D0)), INTENT(in) :: PipeCapacity !Capacity of pipes to transfer water
428 REAL(KIND(1D0)), INTENT(in) :: RunoffToWater !Fraction of surface runoff going to water body
429 REAL(KIND(1D0)), INTENT(in) :: pin !Rain per time interval
430 REAL(KIND(1D0)), INTENT(in) :: addImpervious !Water from impervious surfaces of other grids [mm] for whole surface area
431 REAL(KIND(1D0)), INTENT(in) :: nsh_real !nsh cast as a real for use in calculations
432 REAL(KIND(1D0)), INTENT(in) :: PervFraction ! sum of surface cover fractions for impervious surfaces
433 REAL(KIND(1D0)), INTENT(in) :: addVeg !Water from vegetated surfaces of other grids [mm] for whole surface area
434 REAL(KIND(1D0)), INTENT(in) :: addWaterBody !Water from water surface of other grids [mm] for whole surface area
435 REAL(KIND(1D0)), INTENT(in) :: FlowChange !Difference between the input and output flow in the water body
436
437 ! REAL(KIND(1D0)), INTENT(out) :: runoff_grid !Wetness status of each surface type [mm]
438 ! REAL(KIND(1D0)), INTENT(out) :: state_grid !Wetness status of each surface type [mm]
439 ! REAL(KIND(1D0)), INTENT(out) :: ev_grid !Wetness status of each surface type [mm]
440 ! REAL(KIND(1D0)), INTENT(out) :: surf_chang_grid !Wetness status of each surface type [mm]
441 REAL(KIND(1D0)), INTENT(out) :: runoffAGimpervious_grid !Above ground runoff from impervious surface [mm] for whole surface area
442 REAL(KIND(1D0)), INTENT(out) :: runoffAGveg_grid !Above ground runoff from vegetated surfaces [mm] for whole surface area
443 REAL(KIND(1D0)), INTENT(out) :: runoffPipes_grid !Runoff in pipes [mm] for whole surface area
444 REAL(KIND(1D0)), INTENT(out) :: runoffWaterBody_grid !Above ground runoff from water surface [mm] for whole surface area
445 ! REAL(KIND(1D0)), INTENT(out) :: NWstate_grid !Above ground runoff from water surface [mm] for whole surface area
446
447 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: ev_surf_out !evaporation each surface type [mm]
448 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: state_surf_out !Wetness status of each surface type [mm]
449 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: soilstore_surf_out !Wetness status of each surface type [mm]
450
451 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: runoff_surf !Runoff from each surface type [mm]
452 REAL(KIND(1D0)), DIMENSION(nsurf) :: soilstore !Soil moisture of each surface type [mm]
453 ! REAL(KIND(1D0)), DIMENSION(nsurf) :: chang !Change in state_id [mm]
454 REAL(KIND(1D0)), DIMENSION(2) :: SurplusEvap !Surplus for evaporation in 5 min timestep
455
456 ! =============================
457 ! TS 01 Apr 2022:
458 ! these variables were used as inout variables in the original code for water transfer between grids
459 ! but they are not used in the new code, so they are removed here
460 REAL(KIND(1D0)) :: surplusWaterBody !Extra runoff that goes to water body [mm] as specified by RunoffToWater
461
462 ! REAL(KIND(1D0)) :: p_mm !Inputs to surface water balance
463
464 REAL(KIND(1D0)), DIMENSION(nsurf) :: ev_surf !Evaporation
465
466 !Extra evaporation [mm] from impervious surfaces which cannot happen due to lack of water
467 ! REAL(KIND(1D0)) :: EvPart
468 REAL(KIND(1D0)), PARAMETER :: NotUsed = -55.5
469
470 !Threshold for intense precipitation [mm hr-1]
471 REAL(KIND(1D0)), PARAMETER :: IPThreshold_mmhr = 10 ! NB:this should be an input and can be specified. SG 25 Apr 2018
472
473 ev_surf = ev_surf_in
474 soilstore = soilstore_surf_in
475
476 surplusevap = 0
477 runoffagveg_grid = 0
478 runoffpipes_grid = 0
479 runoffagimpervious_grid = 0
480 runoffwaterbody_grid = 0
481 surpluswaterbody = 0
482
483 ! surf_chang_grid = 0
484 runoff_surf = 0
485 state_surf_out = 0
486 ! state_grid = 0
487 ! NWstate_grid = 0
488 ! print *, 'cal_water_storage_surf: entering'
489
490 DO is = 1, nsurf !For each surface in turn
491 ! print *, 'cal_water_storage_surf: is = ', is
492 IF (sfr_surf(is) > 0) THEN
493
494 !Surface water balance and soil store updates (can modify ev, updates state_id)
495 CALL cal_water_storage( &
496 is, sfr_surf, pipecapacity, runofftowater, pin, & ! input:
497 wu_surf, &
498 drain_surf, addwater_surf, addimpervious, nsh_real, state_surf_in, frac_water2runoff_surf, &
499 pervfraction, addveg, soilstorecap_surf, addwaterbody, flowchange, statelimit_surf, &
500 runoffagimpervious_grid, runoffagveg_grid, runoffpipes_grid, ev_surf(is), soilstore, & ! inout:
501 surpluswaterbody, surplusevap, runoffwaterbody_grid, & ! inout:
502 runoff_surf, state_surf_out) !output:
503 END IF
504
505 END DO !end loop over surfaces
506
507 ! update soilstore
508 soilstore_surf_out = soilstore
509
510 ! update evaporation
511 ev_surf_out = ev_surf
512 ! print *, 'cal_water_storage_surf: exiting'
513

References cal_water_storage(), and allocatearray::nsurf.

Referenced by suews_driver::suews_cal_qe().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ drainage()

subroutine waterdist_module::drainage ( integer, intent(in)  is,
real(kind(1d0)), intent(in)  state_is,
real(kind(1d0)), intent(in)  StorCap,
real(kind(1d0)), intent(in)  DrainEq,
real(kind(1d0)), intent(in)  DrainCoef1,
real(kind(1d0)), intent(in)  DrainCoef2,
real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), intent(out)  drain_is 
)

Definition at line 20 of file suews_phys_waterdist.f95.

29
30 !Calculation of drainage for each land surface.
31 !INPUT: Storage capacity, type of drainage equation used, drainage coefficients
32 ! used in the equation
33 !Modified by HCW 16 Feb 2015
34 ! Removed option of Eq 4 (calculation needs to be checked before re-implementing).
35 ! Code writes an error if calculated drainage exceeds surface state_id (but code continues).
36 ! This may indicate inappropriate drainage equation, storage capacities or model tstep.
37 !Modified by LJ in Aug 2011. Drainage cannot exceed the surface storage.
38 !Modified LJ in 10/2010
39 !------------------------------------------------------------------------------
40
41 IMPLICIT NONE
42 INTEGER, INTENT(in) :: is ! surface type number
43
44 REAL(KIND(1D0)), INTENT(in) :: state_is !Wetness status of surface type "is" [mm]
45 REAL(KIND(1D0)), INTENT(in) :: StorCap !current storage capacity [mm]
46 REAL(KIND(1D0)), INTENT(in) :: DrainCoef1 !Drainage coeff 1 [units depend on choice of eqn]
47 REAL(KIND(1D0)), INTENT(in) :: DrainCoef2 !Drainage coeff 2 [units depend on choice of eqn]
48 REAL(KIND(1D0)), INTENT(in) :: DrainEq !Drainage equation to use
49 REAL(KIND(1D0)), INTENT(in) :: nsh_real !nsh cast as a real for use in calculations
50 REAL(KIND(1D0)), INTENT(out) :: drain_is !Drainage of surface type "is" [mm]
51
52 !If surface is dry, no drainage occurs
53 IF (state_is < 0.000000001) THEN
54 drain_is = 0.0
55 ELSE
56 IF (int(draineq) == 1) THEN !Falk and Niemczynowicz (1978): Drainage equation for paved, buildings and irrigated grass
57
58 IF (state_is < storcap) THEN
59 drain_is = 0 !No drainage if state_id is less than storage capacity
60 ELSE
61 drain_is = (draincoef1*(state_is - storcap)**draincoef2)/nsh_real
62 END IF
63
64 ELSEIF (int(draineq) == 2) THEN !Rutter eqn corrected for c=0, see Eq 9 of Calder & Wright 1986
65 drain_is = (draincoef1*(exp(draincoef2*state_is) - 1))/nsh_real
66 ! N.B. -1 is correct here but brackets are wrong in G&O 1991 Eq 5 & Ja11 Eq 18.
67
68 ELSEIF (int(draineq) == 3) THEN !Falk and Niemczynowicz (1978)
69 drain_is = (draincoef1*(state_is**draincoef2))/nsh_real
70
71 END IF
72
73 ! Check value obtained is physically reasonable
74 ! More water cannot drain than is in the surface state_id
75 ! although high initial rate of drainage tries to drain more water than is in state_id within tstep
76 ! May indicate shorter tstep needed, or a more suitable equation
77 IF (drain_is > state_is) THEN
78 !write(*,*) 'Drainage:', is, drain(is), state_id(is), drain(is)-state_id(is), DrainEq, DrainCoef1, DrainCoef2, nsh_real
79 CALL errorhint(61, 'SUEWS_drain: drain_is > state_is for surface is ', drain_is, state_is, is)
80 drain_is = state_is !All water in state_id is drained (but no more)
81 ELSEIF (drain_is < 0.0001) THEN
82 drain_is = 0
83 END IF
84 END IF
85
86 RETURN
87

References errorhint().

Referenced by suews_driver::suews_cal_water().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ redistributewater()

subroutine waterdist_module::redistributewater ( integer, intent(in)  SnowUse,
real(kind(1d0)), dimension(nsurf + 1, nsurf - 1), intent(in)  WaterDist,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  Drain,
real(kind(1d0)), dimension(nsurf), intent(out)  AddWaterRunoff,
real(kind(1d0)), dimension(nsurf), intent(out)  AddWater 
)

Definition at line 773 of file suews_phys_waterdist.f95.

776 !Drainage moves into different parts defined by WaterDistSS_YYYY.txt. LJ 2010
777 !AddWater(is) is that amount of water that is gained for each surface
778 !Latest update takes snow into account. 22/03/2013 LJ
779 !-------------------------------------------------------------------
780
781 IMPLICIT NONE
782 INTEGER, INTENT(in) :: SnowUse !Snow part used (1) or not used (0)
783
784 REAL(KIND(1D0)), INTENT(in) :: WaterDist(nsurf + 1, nsurf - 1) !Within-grid water distribution to other surfaces and runoff/soil store [-]
785 REAL(KIND(1D0)), INTENT(in) :: sfr_surf(nsurf) !Surface fractions [-]
786 REAL(KIND(1D0)), INTENT(in) :: Drain(nsurf) !Drainage of each surface type [mm]
787
788 REAL(KIND(1D0)), INTENT(out) :: AddWaterRunoff(nsurf) !Fraction of water going to runoff/sub-surface soil (WGWaterDist) [-]
789 REAL(KIND(1D0)), INTENT(out) :: AddWater(nsurf) !Water from other surfaces (WGWaterDist in SUEWS_ReDistributeWater.f95) [mm]
790
791 INTEGER :: ii, jj
792 INTEGER :: NSurfDoNotReceiveDrainage = 0 !Number of surfaces that do not receive drainage water (green roof)
793
794 !Fractions that go to runoff from each surface
795 DO ii = 1, nsurf - 1 !not water in the calculation
796 addwaterrunoff(ii) = waterdist(8, ii)
797 END DO
798 addwaterrunoff(watersurf) = 0
799 addwater = 0
800
801 DO ii = 1, nsurf - nsurfdonotreceivedrainage !go through surfaces from 1 to 7. These gain water through drainage
802 DO jj = 1, nsurf - (nsurfdonotreceivedrainage + 1) !From where surface ii can gain water - can't gain water from itself
803
804 IF (sfr_surf(ii) /= 0) THEN !Water movement takes place only if surface fraction exists
805
806 !No snow calculations!
807 IF (snowuse == 0) THEN
808 addwater(ii) = addwater(ii) + (drain(jj)*sfr_surf(jj)/sfr_surf(ii))*waterdist(ii, jj) !Original
809
810 !Snow included, This needs to be fixed at some point. LJ Mar 2013
811 ELSE
812 addwaterrunoff(jj) = addwaterrunoff(jj) + waterdist(ii, jj) !No receiving surface -> runoff
813 END IF
814
815 ELSE
816 addwaterrunoff(jj) = addwaterrunoff(jj) + waterdist(ii, jj) !If no receiving surface exists,
817 !water fraction goes to AddWaterRunoff
818 END IF
819 END DO
820 END DO
821

References allocatearray::watersurf.

Referenced by suews_driver::suews_cal_water().

Here is the caller graph for this function:

◆ suews_cal_horizontalsoilwater()

subroutine waterdist_module::suews_cal_horizontalsoilwater ( real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilDepth,
real(kind(1d0)), dimension(nsurf), intent(in)  SatHydraulicConduct,
real(kind(1d0)), intent(in)  SurfaceArea,
real(kind(1d0)), intent(in)  NonWaterFraction,
real(kind(1d0)), intent(in)  tstep_real,
real(kind(1d0)), dimension(nsurf), intent(inout)  soilstore_id,
real(kind(1d0)), dimension(nsurf), intent(out)  runoffSoil,
real(kind(1d0)), intent(out)  runoffSoil_per_tstep 
)

Definition at line 948 of file suews_phys_waterdist.f95.

960 !Transfers water in soil stores of land surfaces LJ (2010)
961 !Change the model to use varying hydraulic conductivity instead of constant value LJ (7/2011)
962 !If one of the surface's soildepth is zero, no water movement is considered
963 ! LJ 15/06/2017 Modification: - Moved location of runoffSoil_per_tstep within previous if-loop to avoid dividing with zero with 100% water surface
964 ! HCW 22/02/2017 Modifications: - Minor bug fixed in VWC1/B_r1 comparison - if statements reversed
965 ! HCW 13/08/2014 Modifications: - Order of surfaces reversed (for both is and jj loops)
966 ! - Number of units (e.g. properties) added to distance calculation
967 ! HCW 12/08/2014 Modifications: - Distance changed from m to mm in dI_dt calculation
968 ! - dI_dt [mm s-1] multiplied by no. seconds in timestep -> dI [mm]
969 ! - if MatPot is set to max. value (100000 mm), Km set to 0 mm s-1
970 ! - Provide parameters for residual volumetric soil moisture [m3 m-3]
971 ! (currently hard coded as 0.1 m3 m-3 for testing)
972 !
973 !------------------------------------------------------
974 ! use SUES_data
975 ! use gis_data
976 ! use time
977 ! use allocateArray
978
979 IMPLICIT NONE
980
981 REAL(KIND(1D0)), INTENT(in) :: sfr_surf(nsurf) ! surface fractions
982 REAL(KIND(1D0)), INTENT(in) :: SoilStoreCap(nsurf) !Capacity of soil store for each surface [mm]
983 REAL(KIND(1D0)), INTENT(in) :: SoilDepth(nsurf) !Depth of sub-surface soil store for each surface [mm]
984 REAL(KIND(1D0)), INTENT(in) :: SatHydraulicConduct(nsurf) !Saturated hydraulic conductivity for each soil subsurface [mm s-1]
985 REAL(KIND(1D0)), INTENT(in) :: SurfaceArea !Surface area of the study area [m2]
986 REAL(KIND(1D0)), INTENT(in) :: NonWaterFraction ! sum of surface cover fractions for all except water surfaces
987 REAL(KIND(1D0)), INTENT(in) :: tstep_real !tstep cast as a real for use in calculations
988
989 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: soilstore_id !Soil moisture of each surface type [mm]
990 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: runoffSoil !Soil runoff from each soil sub-surface [mm]
991
992 REAL(KIND(1D0)), INTENT(out) :: runoffSoil_per_tstep !Runoff to deep soil per timestep [mm] (for whole surface, excluding water body)
993
994 INTEGER :: jj, is
995 REAL(KIND(1D0)) :: &
996 dimenwatercon1, dimenwatercon2, &
997 soilmoistcap_vol1, &
998 soilmoist_vol1, &
999 soilmoistcap_vol2, &
1000 soilmoist_vol2, &
1001 b_r1, matpot1, km1, &
1002 b_r2, matpot2, km2, &
1003 distance, kmweight, di, &
1004 di_dt !Water flow between two stores
1005
1006 REAL(KIND(1D0)), PARAMETER :: &
1007 alphavg = 0.0005, & !Set alphavG to match value in van Genuchten (1980) [mm-1]
1008 nunits = 1 !Can change to represent plot/base unit size
1009
1010 ! SoilMoist_vol1,2 = Volumetric soil moisture [m3 m-3]
1011 ! SoilMoistCap_vol1,2 = Volumetric soil moisture capacity [m3 m-3] (from FunctionalTypes)
1012 ! MatPot1,2 = Water potential (i.e. pressure head) of store [mm]
1013 ! DimenWaterCon1,2 = Dimensionless water content, or relative saturation [-]
1014 ! Distance = Distance between two stores [m]
1015 ! B_r1,2 = Residual volumetric soil moisture [m3 m-3]
1016 ! Km1,2 = Hydraulic conductivity of store [mm s-1]
1017 ! KmWeight = Weighted hydraulic conductivity [mm s-1]
1018 ! alphavG = Parameter (could depend on soil texture) [mm-1]
1019 ! dI = Water flow between stores [mm] dI = dI_dt * no. secs in each timestep
1020 ! if dI > 0, first surface gains water, second surface loses water
1021 ! NUnits = Number of repeating units (e.g. properties, blocks) for distance calculation [-]
1022 runoffsoil = 0
1023 runoffsoil_per_tstep = 0
1024
1025 DO is = 1, nsurf - 1 !nsurf-1,1,-1 !Loop through each surface, excluding water surface (runs backwards as of 13/08/2014, HCW)
1026
1027 IF (sfr_surf(is) /= 0 .AND. soilstorecap(is) > 0) THEN !If particular surface area exists
1028 ! and is capable of storing water (SoilStoreCap [mm])
1029 DO jj = is + 1, nsurf - 1 !is-1,1,-1 !Sub-loop through remaining surfaces (runs backwards as of 13/08/2014, HCW)
1030
1031 IF (sfr_surf(jj) /= 0 .AND. soilstorecap(jj) > 0) THEN !If other surface area exists
1032 ! and is capable of storing water
1033
1034 ! ---- For surface 1 -----------------------------------------------------
1035 ! Calculate non-saturated VWC
1036 soilmoistcap_vol1 = soilstorecap(is)/soildepth(is) !Volumetric soil moisture capacity [m3 m-3] (i.e. saturated VWC)
1037 soilmoist_vol1 = soilstore_id(is)/soildepth(is) !Volumetric soil moisture [m3 m-3]
1038
1039 !B_r1=SoilMoistCap_Vol1-SoilMoist_vol1 !Residual soil moisture content [m3 m-3]
1040 b_r1 = 0.1 !HCW 12/08/2014 Temporary fix
1041 ! Need to add residual soil moisture values to FunctionalTypes
1042 !B_r1=VolSoilMoistRes(is) !Residual soil moisture content [m3 m-3]
1043
1044 !Order of if statements reversed HCW 22 Feb 2017
1045 !If soil moisture less than or equal to residual value, set MatPot to max and Km to 0 to suppress water movement
1046 IF (b_r1 >= soilmoist_vol1) THEN
1047 matpot1 = 100000
1048 km1 = 0 !Added by LJ in Nov 2013
1049 ! Otherwise, there should be enough water in the soil to allow horizontal transfer
1050 ELSE
1051 dimenwatercon1 = (soilmoist_vol1 - b_r1)/(soilmoistcap_vol1 - b_r1) !Dimensionless water content [-]
1052
1053 ! If very large or very small, adjust for calculation of MatPot and Km
1054 IF (dimenwatercon1 > 0.99999) THEN
1055 dimenwatercon1 = dimenwatercon1 - 0.0001 !This cannot equal 1
1056 END IF
1057
1058 IF (dimenwatercon1 < 0.00000005) THEN
1059 dimenwatercon1 = dimenwatercon1 + 0.0000001 !Added HCW 22 Feb 2017
1060 END IF
1061
1062 !van Genuchten (1980), with n=2 and m = 1-1/n = 1/2
1063 !Water potential of first store [mm] (van Genuchten 1980, Eq 3 rearranged)
1064 matpot1 = sqrt(1/dimenwatercon1**2 - 1)/alphavg
1065
1066 !Hydraulic conductivity of first store [mm s-1] (van Genuchten 1980, Eq 8)
1067 km1 = sathydraulicconduct(is)*sqrt(dimenwatercon1)*(1 - (1 - dimenwatercon1**2)**0.5)**2
1068
1069 !Check this value (HCW 12/08/2014)
1070 IF (matpot1 > 100000) THEN
1071 matpot1 = 100000 !Max. potential is 100000 mm (van Genuchten 1980)
1072 km1 = 0 !Added by HCW 12/08/2014
1073 END IF
1074
1075 END IF
1076
1077 ! ---- For surface 2 -----------------------------------------------------
1078 ! Calculate non-saturated VWC
1079 soilmoistcap_vol2 = soilstorecap(jj)/soildepth(jj) !Volumetric soil moisture capacity [m3 m-3] (i.e. saturated VWC)
1080 soilmoist_vol2 = soilstore_id(jj)/soildepth(jj) !Volumetric soil moisture [m3 m-3]
1081
1082 !B_r2=SoilMoistCap_Vol2-SoilMoist_vol2 !Residual soil moisture content [m3 m-3]
1083 b_r2 = 0.1 !HCW 12/08/2014 Temporary fix
1084 ! Need to add residual soil moisture values to FunctionalTypes
1085 !B_r2=VolSoilMoistRes(jj) !Residual soil moisture content [m3 m-3]
1086
1087 !If soil moisture below residual value, set MatPot to maximum
1088 IF (b_r2 >= soilmoist_vol2) THEN
1089 matpot2 = 100000
1090 km2 = 0 !Added by LJ in Nov 2013
1091 ELSE
1092 dimenwatercon2 = (soilmoist_vol2 - b_r2)/(soilmoistcap_vol2 - b_r2) !Dimensionless water content [-]
1093
1094 IF (dimenwatercon2 > 0.99999) THEN
1095 dimenwatercon2 = dimenwatercon2 - 0.0001 !This cannot equal 1
1096 END IF
1097
1098 IF (dimenwatercon2 < 0.00000005) THEN
1099 dimenwatercon2 = dimenwatercon2 + 0.0000001 !Added HCW 22 Feb 2017
1100 END IF
1101 ! print*, 'soilstore_id = ', soilstore_id
1102 ! print*, 'SoilStoreCap = ', SoilStoreCap
1103 ! print*, 'SoilDepth = ', SoilDepth
1104 ! print*, 'SoilMoist_vol2 = ', SoilMoist_vol2
1105 ! print*, 'SoilMoistCap_Vol2 = ', SoilMoistCap_Vol2
1106 ! print*, 'is = ', is, ' jj = ', jj, ' DimenWaterCon2 = ', DimenWaterCon2
1107 !van Genuchten (1980), with n=2 and m = 1-1/n = 1/2
1108 !Water potential of second store [mm] (van Genuchten 1980, Eq 3 rearranged)
1109 matpot2 = sqrt(1/dimenwatercon2**2 - 1)/alphavg
1110
1111 !Hydraulic conductivity of second store [mm s-1] (van Genuchten 1980, Eq 8)
1112 km2 = sathydraulicconduct(jj)*sqrt(dimenwatercon2)*(1 - (1 - dimenwatercon2**2)**0.5)**2
1113
1114 IF ((matpot2) > 100000) THEN
1115 matpot2 = 100000 !Max. potential is 100000 mm (van Genuchten 1980)
1116 km2 = 0 !Added by HCW 12/08/2014
1117 END IF
1118
1119 END IF
1120
1121 ! ------------------------------------------------------------------------
1122
1123 !Find distance between the two stores (see Jarvi et al. 2011)
1124 !SurfaceArea in m2 (changed from ha to m2 n SUEWS_Initial), so Distance in m
1125 distance = (sqrt(sfr_surf(is)*surfacearea/nunits) + sqrt(sfr_surf(jj)*surfacearea/nunits))/2
1126
1127 !Calculate areally-weighted hydraulic conductivity [mm s-1]
1128 kmweight = (sfr_surf(is)*km1 + sfr_surf(jj)*km2)/(sfr_surf(is) + sfr_surf(jj))
1129
1130 !Find water flow between the two stores [mm s-1] (Green-Ampt equation, Hillel 1971)
1131 !Multiply Distance by 1000 to convert m to mm (HCW 12/08/2014)
1132 di_dt = -(kmweight)*(-matpot1 + matpot2)/(distance*1000)
1133
1134 !Multiply dI_dt by number of seconds in timestep to convert mm s-1 to mm
1135 !Added by HCW 12/08/2014
1136 di = di_dt*tstep_real !Use dI instead of dI_dt in the following calculations
1137
1138 !Move water (in mm) ------------------------------------------------------
1139 !Water moves only if (i) there is sufficient water to move and (ii) there is space to move it
1140
1141 ! If there is sufficient water in both surfaces, allow movement of dI to occur
1142 IF ((soilstore_id(jj) >= di*sfr_surf(is)/sfr_surf(jj)) .AND. ((soilstore_id(is) + di) >= 0)) THEN
1143 soilstore_id(is) = soilstore_id(is) + di
1144 soilstore_id(jj) = soilstore_id(jj) - di*sfr_surf(is)/sfr_surf(jj) !Check (HCW 13/08/2014) - QUESTION: why adjust for jj and not is?
1145
1146 ! If insufficient water in first surface to move dI, instead move as much as possible
1147 ELSEIF ((soilstore_id(is) + di) < 0) THEN
1148 soilstore_id(jj) = soilstore_id(jj) + soilstore_id(is)*sfr_surf(is)/sfr_surf(jj) !HCW 12/08/2014 switched order of these two lines
1149 soilstore_id(is) = 0 !Check (HCW 13/08/2014) - QUESTION: can SM actually go to zero, or is this inconsistent with SMres?
1150
1151 ! If insufficient water in second surface to move dI, instead move as much as possible
1152 ELSE
1153 soilstore_id(is) = soilstore_id(is) + soilstore_id(jj)*sfr_surf(jj)/sfr_surf(is)
1154 soilstore_id(jj) = 0
1155 END IF
1156
1157 !If soil moisture exceeds capacity, excess goes to soil runoff (first surface)
1158 IF (soilstore_id(is) > soilstorecap(is)) THEN
1159 runoffsoil(is) = runoffsoil(is) + (soilstore_id(is) - soilstorecap(is))
1160 soilstore_id(is) = soilstorecap(is)
1161 !elseif (soilstore_id(is)<0) then !HCW 13/08/2014 commented out as should never be true here anyway...
1162 ! soilstore_id(is)=0 ! ... and if so, need to do more here (i.e. account for other water too)
1163 END IF
1164
1165 !If soil moisture exceeds capacity, excess goes to soil runoff (second surface)
1166 IF (soilstore_id(jj) > soilstorecap(jj)) THEN
1167 runoffsoil(jj) = runoffsoil(jj) + (soilstore_id(jj) - soilstorecap(jj))
1168 soilstore_id(jj) = soilstorecap(jj)
1169 !elseif (soilstore_id(jj)<0) then !HCW 13/08/2014 commented out (as above)
1170 ! soilstore_id(jj)=0
1171 END IF
1172
1173 END IF !end if second surface exists and is capable of storing water
1174
1175 END DO !end jj loop over second surface
1176
1177 runoffsoil_per_tstep = runoffsoil_per_tstep + (runoffsoil(is)*sfr_surf(is)/nonwaterfraction) !Excludes water body. Moved here as otherwise code crashed when NonWaterFraction=0
1178
1179 END IF !end if first surface exists and is capable of storing water
1180
1181 !runoffSoil_per_tstep=runoffSoil_per_tstep+(runoffSoil(is)*sfr_surf(is)/NonWaterFraction) !Excludes water body
1182
1183 END DO !is loop over first surface
1184

Referenced by suews_driver::suews_cal_main().

Here is the caller graph for this function:

◆ suews_cal_soilstate()

subroutine waterdist_module::suews_cal_soilstate ( integer, intent(in)  SMDMethod,
real(kind(1d0)), intent(in)  xsmd,
real(kind(1d0)), intent(in)  NonWaterFraction,
real(kind(1d0)), intent(in)  SoilMoistCap,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap,
real(kind(1d0)), intent(in)  surf_chang_per_tstep,
real(kind(1d0)), dimension(nsurf), intent(in)  soilstore_id,
real(kind(1d0)), dimension(nsurf), intent(in)  soilstoreOld,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(out)  smd,
real(kind(1d0)), dimension(nsurf), intent(out)  smd_nsurf,
real(kind(1d0)), intent(out)  tot_chang_per_tstep,
real(kind(1d0)), intent(out)  SoilState 
)

Definition at line 879 of file suews_phys_waterdist.f95.

884
885 IMPLICIT NONE
886 ! INTEGER, PARAMETER :: nsurf = 7
887
888 INTEGER, INTENT(in) :: SMDMethod
889 REAL(KIND(1D0)), INTENT(in) :: xsmd
890 REAL(KIND(1D0)), INTENT(in) :: NonWaterFraction
891 REAL(KIND(1D0)), INTENT(in) :: SoilMoistCap
892
893 REAL(KIND(1D0)), INTENT(in) :: surf_chang_per_tstep
894 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: soilstore_id !Soil moisture of each surface type [mm]
895 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: soilstoreOld !Soil moisture of each surface type from previous timestep [mm]
896 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
897 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap !Capacity of soil store for each surface [mm]
898
899 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: smd_nsurf !smd for each surface
900 REAL(KIND(1D0)), INTENT(out) :: SoilState !Area-averaged soil moisture [mm] for whole surface
901 REAL(KIND(1D0)), INTENT(out) :: smd !One value for whole surface
902 REAL(KIND(1D0)), INTENT(out) :: tot_chang_per_tstep !Change in surface state_id
903
904 REAL(KIND(1D0)), PARAMETER :: NotUsed = -999
905 REAL(KIND(1D0)), PARAMETER :: NAN = -999
906 INTEGER :: is
907
908 soilstate = 0 !Area-averaged soil moisture [mm] for whole surface
909 IF (nonwaterfraction /= 0) THEN !Fixed for water surfaces only
910 DO is = 1, nsurf - 1 !No water body included
911 soilstate = soilstate + (soilstore_id(is)*sfr_surf(is)/nonwaterfraction)
912 IF (soilstate < 0) THEN
913 CALL errorhint(62, 'SUEWS_Calculations: total SoilState < 0 (just added surface is) ', soilstate, notused, is)
914 ELSEIF (soilstate > soilmoistcap) THEN
915 CALL errorhint(62, 'SUEWS_Calculations: total SoilState > capacity (just added surface is) ', soilstate, notused, is)
916 !SoilMoist_state=SoilMoistCap !What is this LJ 10/2010 - QUESTION: SM exceeds capacity, but where does extra go?HCW 11/2014
917 END IF
918 END DO !end loop over surfaces
919 ! SoilState = DOT_PRODUCT(soilstore_id(1:nsurf - 1), sfr_surf(1:nsurf - 1))/NonWaterFraction
920 ! IF (SoilState < 0) THEN
921 ! CALL ErrorHint(62, 'SUEWS_Calculations: total SoilState < 0 (just added surface is) ', SoilState, NotUsed, is)
922 ! ELSEIF (SoilState > SoilMoistCap) THEN
923 ! CALL ErrorHint(62, 'SUEWS_Calculations: total SoilState > capacity (just added surface is) ', SoilState, NotUsed, is)
924 ! !SoilMoist_state=SoilMoistCap !What is this LJ 10/2010 - QUESTION: SM exceeds capacity, but where does extra go?HCW 11/2014
925 ! ENDIF
926 END IF
927
928 ! Calculate soil moisture deficit
929 smd = soilmoistcap - soilstate !One value for whole surface
930 smd_nsurf = soilstorecap - soilstore_id !smd for each surface
931
932 ! Soil stores can change after horizontal water movements
933 ! Calculate total change in surface and soil state_id
934 tot_chang_per_tstep = surf_chang_per_tstep !Change in surface state_id
935 DO is = 1, (nsurf - 1) !No soil for water surface (so change in soil moisture is zero)
936 tot_chang_per_tstep = tot_chang_per_tstep + ((soilstore_id(is) - soilstoreold(is))*sfr_surf(is)) !Add change in soil state_id
937 END DO
938
939 IF (smdmethod > 0) THEN ! use observed value
940 ! smd_nsurf=NAN
941 smd_nsurf = nan
942 smd = xsmd
943 END IF
944

References errorhint(), and allocatearray::nsurf.

Referenced by suews_driver::suews_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ suews_cal_wateruse()

subroutine waterdist_module::suews_cal_wateruse ( real(kind(1d0)), intent(in)  nsh_real,
real(kind(1d0)), intent(in)  wu_m3,
real(kind(1d0)), intent(in)  SurfaceArea,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(in)  IrrFracPaved,
real(kind(1d0)), intent(in)  IrrFracBldgs,
real(kind(1d0)), intent(in)  IrrFracEveTr,
real(kind(1d0)), intent(in)  IrrFracDecTr,
real(kind(1d0)), intent(in)  IrrFracGrass,
real(kind(1d0)), intent(in)  IrrFracBSoil,
real(kind(1d0)), intent(in)  IrrFracWater,
integer, dimension(3), intent(in)  DayofWeek_id,
real(kind(1d0)), dimension(0:23, 2), intent(in)  WUProfA_24hr,
real(kind(1d0)), dimension(0:23, 2), intent(in)  WUProfM_24hr,
real(kind(1d0)), intent(in)  InternalWaterUse_h,
real(kind(1d0)), dimension(12), intent(in)  HDD_id,
real(kind(1d0)), dimension(9), intent(in)  WUDay_id,
integer, intent(in)  WaterUseMethod,
integer, intent(in)  NSH,
integer, intent(in)  it,
integer, intent(in)  imin,
integer, intent(in)  DLS,
real(kind(1d0)), dimension(nsurf), intent(out)  wu_surf,
real(kind(1d0)), intent(out)  wu_int,
real(kind(1d0)), intent(out)  wu_ext 
)

Definition at line 1189 of file suews_phys_waterdist.f95.

1199 ! Conversion of water use (irrigation)
1200 ! Last modified:
1201 ! TS 30 Nov 2019 - allow external water use for all surfaces
1202 ! (previously only on vegetated surfaces)
1203 ! TS 30 Oct 2018 - fixed a bug in external water use
1204 ! TS 08 Aug 2017 - addded explicit interface
1205 ! LJ 6 Apr 2017 - WUchoice changed to WaterUseMethod
1206 ! TK 14 Mar 2017 - Corrected the variable name WUAreaEveTr_m2 -> WUAreaGrass_m2 (row 35)
1207 ! Corrected conversion from m to mm /1000 -> *1000 (row 47 and 60)
1208 ! LJ 27 Jan 2016 - Removing Tab:s and cleaning the code
1209 ! HCW 12 Feb 2015 - Water use [mm] now inidcates the amount of water supplied for each surface
1210 ! HCW 26 Jan 2015 - Water use [mm] is the same for each surface at the moment and indicates the
1211 ! amount of water supplied for each irrigated area
1212 !
1213 ! To Do:
1214 ! - Add functionality for water on paved surfaces (street cleaning, fountains)
1215
1216 IMPLICIT NONE
1217 ! INTEGER, PARAMETER :: nsurf = 7
1218
1219 REAL(KIND(1D0)), INTENT(in) :: nsh_real
1220 REAL(KIND(1D0)), INTENT(in) :: wu_m3 ! external water input (e.g., irrigation) [m3]
1221 REAL(KIND(1D0)), INTENT(in) :: SurfaceArea !Surface area of the study area [m2]
1222 REAL(KIND(1D0)), INTENT(IN) :: IrrFracPaved !Fraction of paved which are irrigated
1223 REAL(KIND(1D0)), INTENT(IN) :: IrrFracBldgs !Fraction of buildings (e.g., green roofs) which are irrigated
1224 REAL(KIND(1D0)), INTENT(IN) :: IrrFracEveTr !Fraction of evergreen trees which are irrigated
1225 REAL(KIND(1D0)), INTENT(IN) :: IrrFracDecTr !Fraction of deciduous trees which are irrigated
1226 REAL(KIND(1D0)), INTENT(IN) :: IrrFracGrass !Fraction of grass which is irrigated
1227 REAL(KIND(1D0)), INTENT(IN) :: IrrFracBSoil !Fraction of bare soil trees which are irrigated
1228 REAL(KIND(1D0)), INTENT(IN) :: IrrFracWater !Fraction of water which are irrigated
1229 REAL(KIND(1D0)), INTENT(in) :: InternalWaterUse_h !Internal water use [mm h-1]
1230 REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: WUProfA_24hr !Automatic water use profiles at hourly scales
1231 REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: WUProfM_24hr !Manual water use profiles at hourly scales
1232 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf !Surface fractions [-]
1233
1234 REAL(KIND(1D0)), DIMENSION(12), INTENT(in) :: HDD_id !HDD(id-1), Heating Degree Days (see SUEWS_DailyState.f95)
1235 REAL(KIND(1D0)), DIMENSION(9), INTENT(in) :: WUDay_id !WUDay(id-1), Daily water use for EveTr, DecTr, Grass [mm] (see SUEWS_DailyState.f95)
1236
1237 INTEGER, INTENT(in) :: DayofWeek_id(3) !DayofWeek(id) 1 - day of week; 2 - month; 3 - season
1238 INTEGER, INTENT(in) :: WaterUseMethod !Use modelled (0) or observed (1) water use
1239 INTEGER, INTENT(in) :: NSH !Number of timesteps per hour
1240 INTEGER, INTENT(in) :: it !Hour
1241 INTEGER, INTENT(in) :: imin !Minutes
1242 INTEGER, INTENT(in) :: DLS !day lightsavings =1 + 1h) =0
1243
1244 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: wu_surf !external Water use for each surface [mm]
1245 REAL(KIND(1D0)), INTENT(out) :: wu_int !Internal water use for the model timestep [mm] (over whole study area)
1246 REAL(KIND(1D0)), INTENT(out) :: wu_ext !External water use for the model timestep [mm] (over whole study area)
1247
1248 REAL(KIND(1D0)) :: wu_EveTr !Water use for evergreen trees/shrubs [mm]
1249 REAL(KIND(1D0)) :: wu_DecTr !Water use for deciduous trees/shrubs [mm]
1250 REAL(KIND(1D0)) :: wu_Grass !Water use for grass [mm]
1251
1252 REAL(KIND(1D0)), DIMENSION(nsurf) :: WUDay_A_id !modelled Automatic Daily water use for each surface [mm] (see SUEWS_DailyState.f95)
1253 REAL(KIND(1D0)), DIMENSION(nsurf) :: WUDay_M_id !modelled Manual Daily water use for each surface [mm] (see SUEWS_DailyState.f95)
1254 REAL(KIND(1D0)), DIMENSION(nsurf) :: IrrFrac !faction of irrigated part in each surface [-]
1255 REAL(KIND(1D0)), DIMENSION(nsurf) :: WUArea !water use area [m2] for each surface type
1256
1257 ! REAL(KIND(1d0)):: WUAreaEveTr_m2
1258 ! REAL(KIND(1d0)):: WUAreaDecTr_m2
1259 ! REAL(KIND(1d0)):: WUAreaGrass_m2
1260 REAL(KIND(1D0)) :: WUAreaTotal_m2
1261 REAL(KIND(1D0)) :: InternalWaterUse !Internal water use for the model timestep [mm]
1262 REAL(KIND(1D0)) :: flag_WuM = 1
1263 REAL(KIND(1D0)) :: wu !Water use for the model timestep [mm]
1264 INTEGER :: ih !Hour corrected for Daylight savings
1265 INTEGER :: iu !1=weekday OR 2=weekend
1266 INTEGER :: tstep ! timestep in second
1267 REAL(KIND(1D0)), PARAMETER :: NAN = -999.
1268 REAL(KIND(1D0)) :: OverUse
1269 REAL(KIND(1D0)) :: rain_cum_daily ! accumulated daily rainfall
1270
1271 REAL(KIND(1D0)) :: get_Prof_SpecTime_sum
1272
1273 REAL(KIND(1D0)) :: WUProfA_tstep ! automatic water use profile value at tstep
1274 REAL(KIND(1D0)) :: WUProfM_tstep ! mannual water use profile value at tstep
1275
1276 ! NB: set OverUse as 0 as done module_constants, TS 22 Oct 2017
1277 ! and the logic for calculating OverUse to be determined
1278 overuse = 0
1279
1280 ! initialise wu
1281 wu = 0
1282
1283 ! timestep in second
1284 tstep = int(3600/nsh)
1285
1286 ! accumulated daily rainfall
1287 rain_cum_daily = hdd_id(11)
1288
1289 ! Irrigated Fraction of each surface
1290 ! TS: as of 20191130, assuming irrigation fraction as ONE except for vegetated surfaces
1291
1292 ! Irrigated Fraction of each surface
1293 ! TS: 20200409, add irrigation fractions for all surfaces
1294 irrfrac = [irrfracpaved, irrfracbldgs, &
1295 irrfracevetr, irrfracdectr, irrfracgrass, &
1296 irrfracbsoil, irrfracwater]
1297
1298 ! --------------------------------------------------------------------------------
1299 ! If water used is observed and provided in the met forcing file, units are m3
1300 ! Divide observed water use (in m3) by water use area to find water use (in mm)
1301 IF (waterusemethod == 1) THEN !If water use is observed
1302 ! Calculate water use area [m2] for each surface type
1303 ! WUAreaEveTr_m2 = IrrFracEveTr*sfr_surf(ConifSurf)*SurfaceArea
1304 ! WUAreaDecTr_m2 = IrrFracDecTr*sfr_surf(DecidSurf)*SurfaceArea
1305 ! WUAreaGrass_m2 = IrrFracGrass*sfr_surf(GrassSurf)*SurfaceArea
1306 ! WUAreaTotal_m2 = WUAreaEveTr_m2 + WUAreaDecTr_m2 + WUAreaGrass_m2
1307
1308 wuarea = irrfrac*sfr_surf*surfacearea
1309 wuareatotal_m2 = sum(wuarea)
1310
1311 !Set water use [mm] for each surface type to zero initially
1312 wu_evetr = 0
1313 wu_dectr = 0
1314 wu_grass = 0
1315
1316 wu_surf = 0
1317 IF (wu_m3 == nan .OR. wu_m3 == 0) THEN !If no water use
1318 ! wu_m3=0
1319 wu = 0
1320 ELSE !If water use
1321 IF (wuareatotal_m2 > 0) THEN
1322 wu = (wu_m3/wuareatotal_m2*1000) !Water use in mm for the whole irrigated area
1323 ! IF (WUAreaEveTr_m2 > 0) THEN
1324 ! wu_EveTr = wu !Water use for Irr EveTr in mm - these are all the same at the moment
1325 ! wu_EveTr = wu_EveTr*IrrFracEveTr !Water use for EveTr in mm
1326 ! ENDIF
1327 ! IF (WUAreaDecTr_m2 > 0) THEN
1328 ! wu_DecTr = wu !Water use for Irr DecTr in mm - these are all the same at the moment
1329 ! wu_DecTr = wu_DecTr*IrrFracDecTr !Water use for DecTr in mm
1330 ! ENDIF
1331 ! IF (WUAreaGrass_m2 > 0) THEN
1332 ! wu_Grass = wu !Water use for Irr Grass in mm - these are all the same at the moment
1333 ! wu_Grass = wu_Grass*IrrFracGrass !Water use for Grass in mm
1334 ! ENDIF
1335
1336 wu_surf = wu*irrfrac
1337
1338 wu = (wu_m3/surfacearea*1000) !Water use for the whole study area in mm
1339 END IF
1340 END IF
1341
1342 ! --------------------------------------------------------------------------------
1343 ! If water use is modelled, calculate at timestep of model resolution [mm]
1344 ELSEIF (waterusemethod == 0) THEN !If water use is modelled
1345
1346 ! Account for Daylight saving
1347 ih = it - dls
1348 IF (ih < 0) ih = 23
1349
1350 ! Weekday or weekend profile
1351 iu = 1 !Set to 1=weekday
1352 ! IF(DayofWeek(id,1)==1.OR.DayofWeek(id,1)==7) THEN
1353 IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) THEN
1354 iu = 2 !Set to 2=weekend
1355 END IF
1356
1357 !write(*,*) (NSH*(ih+1-1)+imin*NSH/60+1)
1358 wuday_a_id = 0
1359 wuday_a_id(conifsurf) = wuday_id(2)
1360 wuday_a_id(decidsurf) = wuday_id(5)
1361 wuday_a_id(grasssurf) = wuday_id(8)
1362
1363 wuday_m_id = 0
1364 wuday_m_id(conifsurf) = wuday_id(3)
1365 wuday_m_id(decidsurf) = wuday_id(6)
1366 wuday_m_id(grasssurf) = wuday_id(9)
1367
1368 ! ---- Automatic irrigation ----
1369 wuprofa_tstep = get_prof_spectime_sum(ih, imin, 0, wuprofa_24hr(:, iu), tstep)
1370
1371 ! ---- Manual irrigation ----
1372 flag_wum = 1 !Initialize flag_WuM to 1, but if raining, reduce manual fraction of water use
1373 ! If cumulative daily precipitation exceeds 2 mm
1374 IF (rain_cum_daily > 2) THEN !.and.WUDay(id-1,3)>0) then !Commented out HCW 23/01/2015
1375 flag_wum = 0 ! 0 -> No manual irrigation if raining
1376 END IF
1377
1378 ! Add manual to automatic to find total irrigation
1379 wuprofm_tstep = get_prof_spectime_sum(ih, imin, 0, wuprofm_24hr(:, iu), tstep)
1380
1381 ! sum up irrigation amount of automatic and manual approaches
1382 wu_surf = wuprofa_tstep*wuday_a_id + wuprofm_tstep*wuday_m_id*flag_wum
1383 ! apply irrigation fraction: part of land covers are not irrigated
1384 wu_surf = wu_surf*irrfrac
1385
1386 ! Total water use for the whole study area [mm]
1387 ! wu = wu_EveTr*sfr_surf(ConifSurf) + wu_DecTr*sfr_surf(DecidSurf) + wu_Grass*sfr_surf(GrassSurf)
1388 wu = dot_product(wu_surf, sfr_surf)
1389
1390 END IF !End WU_choice
1391 ! --------------------------------------------------------------------------------
1392
1393 ! Internal water use is supplied in SUEWS_Irrigation in mm h-1
1394 ! Convert to mm for the model timestep
1395 internalwateruse = internalwateruse_h/nsh_real
1396
1397 ! Remove InternalWaterUse from the total water use
1398 wu_ext = wu - (internalwateruse + overuse)
1399 ! Check ext_wu cannot be negative
1400 IF (wu_ext < 0) THEN
1401 overuse = abs(wu_ext)
1402 wu_ext = 0
1403 ELSE
1404 overuse = 0
1405 END IF
1406
1407 wu_int = wu - wu_ext
1408
1409 ! Decrease the water use for each surface by the same proportion
1410 IF (wu_ext /= 0 .AND. wu /= 0) THEN
1411 ! wu_EveTr = wu_EveTr*wu_ext/wu
1412 ! wu_DecTr = wu_DecTr*wu_ext/wu
1413 ! wu_Grass = wu_Grass*wu_ext/wu
1414 wu_surf = wu_surf*wu_ext/wu
1415 END IF
1416
real(kind(1d0)) function get_prof_spectime_sum(Hour, Min, Sec, Prof_24h, dt)

References allocatearray::conifsurf, allocatearray::decidsurf, and allocatearray::grasssurf.

Referenced by suews_driver::suews_cal_main().

Here is the caller graph for this function:

◆ suews_update_soilmoist()

subroutine waterdist_module::suews_update_soilmoist ( real(kind(1d0)), intent(in)  NonWaterFraction,
real(kind(1d0)), dimension(nsurf), intent(in)  SoilStoreCap,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  soilstore_id,
real(kind(1d0)), intent(out)  SoilMoistCap,
real(kind(1d0)), intent(out)  SoilState,
real(kind(1d0)), intent(out)  vsmd,
real(kind(1d0)), intent(out)  smd 
)

Definition at line 826 of file suews_phys_waterdist.f95.

831 IMPLICIT NONE
832
833 ! INTEGER,INTENT(in)::nsurf,ConifSurf,DecidSurf,GrassSurf
834 REAL(KIND(1D0)), INTENT(in) :: NonWaterFraction
835 REAL(KIND(1D0)), INTENT(in), DIMENSION(nsurf) :: SoilStoreCap, sfr_surf, soilstore_id
836
837 REAL(KIND(1D0)), INTENT(out) :: SoilMoistCap, SoilState
838 REAL(KIND(1D0)), INTENT(out) :: vsmd, smd
839
840 INTEGER :: is
841 REAL(KIND(1D0)) :: fveg
842
843 soilmoistcap = 0 !Maximum capacity of soil store [mm] for whole surface
844 soilstate = 0 !Area-averaged soil moisture [mm] for whole surface
845
846 IF (nonwaterfraction /= 0) THEN !Soil states only calculated if soil exists. LJ June 2017
847 DO is = 1, nsurf - 1 !No water body included
848 soilmoistcap = soilmoistcap + (soilstorecap(is)*sfr_surf(is)/nonwaterfraction)
849 soilstate = soilstate + (soilstore_id(is)*sfr_surf(is)/nonwaterfraction)
850 END DO
851 END IF
852
853 !If loop removed HCW 26 Feb 2015
854 !if (ir==1) then !Calculate initial smd
855 smd = soilmoistcap - soilstate
856 !endif
857
858 ! Calculate soil moisture for vegetated surfaces only (for use in surface conductance)
859 vsmd = 0
860 IF ((sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf)) > 0) THEN
861
862 fveg = sfr_surf(is)/(sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf))
863 ELSE
864 fveg = 0
865 END IF
866 DO is = conifsurf, grasssurf !Vegetated surfaces only
867 IF (fveg == 0) THEN
868 vsmd = 0
869 ELSE
870 vsmd = vsmd + (soilstorecap(is) - soilstore_id(is))*sfr_surf(is)/fveg
871 END IF
872 !write(*,*) is, vsmd, smd
873 END DO
874

References allocatearray::conifsurf, allocatearray::decidsurf, allocatearray::grasssurf, and allocatearray::nsurf.

Referenced by suews_driver::suews_cal_main().

Here is the caller graph for this function:

◆ updateflood()

subroutine waterdist_module::updateflood ( integer, intent(in)  is,
real(kind(1d0)), dimension(nsurf), intent(in)  runoff,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(in)  PipeCapacity,
real(kind(1d0)), intent(in)  RunoffToWater,
real(kind(1d0)), intent(inout)  runoffAGimpervious,
real(kind(1d0)), intent(inout)  surplusWaterBody,
real(kind(1d0)), intent(inout)  runoffAGveg,
real(kind(1d0)), intent(inout)  runoffPipes 
)

Definition at line 724 of file suews_phys_waterdist.f95.

728
729 IMPLICIT NONE
730
731 INTEGER, INTENT(in) :: is
732 REAL(KIND(1D0)), INTENT(in) :: sfr_surf(nsurf), runoff(nsurf), PipeCapacity, RunoffToWater
733 REAL(KIND(1D0)), INTENT(inout) :: runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes
734
735 ! Add runoff to pipes
736 runoffpipes = runoffpipes + (runoff(is)*sfr_surf(is))
737
738 ! If pipe capacity is full, surface runoff occurs
739 ! N.B. this will happen each loop (replicates pipes filling up)
740 IF (runoffpipes > pipecapacity) THEN
741
742 !------Paved and building surface
743 IF (is == pavsurf .OR. is == bldgsurf) THEN
744 IF (sfr_surf(watersurf) > 0.0000001) THEN
745 ! If there is some water present, the water surface will take some of the flood water (fraction RunoffToWater)
746 ! RunoffToWater is specified in SUEWS_SiteSelect.txt
747 runoffagimpervious = runoffagimpervious + (runoffpipes - pipecapacity)*(1 - runofftowater)
748 surpluswaterbody = surpluswaterbody + (runoffpipes - pipecapacity)*runofftowater
749 ELSE
750 ! Otherwise, all flood water must go to runoff
751 runoffagimpervious = runoffagimpervious + (runoffpipes - pipecapacity)
752 END IF
753 !------other surfaces
754 ELSEIF (is >= conifsurf .AND. is <= bsoilsurf) THEN
755 IF (sfr_surf(watersurf) > 0.0000001) THEN
756 ! If there is some water present, the water surface will take some of the flood water (fraction RunoffToWater)
757 runoffagveg = runoffagveg + (runoffpipes - pipecapacity)*(1 - runofftowater)
758 surpluswaterbody = surpluswaterbody + (runoffpipes - pipecapacity)*runofftowater
759 ELSE
760 ! Otherwise, all flood water must go to runoff
761 runoffagveg = runoffagveg + (runoffpipes - pipecapacity)
762 END IF
763 END IF
764
765 runoffpipes = pipecapacity !Pipes are at their max capacity
766
767 END IF !If runoff exceed pipe capacity
768

References allocatearray::bldgsurf, allocatearray::bsoilsurf, allocatearray::conifsurf, allocatearray::pavsurf, and allocatearray::watersurf.

Referenced by cal_water_storage(), and snow_module::snowcalc().

Here is the caller graph for this function: