SUEWS API Site
Documentation of SUEWS source code
suews_phys_waterdist.f95
Go to the documentation of this file.
2  USE allocatearray, ONLY: nsurf, &
3  pavsurf, bldgsurf, &
6 
7  IMPLICIT NONE
8  ! INTEGER, PARAMETER :: nsurf = 7
9  ! INTEGER, PARAMETER :: PavSurf = 1
10  ! INTEGER, PARAMETER :: BldgSurf = 2
11  ! INTEGER, PARAMETER :: ConifSurf = 3
12  ! INTEGER, PARAMETER :: DecidSurf = 4
13  ! INTEGER, PARAMETER :: GrassSurf = 5
14  ! INTEGER, PARAMETER :: BSoilSurf = 6
15  ! INTEGER, PARAMETER :: WaterSurf = 7
16  ! INTEGER, PARAMETER :: ExcessSurf = 8
17 CONTAINS
18 
19  !------------------------------------------------------------------------------
20  SUBROUTINE drainage( &
21  is, & !input
22  state_is, &
23  StorCap, &
24  DrainEq, &
25  DrainCoef1, &
26  DrainCoef2, &
27  nsh_real, &
28  drain_is)!output
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  ENDIF
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  ENDIF
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  ENDIF
84  ENDIF
85 
86  RETURN
87 
88  END SUBROUTINE drainage
89  !------------------------------------------------------------------------------
90 
91  !--------------Calculation of water storage change of specific land cover------------------------------
92  SUBROUTINE cal_water_storage( &
93  is, sfr, PipeCapacity, RunoffToWater, pin, & ! input:
94  WU_nsurf, &
95  drain, AddWater, addImpervious, nsh_real, stateOld, AddWaterRunoff, &
96  PervFraction, addVeg, SoilStoreCap, addWaterBody, FlowChange, StateLimit, &
97  runoffAGimpervious, surplusWaterBody, & ! inout:
98  runoffAGveg, runoffPipes, ev, soilstore_id, SurplusEvap, runoffWaterBody, &
99  p_mm, chang, runoff, state_id)!output:
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(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! 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)::stateOld!Wetness status of each surface type from previous timestep [mm]
148  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::AddWaterRunoff!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 
152  REAL(KIND(1d0)), INTENT(in)::PipeCapacity!Capacity of pipes to transfer water
153  REAL(KIND(1d0)), INTENT(in)::RunoffToWater!Fraction of surface runoff going to water body
154  REAL(KIND(1d0)), INTENT(in)::pin!Rain per time interval
155  REAL(KIND(1d0)), INTENT(in)::addImpervious!Water from impervious surfaces of other grids [mm] for whole surface area
156  REAL(KIND(1d0)), INTENT(in)::nsh_real!nsh cast as a real for use in calculations
157  REAL(KIND(1d0)), INTENT(in)::PervFraction! sum of surface cover fractions for impervious surfaces
158  REAL(KIND(1d0)), INTENT(in)::addVeg!Water from vegetated surfaces of other grids [mm] for whole surface area
159  REAL(KIND(1d0)), INTENT(in)::addWaterBody!Water from water surface of other grids [mm] for whole surface area
160  REAL(KIND(1d0)), INTENT(in)::FlowChange!Difference between the input and output flow in the water body
161 
162  REAL(KIND(1d0)), INTENT(inout)::runoffAGimpervious!Above ground runoff from impervious surface [mm] for whole surface area
163  REAL(KIND(1d0)), INTENT(inout)::surplusWaterBody!Extra runoff that goes to water body [mm] as specified by RunoffToWater
164  REAL(KIND(1d0)), INTENT(inout)::runoffAGveg!Above ground runoff from vegetated surfaces [mm] for whole surface area
165  REAL(KIND(1d0)), INTENT(inout)::runoffPipes!Runoff in pipes [mm] for whole surface area
166  REAL(KIND(1d0)), INTENT(inout)::ev!Evaporation
167  REAL(KIND(1d0)), INTENT(inout)::runoffWaterBody!Above ground runoff from water surface [mm] for whole surface area
168 
169  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout)::soilstore_id !Soil moisture of each surface type [mm]
170  REAL(KIND(1d0)), DIMENSION(2), INTENT(inout) ::SurplusEvap!Surplus for evaporation in 5 min timestep
171 
172  REAL(KIND(1d0)), INTENT(out)::p_mm!Inputs to surface water balance
173 
174  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::drain !Drainage of each surface type [mm]
175  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::WU_nsurf !external water use of each surface type [mm]
176  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::chang !Change in state_id [mm]
177  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::runoff!Runoff from each surface type [mm]
178  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::state_id !Wetness status of each surface type [mm]
179 
180  !Extra evaporation [mm] from impervious surfaces which cannot happen due to lack of water
181  REAL(KIND(1d0)):: EvPart
182  REAL(KIND(1d0)), PARAMETER:: NotUsed = -55.5
183  REAL(KIND(1d0)), PARAMETER:: IPThreshold_mmhr = 10 ! NB:this should be an input and can be specified. SG 25 Apr 2018
184 
185  !Initialise extra evaporation to zero
186  evpart = 0
187 
188  !Initialise runoff to zero
189  runoff(is) = 0
190 
191  !SurfaceFlood(is) = 0 !!This probably needs to be carried over between timesteps, but reset for now
192 
193  !==================================================================
194  ! Combine water inputs to the current surface
195  ! Add external water use for each surface type
196  ! SELECT CASE (is)
197  ! CASE (ConifSurf)
198  ! p_mm = pin + wu_EveTr
199  ! CASE (DecidSurf)
200  ! p_mm = pin + wu_DecTr
201  ! CASE (GrassSurf)
202  ! p_mm = pin + wu_Grass
203  ! CASE default
204  ! p_mm = pin
205  ! END SELECT
206 
207  ! Combine water inputs to the current surface
208  ! Add external water use for each surface type
209  p_mm = pin + wu_nsurf(is)
210 
211  ! Add water from other surfaces within the same grid (RS2S) ----
212  ! AddWater is the water supplied to the current surface from other surfaces
213  ! i.e. drain*WaterDist (see SUEWS_ReDistributeWater)
214  p_mm = p_mm + addwater(is)
215  !==================================================================
216 
217  !========surface-specific calculation=========================
218  SELECT CASE (is)
219  CASE (pavsurf, bldgsurf)
220  !==== Impervious surfaces (Paved, Buildings) ======================
221  ! Add water from neighbouring grids (RG2G)
222  ! Add to PavSurf only, as water cannot flow onto buildings
223  IF (is == pavsurf) THEN
224  IF (sfr(pavsurf) /= 0) THEN ! If loop added HCW 08 Dec 2015
225  p_mm = p_mm + addimpervious/sfr(pavsurf)
226  ENDIF
227  ENDIF
228 
229  ! Calculate change in surface state_id (inputs - outputs)
230  chang(is) = p_mm - (drain(is) + ev)
231 
232  ! If p_mm is too large, excess goes to runoff (i.e. the rate of water supply is too fast)
233  ! and does not affect state_id
234  IF (p_mm > ipthreshold_mmhr/nsh_real) THEN
235  runoff(is) = runoff(is) + (p_mm - ipthreshold_mmhr/nsh_real)
236  chang(is) = ipthreshold_mmhr/nsh_real - (drain(is) + ev)
237  ENDIF
238 
239  ! Calculate updated state_id using chang
240  state_id(is) = stateold(is) + chang(is)
241 
242  ! Check state_id is within physical limits between zero (dry) and max. storage capacity
243  IF (state_id(is) < 0.0) THEN ! Cannot have a negative surface state_id
244  ! If there is not sufficient water on the surface, then don't allow this evaporation to happen
245  ! Allow evaporation only until surface is dry (state_id(is)=0); additional evaporation -> evaporation surplus
246  surplusevap(is) = abs(state_id(is)) !Surplus evaporation is that which tries to evaporate non-existent water
247  ev = ev - surplusevap(is) !Limit evaporation according to water availability
248  state_id(is) = 0.0 !Now surface is dry
249  ! elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
250  ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
251  ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
252  ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
253  ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
254  ENDIF
255 
256  ! Recalculate change in surface state_id from difference with previous timestep
257  chang(is) = state_id(is) - stateold(is)
258 
259  ! Runoff -------------------------------------------------------
260  ! For impervious surfaces, some of drain(is) becomes runoff
261  runoff(is) = runoff(is) + drain(is)*addwaterrunoff(is) !Drainage (that is not flowing to other surfaces) goes to runoff
262 
263  !So, up to this point, runoff(is) can have contributions if
264  ! p_mm > ipthreshold (water input too fast)
265  ! state_id > StoreDrainPrm(6,is) (net water exceeds storage capacity)
266  ! WaterDist specifies some fraction of drain(is) -> runoff
267 
268  CASE (conifsurf:bsoilsurf)
269  !==== For Conif, Decid, Grass, BSoil surfaces ==================
270  ! Transfer evaporation surplus from impervious surfaces to pervious surfaces
271  evpart = merge( &
272  dot_product(surplusevap(pavsurf:bldgsurf), sfr(pavsurf:bldgsurf)/pervfraction), &
273  0d0, &
274  pervfraction /= 0)
275 
276  ! Add surplus evaporation to ev for pervious surfaces
277  ev = ev + evpart
278 
279  ! ---- Add water from neighbouring grids (RG2G) ----
280  ! Add to Grass and BSoil only, as water cannot flow onto trees
281  IF (is == grasssurf .OR. is == bsoilsurf) THEN
282  IF ((sfr(grasssurf) + sfr(bsoilsurf)) /= 0) THEN
283  p_mm = p_mm + addveg/(sfr(grasssurf) + sfr(bsoilsurf))
284  ENDIF
285  ENDIF
286 
287  ! Calculate change in surface state_id (inputs - outputs)
288  chang(is) = p_mm - (drain(is) + ev)
289 
290  ! If p_mm is too large, excess goes to runoff (i.e. the rate of water supply is too fast)
291  ! and does not affect state_id
292  IF (p_mm > ipthreshold_mmhr/nsh_real) THEN
293  runoff(is) = runoff(is) + (p_mm - ipthreshold_mmhr/nsh_real)
294  chang(is) = ipthreshold_mmhr/nsh_real - (drain(is) + ev)
295  ENDIF
296 
297  ! Calculate updated state_id using chang
298  state_id(is) = stateold(is) + chang(is)
299 
300  ! Check state_id is within physical limits between zero (dry) and max. storage capacity
301  IF (state_id(is) < 0.0) THEN ! Cannot have a negative surface state_id
302  ! If there is not sufficient water on the surface, then remove water from soilstore
303  ! Allow evaporation until soilstore_id is depleted and surface is dry
304  IF ((soilstore_id(is) + state_id(is)) >= 0) THEN
305  soilstore_id(is) = soilstore_id(is) + state_id(is)
306  state_id(is) = 0.0
307  ! If there is not sufficient water on the surface or soilstore, then don't allow this evaporation to happen
308  ELSE
309  ev = ev - abs(state_id(is)) !Limit evaporation according to water availability
310  state_id(is) = 0.0 !Now surface is dry
311  ENDIF
312 
313  !elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
314  ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
315  ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
316  ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
317  ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
318  ENDIF
319 
320  ! Recalculate change in surface state_id from difference with previous timestep
321  chang(is) = state_id(is) - stateold(is)
322 
323  !Where should this go? Used to be before previous part!!
324  ! soilstore_id -------------------------------------------------
325  ! For pervious surfaces (not water), some of drain(is) goes to soil storage
326  ! Drainage (that is not flowing to other surfaces) goes to soil storages
327  soilstore_id(is) = soilstore_id(is) + drain(is)*addwaterrunoff(is)
328 
329  ! If soilstore is full, the excess will go to runoff
330  IF (soilstore_id(is) > soilstorecap(is)) THEN ! TODO: this should also go to flooding of some sort
331  runoff(is) = runoff(is) + (soilstore_id(is) - soilstorecap(is))
332  soilstore_id(is) = soilstorecap(is)
333  ELSEIF (soilstore_id(is) < 0) THEN !! QUESTION: But where does this lack of water go? !!Can this really happen here?
334  CALL errorhint(62, 'SUEWS_store: soilstore_id(is) < 0 ', soilstore_id(is), notused, is)
335  ! Code this properly - soilstore_id(is) < 0 shouldn't happen given the above loops
336  !soilstore_id(is)=0 !Groundwater / deeper soil should kick in
337  ENDIF
338 
339  CASE (watersurf)
340  IF (sfr(watersurf) /= 0) THEN
341 
342  ! ---- Add water from neighbouring grids (RG2G) ----
343  p_mm = p_mm + addwaterbody/sfr(watersurf)
344 
345  ! Calculate change in surface state_id (inputs - outputs)
346  ! No drainage for water surface
347  ! FlowChange is the difference in input and output flows [mm hr-1]
348  chang(is) = p_mm + flowchange/nsh_real - ev
349 
350  ! Calculate updated state_id using chang
351  state_id(is) = stateold(is) + chang(is)
352 
353  ! Check state_id is within physical limits between zero (dry) and max. storage capacity
354  IF (state_id(is) < 0.0) THEN ! Cannot have a negative surface state_id
355  ! If there is not sufficient water on the surface, then don't allow this evaporation to happen
356  ev = ev - abs(state_id(is)) !Limit evaporation according to water availability
357  state_id(is) = 0.0 !Now surface is dry
358  !elseif (state_id(is)>StoreDrainPrm(6,is)) then !!This should perhaps be StateLimit(is)
359  ! !! If state_id exceeds the storage capacity, then the excess goes to surface flooding
360  ! !SurfaceFlood(is)=SurfaceFlood(is)+(state_id(is)-StoreDrainPrm(6,is)) !!Need to deal with this properly
361  ! runoff(is)=runoff(is)+(state_id(is)-StoreDrainPrm(6,is)) !!needs to go to flooding
362  ! state_id(is)=StoreDrainPrm(6,is) !Now surface state_id is at max (storage) capacity
363  ENDIF
364 
365  ! Recalculate change in surface state_id from difference with previous timestep
366  chang(is) = state_id(is) - stateold(is)
367 
368  ! If state_id exceeds limit, then excess goes to runoff (currently applies to water StoreDrainPrm only)
369  IF (state_id(watersurf) > statelimit(watersurf)) THEN
370  runoff(watersurf) = runoff(watersurf) + (state_id(watersurf) - statelimit(watersurf))
371  state_id(watersurf) = statelimit(watersurf)
372  runoffwaterbody = runoffwaterbody + runoff(watersurf)*sfr(watersurf)
373  ELSE
374  state_id(watersurf) = state_id(watersurf) + surpluswaterbody
375  IF (state_id(watersurf) > statelimit(watersurf)) THEN
376  runoffwaterbody = runoffwaterbody + (state_id(watersurf) - statelimit(watersurf))*sfr(watersurf)
377  state_id(watersurf) = statelimit(watersurf)
378  ENDIF
379  ENDIF
380 
381  ! Recalculate change in surface state_id from difference with previous timestep
382  chang(is) = state_id(is) - stateold(is)
383  ENDIF
384  END SELECT
385  !==================================================================
386 
387  !==== RUNOFF ======================================================
388  ! TODO: to consider areas here - SurfaceArea may vary between grids too
389  ! - also implement where water for next surface is calculated (RunoffFromGrid subroutine)
390  ! Calculations of the piperunoff exceedensances moved to separate subroutine so that from snow same
391  ! calculations can be made. LJ in May 2015
392 
393  IF (is < watersurf) THEN !Not for water body
394  ! CALL updateFlood
395  CALL updateflood( &
396  is, runoff, &! input:
397  sfr, pipecapacity, runofftowater, &
398  runoffagimpervious, surpluswaterbody, runoffagveg, runoffpipes)! inout:
399  ENDIF
400 
401  END SUBROUTINE cal_water_storage
402  !------------------------------------------------------------------------------
403 
404  !------------------------------------------------------------------------------
405  SUBROUTINE updateflood( &
406  is, runoff, &! input:
407  sfr, PipeCapacity, RunoffToWater, &
408  runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes)! inout:
410  IMPLICIT NONE
411 
412  INTEGER, INTENT(in) :: is
413  REAL(KIND(1d0)), INTENT(in) :: sfr(nsurf), runoff(nsurf), PipeCapacity, RunoffToWater
414  REAL(KIND(1d0)), INTENT(inout) :: runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes
415 
416  ! Add runoff to pipes
417  runoffpipes = runoffpipes + (runoff(is)*sfr(is))
418 
419  ! If pipe capacity is full, surface runoff occurs
420  ! N.B. this will happen each loop (replicates pipes filling up)
421  IF (runoffpipes > pipecapacity) THEN
422 
423  !------Paved and building surface
424  IF (is == pavsurf .OR. is == bldgsurf) THEN
425  IF (sfr(watersurf) > 0.0000001) THEN
426  ! If there is some water present, the water surface will take some of the flood water (fraction RunoffToWater)
427  ! RunoffToWater is specified in SUEWS_SiteSelect.txt
428  runoffagimpervious = runoffagimpervious + (runoffpipes - pipecapacity)*(1 - runofftowater)
429  surpluswaterbody = surpluswaterbody + (runoffpipes - pipecapacity)*runofftowater
430  ELSE
431  ! Otherwise, all flood water must go to runoff
432  runoffagimpervious = runoffagimpervious + (runoffpipes - pipecapacity)
433  ENDIF
434  !------other surfaces
435  ELSEIF (is >= conifsurf .AND. is <= bsoilsurf) THEN
436  IF (sfr(watersurf) > 0.0000001) THEN
437  ! If there is some water present, the water surface will take some of the flood water (fraction RunoffToWater)
438  runoffagveg = runoffagveg + (runoffpipes - pipecapacity)*(1 - runofftowater)
439  surpluswaterbody = surpluswaterbody + (runoffpipes - pipecapacity)*runofftowater
440  ELSE
441  ! Otherwise, all flood water must go to runoff
442  runoffagveg = runoffagveg + (runoffpipes - pipecapacity)
443  ENDIF
444  ENDIF
445 
446  runoffpipes = pipecapacity !Pipes are at their max capacity
447 
448  ENDIF !If runoff exceed pipe capacity
449 
450  END SUBROUTINE updateflood
451  !------------------------------------------------------------------------------
452 
453  !------------------------------------------------------------------------------
454  SUBROUTINE redistributewater( &
455  snowUse, WaterDist, sfr, Drain, &! input:
456  AddWaterRunoff, AddWater)! output:
457  !Drainage moves into different parts defined by WaterDistSS_YYYY.txt. LJ 2010
458  !AddWater(is) is that amount of water that is gained for each surface
459  !Latest update takes snow into account. 22/03/2013 LJ
460  !-------------------------------------------------------------------
461 
462  IMPLICIT NONE
463  INTEGER, INTENT(in)::snowUse!Snow part used (1) or not used (0)
464 
465  REAL(KIND(1d0)), INTENT(in)::WaterDist(nsurf + 1, nsurf - 1) !Within-grid water distribution to other surfaces and runoff/soil store [-]
466  REAL(KIND(1d0)), INTENT(in)::sfr(nsurf) !Surface fractions [-]
467  REAL(KIND(1d0)), INTENT(in)::Drain(nsurf) !Drainage of each surface type [mm]
468 
469  REAL(KIND(1d0)), INTENT(out)::AddWaterRunoff(nsurf)!Fraction of water going to runoff/sub-surface soil (WGWaterDist) [-]
470  REAL(KIND(1d0)), INTENT(out)::AddWater(nsurf) !Water from other surfaces (WGWaterDist in SUEWS_ReDistributeWater.f95) [mm]
471 
472  INTEGER::ii, jj
473  INTEGER::NSurfDoNotReceiveDrainage = 0!Number of surfaces that do not receive drainage water (green roof)
474 
475  !Fractions that go to runoff from each surface
476  DO ii = 1, nsurf - 1 !not water in the calculation
477  addwaterrunoff(ii) = waterdist(8, ii)
478  ENDDO
479  addwaterrunoff(watersurf) = 0
480  addwater = 0
481 
482  DO ii = 1, nsurf - nsurfdonotreceivedrainage !go through surfaces from 1 to 7. These gain water through drainage
483  DO jj = 1, nsurf - (nsurfdonotreceivedrainage + 1) !From where surface ii can gain water - can't gain water from itself
484 
485  IF (sfr(ii) /= 0) THEN !Water movement takes place only if surface fraction exists
486 
487  !No snow calculations!
488  IF (snowuse == 0) THEN
489  addwater(ii) = addwater(ii) + (drain(jj)*sfr(jj)/sfr(ii))*waterdist(ii, jj) !Original
490 
491  !Snow included, This needs to be fixed at some point. LJ Mar 2013
492  ELSE
493  addwaterrunoff(jj) = addwaterrunoff(jj) + waterdist(ii, jj) !No receiving surface -> runoff
494  ENDIF
495 
496  ELSE
497  addwaterrunoff(jj) = addwaterrunoff(jj) + waterdist(ii, jj) !If no receiving surface exists,
498  !water fraction goes to AddWaterRunoff
499  ENDIF
500  ENDDO
501  ENDDO
502 
503  END SUBROUTINE redistributewater
504  !------------------------------------------------------------------------------
505 
506  !------------------------------------------------------------------------------
507  SUBROUTINE suews_update_soilmoist( &
508  NonWaterFraction, &!input
509  SoilStoreCap, sfr, soilstore_id, &
510  SoilMoistCap, SoilState, &!output
511  vsmd, smd)
512  IMPLICIT NONE
513 
514  ! INTEGER,INTENT(in)::nsurf,ConifSurf,DecidSurf,GrassSurf
515  REAL(KIND(1d0)), INTENT(in)::NonWaterFraction
516  REAL(KIND(1d0)), INTENT(in), DIMENSION(nsurf)::SoilStoreCap, sfr, soilstore_id
517 
518  REAL(KIND(1d0)), INTENT(out)::SoilMoistCap, SoilState
519  REAL(KIND(1d0)), INTENT(out)::vsmd, smd
520 
521  INTEGER :: is
522 
523  soilmoistcap = 0 !Maximum capacity of soil store [mm] for whole surface
524  soilstate = 0 !Area-averaged soil moisture [mm] for whole surface
525 
526  IF (nonwaterfraction /= 0) THEN !Soil states only calculated if soil exists. LJ June 2017
527  DO is = 1, nsurf - 1 !No water body included
528  soilmoistcap = soilmoistcap + (soilstorecap(is)*sfr(is)/nonwaterfraction)
529  soilstate = soilstate + (soilstore_id(is)*sfr(is)/nonwaterfraction)
530  ENDDO
531  ENDIF
532 
533  !If loop removed HCW 26 Feb 2015
534  !if (ir==1) then !Calculate initial smd
535  smd = soilmoistcap - soilstate
536  !endif
537 
538  ! Calculate soil moisture for vegetated surfaces only (for use in surface conductance)
539  vsmd = 0
540  DO is = conifsurf, grasssurf !Vegetated surfaces only
541  IF (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) == 0) THEN
542  vsmd = 0
543  ELSE
544  vsmd = vsmd + (soilstorecap(is) - soilstore_id(is))*sfr(is)/(sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf))
545  END IF
546  !write(*,*) is, vsmd, smd
547  ENDDO
548 
549  END SUBROUTINE suews_update_soilmoist
550  !------------------------------------------------------------------------------
551 
552  !========== Calculate soil moisture of a whole grid ============
553  SUBROUTINE suews_cal_soilstate( &
554  SMDMethod, xsmd, NonWaterFraction, SoilMoistCap, &!input
555  SoilStoreCap, surf_chang_per_tstep, &
556  soilstore_id, soilstoreOld, sfr, &
557  smd, smd_nsurf, tot_chang_per_tstep, SoilState)!output
559  IMPLICIT NONE
560  ! INTEGER, PARAMETER :: nsurf = 7
561 
562  INTEGER, INTENT(in) ::SMDMethod
563  REAL(KIND(1d0)), INTENT(in)::xsmd
564  REAL(KIND(1d0)), INTENT(in)::NonWaterFraction
565  REAL(KIND(1d0)), INTENT(in)::SoilMoistCap
566 
567  REAL(KIND(1d0)), INTENT(in)::surf_chang_per_tstep
568  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::soilstore_id !Soil moisture of each surface type [mm]
569  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::soilstoreOld !Soil moisture of each surface type from previous timestep [mm]
570  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
571  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SoilStoreCap !Capacity of soil store for each surface [mm]
572 
573  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::smd_nsurf !smd for each surface
574  REAL(KIND(1d0)), INTENT(out)::SoilState !Area-averaged soil moisture [mm] for whole surface
575  REAL(KIND(1d0)), INTENT(out)::smd !One value for whole surface
576  REAL(KIND(1d0)), INTENT(out)::tot_chang_per_tstep !Change in surface state_id
577 
578  REAL(KIND(1d0)), PARAMETER::NotUsed = -999
579  REAL(KIND(1d0)), PARAMETER::NAN = -999
580  INTEGER :: is
581 
582  soilstate = 0 !Area-averaged soil moisture [mm] for whole surface
583  IF (nonwaterfraction /= 0) THEN !Fixed for water surfaces only
584  DO is = 1, nsurf - 1 !No water body included
585  soilstate = soilstate + (soilstore_id(is)*sfr(is)/nonwaterfraction)
586  IF (soilstate < 0) THEN
587  CALL errorhint(62, 'SUEWS_Calculations: total SoilState < 0 (just added surface is) ', soilstate, notused, is)
588  ELSEIF (soilstate > soilmoistcap) THEN
589  CALL errorhint(62, 'SUEWS_Calculations: total SoilState > capacity (just added surface is) ', soilstate, notused, is)
590  !SoilMoist_state=SoilMoistCap !What is this LJ 10/2010 - QUESTION: SM exceeds capacity, but where does extra go?HCW 11/2014
591  ENDIF
592  ENDDO !end loop over surfaces
593  ! SoilState = DOT_PRODUCT(soilstore_id(1:nsurf - 1), sfr(1:nsurf - 1))/NonWaterFraction
594  ! IF (SoilState < 0) THEN
595  ! CALL ErrorHint(62, 'SUEWS_Calculations: total SoilState < 0 (just added surface is) ', SoilState, NotUsed, is)
596  ! ELSEIF (SoilState > SoilMoistCap) THEN
597  ! CALL ErrorHint(62, 'SUEWS_Calculations: total SoilState > capacity (just added surface is) ', SoilState, NotUsed, is)
598  ! !SoilMoist_state=SoilMoistCap !What is this LJ 10/2010 - QUESTION: SM exceeds capacity, but where does extra go?HCW 11/2014
599  ! ENDIF
600  ENDIF
601 
602  ! Calculate soil moisture deficit
603  smd = soilmoistcap - soilstate !One value for whole surface
604  smd_nsurf = soilstorecap - soilstore_id !smd for each surface
605 
606  ! Soil stores can change after horizontal water movements
607  ! Calculate total change in surface and soil state_id
608  tot_chang_per_tstep = surf_chang_per_tstep !Change in surface state_id
609  DO is = 1, (nsurf - 1) !No soil for water surface (so change in soil moisture is zero)
610  tot_chang_per_tstep = tot_chang_per_tstep + ((soilstore_id(is) - soilstoreold(is))*sfr(is)) !Add change in soil state_id
611  ENDDO
612 
613  IF (smdmethod > 0) THEN ! use observed value
614  ! smd_nsurf=NAN
615  smd_nsurf = nan
616  smd = xsmd
617  ENDIF
618 
619  END SUBROUTINE suews_cal_soilstate
620  !===================================================================================
621 
622  SUBROUTINE suews_cal_horizontalsoilwater( &
623  sfr, &! input: ! surface fractions
624  SoilStoreCap, &!Capacity of soil store for each surface [mm]
625  SoilDepth, &!Depth of sub-surface soil store for each surface [mm]
626  SatHydraulicConduct, &!Saturated hydraulic conductivity for each soil subsurface [mm s-1]
627  SurfaceArea, &!Surface area of the study area [m2]
628  NonWaterFraction, &! sum of surface cover fractions for all except water surfaces
629  tstep_real, & !tstep cast as a real for use in calculations
630  soilstore_id, &! inout: !Soil moisture of each surface type [mm]
631  runoffSoil, &!Soil runoff from each soil sub-surface [mm]
632  runoffSoil_per_tstep &! output:!Runoff to deep soil per timestep [mm] (for whole surface, excluding water body)
633  )
634  !Transfers water in soil stores of land surfaces LJ (2010)
635  !Change the model to use varying hydraulic conductivity instead of constant value LJ (7/2011)
636  !If one of the surface's soildepth is zero, no water movement is considered
637  ! LJ 15/06/2017 Modification: - Moved location of runoffSoil_per_tstep within previous if-loop to avoid dividing with zero with 100% water surface
638  ! HCW 22/02/2017 Modifications: - Minor bug fixed in VWC1/B_r1 comparison - if statements reversed
639  ! HCW 13/08/2014 Modifications: - Order of surfaces reversed (for both is and jj loops)
640  ! - Number of units (e.g. properties) added to distance calculation
641  ! HCW 12/08/2014 Modifications: - Distance changed from m to mm in dI_dt calculation
642  ! - dI_dt [mm s-1] multiplied by no. seconds in timestep -> dI [mm]
643  ! - if MatPot is set to max. value (100000 mm), Km set to 0 mm s-1
644  ! - Provide parameters for residual volumetric soil moisture [m3 m-3]
645  ! (currently hard coded as 0.1 m3 m-3 for testing)
646  !
647  !------------------------------------------------------
648  ! use SUES_data
649  ! use gis_data
650  ! use time
651  ! use allocateArray
652 
653  IMPLICIT NONE
654 
655  REAL(KIND(1d0)), INTENT(in) ::sfr(nsurf)! surface fractions
656  REAL(KIND(1d0)), INTENT(in) ::SoilStoreCap(nsurf)!Capacity of soil store for each surface [mm]
657  REAL(KIND(1d0)), INTENT(in) ::SoilDepth(nsurf)!Depth of sub-surface soil store for each surface [mm]
658  REAL(KIND(1d0)), INTENT(in) ::SatHydraulicConduct(nsurf)!Saturated hydraulic conductivity for each soil subsurface [mm s-1]
659  REAL(KIND(1d0)), INTENT(in) ::SurfaceArea!Surface area of the study area [m2]
660  REAL(KIND(1d0)), INTENT(in) ::NonWaterFraction! sum of surface cover fractions for all except water surfaces
661  REAL(KIND(1d0)), INTENT(in) ::tstep_real !tstep cast as a real for use in calculations
662 
663  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout) ::soilstore_id!Soil moisture of each surface type [mm]
664  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(inout) ::runoffSoil!Soil runoff from each soil sub-surface [mm]
665 
666  REAL(KIND(1d0)), INTENT(out) :: runoffSoil_per_tstep!Runoff to deep soil per timestep [mm] (for whole surface, excluding water body)
667 
668  INTEGER::jj, is
669  REAL(KIND(1d0)):: &
670  DimenWaterCon1, DimenWaterCon2, &
671  SoilMoistCap_Vol1, &
672  SoilMoist_vol1, &
673  SoilMoistCap_Vol2, &
674  SoilMoist_vol2, &
675  B_r1, MatPot1, Km1, &
676  B_r2, MatPot2, Km2, &
677  Distance, KmWeight, dI, &
678  dI_dt!Water flow between two stores
679 
680  REAL(KIND(1d0)), PARAMETER:: &
681  alphavG = 0.0005, & !Set alphavG to match value in van Genuchten (1980) [mm-1]
682  nunits = 1 !Can change to represent plot/base unit size
683 
684  ! SoilMoist_vol1,2 = Volumetric soil moisture [m3 m-3]
685  ! SoilMoistCap_vol1,2 = Volumetric soil moisture capacity [m3 m-3] (from FunctionalTypes)
686  ! MatPot1,2 = Water potential (i.e. pressure head) of store [mm]
687  ! DimenWaterCon1,2 = Dimensionless water content, or relative saturation [-]
688  ! Distance = Distance between two stores [m]
689  ! B_r1,2 = Residual volumetric soil moisture [m3 m-3]
690  ! Km1,2 = Hydraulic conductivity of store [mm s-1]
691  ! KmWeight = Weighted hydraulic conductivity [mm s-1]
692  ! alphavG = Parameter (could depend on soil texture) [mm-1]
693  ! dI = Water flow between stores [mm] dI = dI_dt * no. secs in each timestep
694  ! if dI > 0, first surface gains water, second surface loses water
695  ! NUnits = Number of repeating units (e.g. properties, blocks) for distance calculation [-]
696 
697  runoffsoil_per_tstep = 0
698 
699  DO is = 1, nsurf - 1 !nsurf-1,1,-1 !Loop through each surface, excluding water surface (runs backwards as of 13/08/2014, HCW)
700 
701  IF (sfr(is) /= 0 .AND. soilstorecap(is) > 0) THEN !If particular surface area exists
702  ! and is capable of storing water (SoilStoreCap [mm])
703  DO jj = is + 1, nsurf - 1 !is-1,1,-1 !Sub-loop through remaining surfaces (runs backwards as of 13/08/2014, HCW)
704 
705  IF (sfr(jj) /= 0 .AND. soilstorecap(jj) > 0) THEN !If other surface area exists
706  ! and is capable of storing water
707 
708  ! ---- For surface 1 -----------------------------------------------------
709  ! Calculate non-saturated VWC
710  soilmoistcap_vol1 = soilstorecap(is)/soildepth(is) !Volumetric soil moisture capacity [m3 m-3] (i.e. saturated VWC)
711  soilmoist_vol1 = soilstore_id(is)/soildepth(is) !Volumetric soil moisture [m3 m-3]
712 
713  !B_r1=SoilMoistCap_Vol1-SoilMoist_vol1 !Residual soil moisture content [m3 m-3]
714  b_r1 = 0.1 !HCW 12/08/2014 Temporary fix
715  ! Need to add residual soil moisture values to FunctionalTypes
716  !B_r1=VolSoilMoistRes(is) !Residual soil moisture content [m3 m-3]
717 
718  !Order of if statements reversed HCW 22 Feb 2017
719  !If soil moisture less than or equal to residual value, set MatPot to max and Km to 0 to suppress water movement
720  IF (b_r1 >= soilmoist_vol1) THEN
721  matpot1 = 100000
722  km1 = 0 !Added by LJ in Nov 2013
723  ! Otherwise, there should be enough water in the soil to allow horizontal transfer
724  ELSE
725  dimenwatercon1 = (soilmoist_vol1 - b_r1)/(soilmoistcap_vol1 - b_r1) !Dimensionless water content [-]
726 
727  ! If very large or very small, adjust for calculation of MatPot and Km
728  IF (dimenwatercon1 > 0.99999) THEN
729  dimenwatercon1 = dimenwatercon1 - 0.0001 !This cannot equal 1
730  ENDIF
731 
732  IF (dimenwatercon1 < 0.00000005) THEN
733  dimenwatercon1 = dimenwatercon1 + 0.0000001 !Added HCW 22 Feb 2017
734  ENDIF
735 
736  !van Genuchten (1980), with n=2 and m = 1-1/n = 1/2
737  !Water potential of first store [mm] (van Genuchten 1980, Eq 3 rearranged)
738  matpot1 = sqrt(1/dimenwatercon1**2 - 1)/alphavg
739 
740  !Hydraulic conductivity of first store [mm s-1] (van Genuchten 1980, Eq 8)
741  km1 = sathydraulicconduct(is)*sqrt(dimenwatercon1)*(1 - (1 - dimenwatercon1**2)**0.5)**2
742 
743  !Check this value (HCW 12/08/2014)
744  IF (matpot1 > 100000) THEN
745  matpot1 = 100000 !Max. potential is 100000 mm (van Genuchten 1980)
746  km1 = 0 !Added by HCW 12/08/2014
747  ENDIF
748 
749  ENDIF
750 
751  ! ---- For surface 2 -----------------------------------------------------
752  ! Calculate non-saturated VWC
753  soilmoistcap_vol2 = soilstorecap(jj)/soildepth(jj) !Volumetric soil moisture capacity [m3 m-3] (i.e. saturated VWC)
754  soilmoist_vol2 = soilstore_id(jj)/soildepth(jj) !Volumetric soil moisture [m3 m-3]
755 
756  !B_r2=SoilMoistCap_Vol2-SoilMoist_vol2 !Residual soil moisture content [m3 m-3]
757  b_r2 = 0.1 !HCW 12/08/2014 Temporary fix
758  ! Need to add residual soil moisture values to FunctionalTypes
759  !B_r2=VolSoilMoistRes(jj) !Residual soil moisture content [m3 m-3]
760 
761  !If soil moisture below residual value, set MatPot to maximum
762  IF (b_r2 >= soilmoist_vol2) THEN
763  matpot2 = 100000
764  km2 = 0 !Added by LJ in Nov 2013
765  ELSE
766  dimenwatercon2 = (soilmoist_vol2 - b_r2)/(soilmoistcap_vol2 - b_r2) !Dimensionless water content [-]
767 
768  IF (dimenwatercon2 > 0.99999) THEN
769  dimenwatercon2 = dimenwatercon2 - 0.0001 !This cannot equal 1
770  ENDIF
771 
772  IF (dimenwatercon2 < 0.00000005) THEN
773  dimenwatercon2 = dimenwatercon2 + 0.0000001 !Added HCW 22 Feb 2017
774  ENDIF
775 
776  !van Genuchten (1980), with n=2 and m = 1-1/n = 1/2
777  !Water potential of second store [mm] (van Genuchten 1980, Eq 3 rearranged)
778  matpot2 = sqrt(1/dimenwatercon2**2 - 1)/alphavg
779 
780  !Hydraulic conductivity of second store [mm s-1] (van Genuchten 1980, Eq 8)
781  km2 = sathydraulicconduct(jj)*sqrt(dimenwatercon2)*(1 - (1 - dimenwatercon2**2)**0.5)**2
782 
783  IF ((matpot2) > 100000) THEN
784  matpot2 = 100000 !Max. potential is 100000 mm (van Genuchten 1980)
785  km2 = 0 !Added by HCW 12/08/2014
786  ENDIF
787 
788  ENDIF
789 
790  ! ------------------------------------------------------------------------
791 
792  !Find distance between the two stores (see Jarvi et al. 2011)
793  !SurfaceArea in m2 (changed from ha to m2 n SUEWS_Initial), so Distance in m
794  distance = (sqrt(sfr(is)*surfacearea/nunits) + sqrt(sfr(jj)*surfacearea/nunits))/2
795 
796  !Calculate areally-weighted hydraulic conductivity [mm s-1]
797  kmweight = (sfr(is)*km1 + sfr(jj)*km2)/(sfr(is) + sfr(jj))
798 
799  !Find water flow between the two stores [mm s-1] (Green-Ampt equation, Hillel 1971)
800  !Multiply Distance by 1000 to convert m to mm (HCW 12/08/2014)
801  di_dt = -(kmweight)*(-matpot1 + matpot2)/(distance*1000)
802 
803  !Multiply dI_dt by number of seconds in timestep to convert mm s-1 to mm
804  !Added by HCW 12/08/2014
805  di = di_dt*tstep_real !Use dI instead of dI_dt in the following calculations
806 
807  !Move water (in mm) ------------------------------------------------------
808  !Water moves only if (i) there is sufficient water to move and (ii) there is space to move it
809 
810  ! If there is sufficient water in both surfaces, allow movement of dI to occur
811  IF ((soilstore_id(jj) >= di*sfr(is)/sfr(jj)) .AND. ((soilstore_id(is) + di) >= 0)) THEN
812  soilstore_id(is) = soilstore_id(is) + di
813  soilstore_id(jj) = soilstore_id(jj) - di*sfr(is)/sfr(jj) !Check (HCW 13/08/2014) - QUESTION: why adjust for jj and not is?
814 
815  ! If insufficient water in first surface to move dI, instead move as much as possible
816  ELSEIF ((soilstore_id(is) + di) < 0) THEN
817  soilstore_id(jj) = soilstore_id(jj) + soilstore_id(is)*sfr(is)/sfr(jj) !HCW 12/08/2014 switched order of these two lines
818  soilstore_id(is) = 0 !Check (HCW 13/08/2014) - QUESTION: can SM actually go to zero, or is this inconsistent with SMres?
819 
820  ! If insufficient water in second surface to move dI, instead move as much as possible
821  ELSE
822  soilstore_id(is) = soilstore_id(is) + soilstore_id(jj)*sfr(jj)/sfr(is)
823  soilstore_id(jj) = 0
824  ENDIF
825 
826  !If soil moisture exceeds capacity, excess goes to soil runoff (first surface)
827  IF (soilstore_id(is) > soilstorecap(is)) THEN
828  runoffsoil(is) = runoffsoil(is) + (soilstore_id(is) - soilstorecap(is))
829  soilstore_id(is) = soilstorecap(is)
830  !elseif (soilstore_id(is)<0) then !HCW 13/08/2014 commented out as should never be true here anyway...
831  ! soilstore_id(is)=0 ! ... and if so, need to do more here (i.e. account for other water too)
832  ENDIF
833 
834  !If soil moisture exceeds capacity, excess goes to soil runoff (second surface)
835  IF (soilstore_id(jj) > soilstorecap(jj)) THEN
836  runoffsoil(jj) = runoffsoil(jj) + (soilstore_id(jj) - soilstorecap(jj))
837  soilstore_id(jj) = soilstorecap(jj)
838  !elseif (soilstore_id(jj)<0) then !HCW 13/08/2014 commented out (as above)
839  ! soilstore_id(jj)=0
840  ENDIF
841 
842  ENDIF !end if second surface exists and is capable of storing water
843 
844  ENDDO !end jj loop over second surface
845 
846  runoffsoil_per_tstep = runoffsoil_per_tstep + (runoffsoil(is)*sfr(is)/nonwaterfraction) !Excludes water body. Moved here as otherwise code crashed when NonWaterFraction=0
847 
848  ENDIF !end if first surface exists and is capable of storing water
849 
850  !runoffSoil_per_tstep=runoffSoil_per_tstep+(runoffSoil(is)*sfr(is)/NonWaterFraction) !Excludes water body
851 
852  ENDDO !is loop over first surface
853 
854  END SUBROUTINE suews_cal_horizontalsoilwater
855  !===================================================================================
856 
857  !===================================================================================
858  SUBROUTINE suews_cal_wateruse( &
859  nsh_real, & ! input:
860  wu_m3, SurfaceArea, sfr, &
861  IrrFracPaved, IrrFracBldgs, &
862  IrrFracEveTr, IrrFracDecTr, IrrFracGrass, &
863  IrrFracBSoil, IrrFracWater, &
864  DayofWeek_id, WUProfA_24hr, WUProfM_24hr, &
865  InternalWaterUse_h, HDD_id, WUDay_id, &
866  WaterUseMethod, NSH, it, imin, DLS, &
867  wu_nsurf, wu_int, wu_ext)! output:
868  ! Conversion of water use (irrigation)
869  ! Last modified:
870  ! TS 30 Nov 2019 - allow external water use for all surfaces
871  ! (previously only on vegetated surfaces)
872  ! TS 30 Oct 2018 - fixed a bug in external water use
873  ! TS 08 Aug 2017 - addded explicit interface
874  ! LJ 6 Apr 2017 - WUchoice changed to WaterUseMethod
875  ! TK 14 Mar 2017 - Corrected the variable name WUAreaEveTr_m2 -> WUAreaGrass_m2 (row 35)
876  ! Corrected conversion from m to mm /1000 -> *1000 (row 47 and 60)
877  ! LJ 27 Jan 2016 - Removing Tab:s and cleaning the code
878  ! HCW 12 Feb 2015 - Water use [mm] now inidcates the amount of water supplied for each surface
879  ! HCW 26 Jan 2015 - Water use [mm] is the same for each surface at the moment and indicates the
880  ! amount of water supplied for each irrigated area
881  !
882  ! To Do:
883  ! - Add functionality for water on paved surfaces (street cleaning, fountains)
884 
885  IMPLICIT NONE
886  ! INTEGER, PARAMETER :: nsurf = 7
887 
888  REAL(KIND(1d0)), INTENT(in)::nsh_real
889  REAL(KIND(1d0)), INTENT(in)::wu_m3 ! external water input (e.g., irrigation) [m3]
890  REAL(KIND(1d0)), INTENT(in)::SurfaceArea !Surface area of the study area [m2]
891  REAL(KIND(1D0)), INTENT(IN)::IrrFracPaved!Fraction of paved which are irrigated
892  REAL(KIND(1D0)), INTENT(IN)::IrrFracBldgs!Fraction of buildings (e.g., green roofs) which are irrigated
893  REAL(KIND(1D0)), INTENT(IN)::IrrFracEveTr!Fraction of evergreen trees which are irrigated
894  REAL(KIND(1D0)), INTENT(IN)::IrrFracDecTr!Fraction of deciduous trees which are irrigated
895  REAL(KIND(1D0)), INTENT(IN)::IrrFracGrass!Fraction of grass which is irrigated
896  REAL(KIND(1D0)), INTENT(IN)::IrrFracBSoil!Fraction of bare soil trees which are irrigated
897  REAL(KIND(1D0)), INTENT(IN)::IrrFracWater!Fraction of water which are irrigated
898  REAL(KIND(1d0)), INTENT(in)::InternalWaterUse_h !Internal water use [mm h-1]
899  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::WUProfA_24hr !Automatic water use profiles at hourly scales
900  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::WUProfM_24hr !Manual water use profiles at hourly scales
901  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr!Surface fractions [-]
902 
903  REAL(KIND(1d0)), DIMENSION(12), INTENT(in)::HDD_id !HDD(id-1), Heating Degree Days (see SUEWS_DailyState.f95)
904  REAL(KIND(1d0)), DIMENSION(9), INTENT(in)::WUDay_id!WUDay(id-1), Daily water use for EveTr, DecTr, Grass [mm] (see SUEWS_DailyState.f95)
905 
906  INTEGER, INTENT(in)::DayofWeek_id(3)!DayofWeek(id) 1 - day of week; 2 - month; 3 - season
907  INTEGER, INTENT(in)::WaterUseMethod !Use modelled (0) or observed (1) water use
908  INTEGER, INTENT(in)::NSH!Number of timesteps per hour
909  INTEGER, INTENT(in)::it !Hour
910  INTEGER, INTENT(in)::imin !Minutes
911  INTEGER, INTENT(in)::DLS !day lightsavings =1 + 1h) =0
912 
913  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::wu_nsurf !external Water use for each surface [mm]
914  REAL(KIND(1d0)), INTENT(out)::wu_int !Internal water use for the model timestep [mm] (over whole study area)
915  REAL(KIND(1d0)), INTENT(out)::wu_ext !External water use for the model timestep [mm] (over whole study area)
916 
917  REAL(KIND(1d0))::wu_EveTr !Water use for evergreen trees/shrubs [mm]
918  REAL(KIND(1d0))::wu_DecTr !Water use for deciduous trees/shrubs [mm]
919  REAL(KIND(1d0))::wu_Grass !Water use for grass [mm]
920 
921  REAL(KIND(1d0)), DIMENSION(nsurf)::WUDay_A_id!modelled Automatic Daily water use for each surface [mm] (see SUEWS_DailyState.f95)
922  REAL(KIND(1d0)), DIMENSION(nsurf)::WUDay_M_id!modelled Manual Daily water use for each surface [mm] (see SUEWS_DailyState.f95)
923  REAL(KIND(1d0)), DIMENSION(nsurf)::IrrFrac !faction of irrigated part in each surface [-]
924  REAL(KIND(1d0)), DIMENSION(nsurf)::WUArea !water use area [m2] for each surface type
925 
926  ! REAL(KIND(1d0)):: WUAreaEveTr_m2
927  ! REAL(KIND(1d0)):: WUAreaDecTr_m2
928  ! REAL(KIND(1d0)):: WUAreaGrass_m2
929  REAL(KIND(1d0)):: WUAreaTotal_m2
930  REAL(KIND(1d0)):: InternalWaterUse !Internal water use for the model timestep [mm]
931  REAL(KIND(1d0)):: flag_WuM = 1
932  REAL(KIND(1d0)):: wu !Water use for the model timestep [mm]
933  INTEGER:: ih !Hour corrected for Daylight savings
934  INTEGER:: iu !1=weekday OR 2=weekend
935  INTEGER :: tstep ! timestep in second
936  REAL(KIND(1d0)), PARAMETER::NAN = -999.
937  REAL(KIND(1d0)):: OverUse
938  REAL(KIND(1d0)):: rain_cum_daily ! accumulated daily rainfall
939 
940  REAL(KIND(1d0)):: get_Prof_SpecTime_sum
941 
942  REAL(KIND(1d0)):: WUProfA_tstep ! automatic water use profile value at tstep
943  REAL(KIND(1d0)):: WUProfM_tstep ! mannual water use profile value at tstep
944 
945  ! NB: set OverUse as 0 as done module_constants, TS 22 Oct 2017
946  ! and the logic for calculating OverUse to be determined
947  overuse = 0
948 
949  ! initialise wu
950  wu = 0
951 
952  ! timestep in second
953  tstep = int(3600/nsh)
954 
955  ! accumulated daily rainfall
956  rain_cum_daily = hdd_id(11)
957 
958  ! Irrigated Fraction of each surface
959  ! TS: as of 20191130, assuming irrigation fraction as ONE except for vegetated surfaces
960  ! IrrFrac = 1
961  ! IrrFrac(ConifSurf) = IrrFracEveTr
962  ! IrrFrac(DecidSurf) = IrrFracDecTr
963  ! IrrFrac(GrassSurf) = IrrFracGrass
964 
965  ! Irrigated Fraction of each surface
966  ! TS: 20200409, add irrigation fractions for all surfaces
967  irrfrac = [irrfracpaved, irrfracbldgs, &
968  irrfracevetr, irrfracdectr, irrfracgrass, &
969  irrfracbsoil, irrfracwater]
970 
971  ! --------------------------------------------------------------------------------
972  ! If water used is observed and provided in the met forcing file, units are m3
973  ! Divide observed water use (in m3) by water use area to find water use (in mm)
974  IF (waterusemethod == 1) THEN !If water use is observed
975  ! Calculate water use area [m2] for each surface type
976  ! WUAreaEveTr_m2 = IrrFracEveTr*sfr(ConifSurf)*SurfaceArea
977  ! WUAreaDecTr_m2 = IrrFracDecTr*sfr(DecidSurf)*SurfaceArea
978  ! WUAreaGrass_m2 = IrrFracGrass*sfr(GrassSurf)*SurfaceArea
979  ! WUAreaTotal_m2 = WUAreaEveTr_m2 + WUAreaDecTr_m2 + WUAreaGrass_m2
980 
981  wuarea = irrfrac*sfr*surfacearea
982  wuareatotal_m2 = sum(wuarea)
983 
984  !Set water use [mm] for each surface type to zero initially
985  wu_evetr = 0
986  wu_dectr = 0
987  wu_grass = 0
988 
989  wu_nsurf = 0
990  IF (wu_m3 == nan .OR. wu_m3 == 0) THEN !If no water use
991  ! wu_m3=0
992  wu = 0
993  ELSE !If water use
994  IF (wuareatotal_m2 > 0) THEN
995  wu = (wu_m3/wuareatotal_m2*1000) !Water use in mm for the whole irrigated area
996  ! IF (WUAreaEveTr_m2 > 0) THEN
997  ! wu_EveTr = wu !Water use for Irr EveTr in mm - these are all the same at the moment
998  ! wu_EveTr = wu_EveTr*IrrFracEveTr !Water use for EveTr in mm
999  ! ENDIF
1000  ! IF (WUAreaDecTr_m2 > 0) THEN
1001  ! wu_DecTr = wu !Water use for Irr DecTr in mm - these are all the same at the moment
1002  ! wu_DecTr = wu_DecTr*IrrFracDecTr !Water use for DecTr in mm
1003  ! ENDIF
1004  ! IF (WUAreaGrass_m2 > 0) THEN
1005  ! wu_Grass = wu !Water use for Irr Grass in mm - these are all the same at the moment
1006  ! wu_Grass = wu_Grass*IrrFracGrass !Water use for Grass in mm
1007  ! ENDIF
1008 
1009  wu_nsurf = wu*irrfrac
1010 
1011  wu = (wu_m3/surfacearea*1000) !Water use for the whole study area in mm
1012  ENDIF
1013  ENDIF
1014 
1015  ! --------------------------------------------------------------------------------
1016  ! If water use is modelled, calculate at timestep of model resolution [mm]
1017  ELSEIF (waterusemethod == 0) THEN !If water use is modelled
1018 
1019  ! Account for Daylight saving
1020  ih = it - dls
1021  IF (ih < 0) ih = 23
1022 
1023  ! Weekday or weekend profile
1024  iu = 1 !Set to 1=weekday
1025  ! IF(DayofWeek(id,1)==1.OR.DayofWeek(id,1)==7) THEN
1026  IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) THEN
1027  iu = 2 !Set to 2=weekend
1028  ENDIF
1029 
1030  !write(*,*) (NSH*(ih+1-1)+imin*NSH/60+1)
1031  wuday_a_id = 0
1032  wuday_a_id(conifsurf) = wuday_id(2)
1033  wuday_a_id(decidsurf) = wuday_id(5)
1034  wuday_a_id(grasssurf) = wuday_id(8)
1035 
1036  wuday_m_id = 0
1037  wuday_m_id(conifsurf) = wuday_id(3)
1038  wuday_m_id(decidsurf) = wuday_id(6)
1039  wuday_m_id(grasssurf) = wuday_id(9)
1040 
1041  ! ---- Automatic irrigation ----
1042  ! wu_EveTr = WUProfA_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(2) !Automatic evergreen trees
1043  ! wu_DecTr = WUProfA_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(5) !Automatic deciduous trees
1044  ! wu_Grass = WUProfA_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(8) !Automatic grass
1045  wuprofa_tstep = get_prof_spectime_sum(ih, imin, 0, wuprofa_24hr(:, iu), tstep)
1046  ! wu_EveTr = get_Prof_SpecTime_sum(ih, imin, 0, WUProfA_24hr(:, iu), tstep)*WUDay_id(2) !Automatic evergreen trees
1047  ! wu_DecTr = get_Prof_SpecTime_sum(ih, imin, 0, WUProfA_24hr(:, iu), tstep)*WUDay_id(5) !Automatic deciduous trees
1048  ! wu_Grass = get_Prof_SpecTime_sum(ih, imin, 0, WUProfA_24hr(:, iu), tstep)*WUDay_id(8) !Automatic grass
1049 
1050  ! PRINT*, ''
1051  ! PRINT*, 'WUDay_id(2) ',WUDay_id(2)
1052  ! PRINT*, 'profile ',get_Prof_SpecTime_sum(ih,imin,0,WUProfA_24hr(:,iu),tstep)
1053  ! PRINT*, 'manual:'
1054  ! PRINT*, 'wu_EveTr',wu_EveTr
1055  ! PRINT*, 'wu_DecTr',wu_DecTr
1056  ! PRINT*, 'wu_Grass',wu_Grass
1057 
1058  ! ---- Manual irrigation ----
1059  flag_wum = 1 !Initialize flag_WuM to 1, but if raining, reduce manual fraction of water use
1060  ! If cumulative daily precipitation exceeds 2 mm
1061  IF (rain_cum_daily > 2) THEN !.and.WUDay(id-1,3)>0) then !Commented out HCW 23/01/2015
1062  flag_wum = 0 ! 0 -> No manual irrigation if raining
1063  ENDIF
1064 
1065  ! Add manual to automatic to find total irrigation
1066  ! wu_EveTr = wu_EveTr + (WuFr*WUProfM_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(3)) !Manual evergreen trees
1067  ! wu_DecTr = wu_DecTr + (WuFr*WUProfM_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(6)) !Manual deciduous trees
1068  ! wu_Grass = wu_Grass + (WuFr*WUProfM_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)*WUDay_id(9)) !Manual grass
1069  wuprofm_tstep = get_prof_spectime_sum(ih, imin, 0, wuprofm_24hr(:, iu), tstep)
1070 
1071  ! sum up irrigation amount of automatic and manual approaches
1072  wu_nsurf = wuprofa_tstep*wuday_a_id + wuprofm_tstep*wuday_m_id*flag_wum
1073  ! apply irrigation fraction: part of land covers are not irrigated
1074  wu_nsurf = wu_nsurf*irrfrac
1075 
1076  ! PRINT*, 'auto:'
1077  ! PRINT*, 'wu_EveTr',wu_EveTr
1078  ! PRINT*, 'wu_DecTr',wu_DecTr
1079  ! PRINT*, 'wu_Grass',wu_Grass
1080  ! Added HCW 12 Feb 2015.
1081  !wu_EveTr=wu_EveTr*sfr(ConifSurf)*IrrFracEveTr !Water use for EveTr [mm]
1082  !wu_DecTr=wu_DecTr*sfr(DecidSurf)*IrrFracDecTr !Water use for DecTr [mm]
1083  !wu_Grass=wu_Grass*sfr(GrassSurf)*IrrFracGrass !Water use for Grass [mm]
1084  ! wu_EveTr = wu_EveTr*IrrFracEveTr !Water use for EveTr [mm]
1085  ! wu_DecTr = wu_DecTr*IrrFracDecTr !Water use for DecTr [mm]
1086  ! wu_Grass = wu_Grass*IrrFracGrass !Water use for Grass [mm]
1087 
1088  ! PRINT*, 'auto:'
1089  ! PRINT*, 'IrrFracEveTr',IrrFracEveTr
1090  ! PRINT*, 'IrrFracDecTr',IrrFracDecTr
1091  ! PRINT*, 'IrrFracGrass',IrrFracGrass
1092 
1093  ! Total water use for the whole study area [mm]
1094  ! wu = wu_EveTr*sfr(ConifSurf) + wu_DecTr*sfr(DecidSurf) + wu_Grass*sfr(GrassSurf)
1095  wu = dot_product(wu_nsurf, sfr)
1096 
1097  ENDIF !End WU_choice
1098  ! --------------------------------------------------------------------------------
1099 
1100  ! Internal water use is supplied in SUEWS_Irrigation in mm h-1
1101  ! Convert to mm for the model timestep
1102  internalwateruse = internalwateruse_h/nsh_real
1103 
1104  ! Remove InternalWaterUse from the total water use
1105  wu_ext = wu - (internalwateruse + overuse)
1106  ! Check ext_wu cannot be negative
1107  IF (wu_ext < 0) THEN
1108  overuse = abs(wu_ext)
1109  wu_ext = 0
1110  ELSE
1111  overuse = 0
1112  ENDIF
1113 
1114  wu_int = wu - wu_ext
1115 
1116  ! Decrease the water use for each surface by the same proportion
1117  IF (wu_ext /= 0 .AND. wu /= 0) THEN
1118  ! wu_EveTr = wu_EveTr*wu_ext/wu
1119  ! wu_DecTr = wu_DecTr*wu_ext/wu
1120  ! wu_Grass = wu_Grass*wu_ext/wu
1121  wu_nsurf = wu_nsurf*wu_ext/wu
1122  ENDIF
1123 
1124  END SUBROUTINE suews_cal_wateruse
1125  !===================================================================================
1126 
1127 END MODULE waterdist_module
subroutine suews_update_soilmoist(NonWaterFraction, SoilStoreCap, sfr, soilstore_id, SoilMoistCap, SoilState, vsmd, smd)
subroutine suews_cal_soilstate(SMDMethod, xsmd, NonWaterFraction, SoilMoistCap, SoilStoreCap, surf_chang_per_tstep, soilstore_id, soilstoreOld, sfr, smd, smd_nsurf, tot_chang_per_tstep, SoilState)
subroutine suews_cal_wateruse(nsh_real, wu_m3, SurfaceArea, sfr, 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_nsurf, wu_int, wu_ext)
subroutine redistributewater(snowUse, WaterDist, sfr, Drain, AddWaterRunoff, AddWater)
integer, parameter bsoilsurf
integer, parameter nsurf
integer, parameter conifsurf
subroutine updateflood(is, runoff, sfr, PipeCapacity, RunoffToWater, runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes)
integer, parameter grasssurf
subroutine suews_cal_horizontalsoilwater(sfr, SoilStoreCap, SoilDepth, SatHydraulicConduct, SurfaceArea, NonWaterFraction, tstep_real, soilstore_id, runoffSoil, runoffSoil_per_tstep)
subroutine cal_water_storage(is, sfr, PipeCapacity, RunoffToWater, pin, WU_nsurf, drain, AddWater, addImpervious, nsh_real, stateOld, AddWaterRunoff, PervFraction, addVeg, SoilStoreCap, addWaterBody, FlowChange, StateLimit, runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes, ev, soilstore_id, SurplusEvap, runoffWaterBody, p_mm, chang, runoff, state_id)
integer, parameter excesssurf
integer, parameter decidsurf
integer, parameter pavsurf
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
integer, parameter bldgsurf
integer, parameter watersurf
subroutine drainage(is, state_is, StorCap, DrainEq, DrainCoef1, DrainCoef2, nsh_real, drain_is)