SUEWS API Site
Documentation of SUEWS source code
suews_phys_dailystate.f95
Go to the documentation of this file.
2 USE allocatearray, ONLY: &
4
5 IMPLICIT NONE
6 ! INTEGER,PARAMETER::ndays=366
7 ! INTEGER,PARAMETER::nvegsurf=3
8 ! INTEGER,PARAMETER::ncolumnsDataOutDailyState=46
9
10CONTAINS
11
12 ! Calculation of daily state variables
13 ! Responds to what has happened in the past (temperature, rainfall, etc)
14 ! Updates each time step, but for many variables, correct values are calculated only at the end of each day!
15 ! --> for these variables, the rest of the code MUST use values from the previous day
16 ! N.B. Some of this code is repeated in SUEWS_Initial
17 ! --> so if changes are made here, SUEWS_Initial may also need to be updated accordingly
18 ! N.B. Currently, daily variables are calculated using 00:00-23:55 timestamps (for 5-min resolution); should use 00:05-00:00
19 !
20 ! Last modified:
21 ! TS 09 Jul 2018 - Modified HDD array to hold values for actual calculation
22 ! TS 18 Sep 2017 - Added explicit interface
23 ! TS 07 Jun 2017 - Improve the format of output with more friendly alignment
24 ! HCW 04 Jul 2016 - GridID can now be up to 10 digits long
25 ! HCW 25 May 2016 - Added extra columns to daily state file (albedo for EveTr and Grass)
26 ! HCW 24 May 2016 - Bug fixed in naming of SUEWS_cal_DailyState file (now uses GridIDmatrix(Gridiv) rather than Gridiv)
27 ! LJ 27 Jan 2016 - Removal of tabs
28 ! HCW 20 Aug 2015 - Sign of the porosity change corrected so that porosity is greatest when LAI is smallest
29 ! HCW 03 Jul 2015 - Increased output resolution of P/day in SUEWS_cal_DailyState file to avoid rounding errors.
30 ! Albedo of EveTr and Grass now adjusted based on change in LAI for EveTr and Grass
31 ! (rather than DecTr)
32 ! HCW 29 Jun 2015 - Added albChange for EveTr and Grass surfaces
33 ! HCW 11 Jun 2015 - Bug fix from 05 Jun now fixed in a different way -
34 ! DecidCap is now treated the same as DecidAlb so should be able to cope with multiple grids.
35 ! HCW 05 Jun 2015 - Bug fix - set all current storage capacities (StoreDrainPrm(6,)) to min. value, then set for DecTr
36 ! LJ 11 Mar 2015 - Removed switch as no longer necessary
37 ! HCW 06 Mar 2015 - iy used instead of year which does not have a value here
38 ! HCW 20 Feb 2015 - Added StoreDrainPrm(6,is) for the current storage capacity
39 ! Updated and corrected SUEWS_cal_DailyState output file
40 ! LJ 05 Feb 2015 - SUEWS_cal_DailyState saving fixed. Now header is printed and the file closed and opened as suggested.
41 ! N.B. Bug in daily Precip - needs fixing!!! - HCW thinks this is fixed 20 Feb 2015
42 ! HCW 26 Jan 2015 - sfr_surf and IrrFracs deleted from WUDay calculations, so that WUDay is not spread over
43 ! the total area
44 ! HCW 23 Jan 2015 - WUDay now has 9 columns (EveTr, DecTr, Grass; automatic, manual, total)
45 ! HCW 27 Nov 2014 - Handles values for different grids (Gridiv & ir arguments)
46 ! Added the calculation of surface temperature
47 ! LJ 22 Feb 2013 - Snow albedo aging and calculation of snow density added,
48 ! LJ 22 Jul 2013 - Calculation of LAI senescence from previous day length added
49 ! sg feb 2012 - rewritten from LUMPS_LAI so done in real time
50 !
51 ! To Do
52 ! - Account for change of year in 5-day running mean
53 ! - Check LAI calcs (N/S hemisphere similarities; use of day length)
54 ! - Take out doy limits (140,170, etc) and code as parameters
55 ! - Could add different coefficients (Ie_m, Ie_a) for each vegetation type
56 !==============================================================================
58 iy, id, it, imin, isec, tstep, tstep_prev, dt_since_start, DayofWeek_id, & !input
59 Tmin_id_prev, Tmax_id_prev, lenDay_id_prev, &
60 BaseTMethod, &
61 WaterUseMethod, Ie_start, Ie_end, &
62 LAICalcYes, LAIType, &
63 nsh_real, avkdn, Temp_C, Precip, BaseT_HC, &
64 BaseT_Heating, BaseT_Cooling, &
65 lat, Faut, LAI_obs, &
66 AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, &
67 AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, &
68 CapMax_dec, CapMin_dec, PorMax_dec, PorMin_dec, &
69 Ie_a, Ie_m, DayWatPer, DayWat, &
70 BaseT, BaseTe, GDDFull, SDDFull, LAIMin, LAIMax, LAIPower, &
71 DecidCap_id_prev, StoreDrainPrm_prev, LAI_id_prev, GDD_id_prev, SDD_id_prev, &
72 albDecTr_id_prev, albEveTr_id_prev, albGrass_id_prev, porosity_id_prev, & !input
73 HDD_id_prev, & !input
74 state_id, soilstore_id, SoilStoreCap, H_maintain, & !input
75 HDD_id_next, & !output
76 Tmin_id_next, Tmax_id_next, lenDay_id_next, &
77 albDecTr_id_next, albEveTr_id_next, albGrass_id_next, porosity_id_next, & !output
78 DecidCap_id_next, StoreDrainPrm_next, LAI_id_next, GDD_id_next, SDD_id_next, deltaLAI, WUDay_id) !output
79
80 ! USE Snow_module, ONLY: SnowUpdate
81 USE datetime_module, ONLY: datetime, timedelta
82
83 IMPLICIT NONE
84
85 INTEGER, INTENT(IN) :: iy
86 INTEGER, INTENT(IN) :: id
87 INTEGER, INTENT(IN) :: it
88 INTEGER, INTENT(IN) :: imin
89 INTEGER, INTENT(IN) :: isec
90 INTEGER, INTENT(IN) :: tstep
91 INTEGER, INTENT(IN) :: tstep_prev
92 INTEGER, INTENT(IN) :: dt_since_start
93
94 INTEGER, INTENT(IN) :: WaterUseMethod
95 INTEGER, INTENT(IN) :: BaseTMethod
96 INTEGER, INTENT(IN) :: Ie_start !Starting time of water use (DOY)
97 INTEGER, INTENT(IN) :: Ie_end !Ending time of water use (DOY)
98 INTEGER, INTENT(IN) :: LAICalcYes
99
100 INTEGER, DIMENSION(nvegsurf), INTENT(IN) :: LAIType !LAI equation to use: original (0) or new (1)
101
102 REAL(KIND(1D0)), INTENT(IN) :: nsh_real
103 REAL(KIND(1D0)), INTENT(IN) :: avkdn
104 REAL(KIND(1D0)), INTENT(IN) :: Temp_C
105 REAL(KIND(1D0)), INTENT(IN) :: Precip
106 REAL(KIND(1D0)), INTENT(IN) :: BaseT_HC
107 REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) :: BaseT_Heating
108 REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) :: BaseT_Cooling
109 REAL(KIND(1D0)), INTENT(IN) :: lat
110 REAL(KIND(1D0)), INTENT(IN) :: Faut
111 REAL(KIND(1D0)), INTENT(IN) :: LAI_obs
112 ! REAL(KIND(1D0)), INTENT(IN)::tau_a
113 ! REAL(KIND(1D0)), INTENT(IN)::tau_f
114 ! REAL(KIND(1D0)), INTENT(IN)::tau_r
115 ! REAL(KIND(1D0)), INTENT(IN)::SnowDensMax
116 ! REAL(KIND(1D0)), INTENT(IN)::SnowDensMin
117 ! REAL(KIND(1D0)), INTENT(in)::SnowAlbMax
118 ! REAL(KIND(1D0)), INTENT(IN)::SnowAlbMin
119 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_DecTr
120 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_EveTr
121 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_Grass
122 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_DecTr
123 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_EveTr
124 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_Grass
125 REAL(KIND(1D0)), INTENT(IN) :: CapMax_dec
126 REAL(KIND(1D0)), INTENT(IN) :: CapMin_dec
127 REAL(KIND(1D0)), INTENT(IN) :: PorMax_dec
128 REAL(KIND(1D0)), INTENT(IN) :: PorMin_dec
129 ! REAL(KIND(1d0)),INTENT(IN) ::VegPhenLumps
130
131 REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) :: Ie_a
132 REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) :: Ie_m !Coefficients for automatic and manual irrigation models
133 REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) :: DayWatPer !% of houses following daily water
134 REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) :: DayWat !Days of watering allowed
135
136 ! ponding-water related
137 REAL(KIND(1D0)), INTENT(IN) :: H_maintain ! ponding water depth to maintain [mm]
138 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: state_id ! surface wetness [mm]
139 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: soilstore_id ! soil water store [mm]
140 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap !Capacity of soil store for each surface [mm]
141
142 ! REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(IN) ::SnowPack
143 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: BaseT !Base temperature for growing degree days [degC]
144 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: BaseTe !Base temperature for senescence degree days [degC]
145 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: GDDFull !Growing degree days needed for full capacity [degC]
146 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: SDDFull !Senescence degree days needed to initiate leaf off [degC]
147 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMin !Min LAI [m2 m-2]
148 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMax !Max LAI [m2 m-2]
149 REAL(KIND(1D0)), DIMENSION(4, nvegsurf), INTENT(IN) :: LAIPower !Coeffs for LAI equation: 1,2 - leaf growth; 3,4 - leaf off
150
151 ! REAL(KIND(1d0)), INTENT(INOUT)::SnowAlb
152
153 ! Growing Degree Days
154 REAL(KIND(1D0)), DIMENSION(3) :: GDD_id ! Growing Degree Days (see SUEWS_DailyState.f95)
155 REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) :: GDD_id_prev ! Growing Degree Days (see SUEWS_DailyState.f95)
156 REAL(KIND(1D0)), DIMENSION(3), INTENT(OUT) :: GDD_id_next ! Growing Degree Days (see SUEWS_DailyState.f95)
157
158 ! Senescence Degree Days
159 REAL(KIND(1D0)), DIMENSION(3) :: SDD_id ! Senescence Degree Days (see SUEWS_DailyState.f95)
160 REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) :: SDD_id_prev ! Senescence Degree Days (see SUEWS_DailyState.f95)
161 REAL(KIND(1D0)), DIMENSION(3), INTENT(OUT) :: SDD_id_next ! Senescence Degree Days (see SUEWS_DailyState.f95)
162
163 ! Daily min temp [degC]
164 REAL(KIND(1D0)) :: Tmin_id
165 REAL(KIND(1D0)), INTENT(IN) :: Tmin_id_prev
166 REAL(KIND(1D0)), INTENT(out) :: Tmin_id_next
167
168 ! Daily max temp [degC]
169 REAL(KIND(1D0)) :: Tmax_id
170 REAL(KIND(1D0)), INTENT(in) :: Tmax_id_prev
171 REAL(KIND(1D0)), INTENT(out) :: Tmax_id_next
172
173 ! Daytime hours [h]
174 REAL(KIND(1D0)) :: lenDay_id
175 REAL(KIND(1D0)), INTENT(IN) :: lenDay_id_prev
176 REAL(KIND(1D0)), INTENT(out) :: lenDay_id_next
177
178 ! LAI for each veg surface [m2 m-2]
179 REAL(KIND(1D0)), DIMENSION(3) :: LAI_id ! LAI for each veg surface [m2 m-2]
180 REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) :: LAI_id_prev ! LAI for each veg surface [m2 m-2]
181 REAL(KIND(1D0)), DIMENSION(3), INTENT(OUT) :: LAI_id_next ! LAI for each veg surface [m2 m-2]
182
183 ! ------------- Key to daily arrays ----------------------------------------------
184 ! TS, 27 Dec 2018: updated the annotation for 2018b and WRF-SUEWS coupling
185
186 ! Heating Degree Days
187 REAL(KIND(1D0)), DIMENSION(12) :: HDD_id ! Heating Degree Days (see SUEWS_DailyState.f95)
188 REAL(KIND(1D0)), DIMENSION(12), INTENT(IN) :: HDD_id_prev ! Heating Degree Days (see SUEWS_DailyState.f95)
189 REAL(KIND(1D0)), DIMENSION(12), INTENT(OUT) :: HDD_id_next ! Heating Degree Days (see SUEWS_DailyState.f95)
190 ! HDD_id:
191 ! first half used for update through the day
192 ! HDD_id(1) ---- Heating [degC]: used for accumulation during calculation
193 ! HDD_id(2) ---- Cooling [degC]: used for accumulation during calculation
194 ! HDD_id(3) ---- Daily mean temp [degC]: used for accumulation during calculation
195 ! HDD_id(4) ---- 5-day running mean temp [degC]: used for actual calculation
196 ! HDD_id(5) ---- Daily precip total [mm]
197 ! HDD_id(6) ---- Days since rain [d]
198 ! second half used for storage of the first half for the prevous day
199 ! HDD_id(6+1) ---- Heating [degC]: used for accumulation during calculation
200 ! HDD_id(6+2) ---- Cooling [degC]: used for accumulation during calculation
201 ! HDD_id(6+3) ---- Daily mean temp [degC]: used for accumulation during calculation
202 ! HDD_id(6+4) ---- 5-day running mean temp [degC]: used for actual calculation
203 ! HDD_id(6+5) ---- Daily precip total [mm]
204 ! HDD_id(6+6) ---- Days since rain [d]
205 ! --------------------------------------------------------------------------------
206
207 ! --------------------------------------------------------------------------------
208 !Daily water use for EveTr, DecTr, Grass [mm] (see SUEWS_DailyState.f95)
209 REAL(KIND(1D0)), DIMENSION(9), INTENT(OUT) :: WUDay_id ! Water use related array
210 ! WUDay_id:
211 ! WUDay_id(1) - Daily water use total for Irr EveTr (automatic+manual) [mm]
212 ! WUDay_id(2) - Automatic irrigation for Irr EveTr [mm]
213 ! WUDay_id(3) - Manual irrigation for Irr EveTr [mm]
214 ! WUDay_id(4) - Daily water use total for Irr DecTr (automatic+manual) [mm]
215 ! WUDay_id(5) - Automatic irrigation for Irr DecTr [mm]
216 ! WUDay_id(6) - Manual irrigation for Irr DecTr [mm]
217 ! WUDay_id(7) - Daily water use total for Irr Grass (automatic+manual) [mm]
218 ! WUDay_id(8) - Automatic irrigation for Irr Grass [mm]
219 ! WUDay_id(9) - Manual irrigation for Irr Grass [mm]
220 ! --------------------------------------------------------------------------------
221
222 ! REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(INOUT)::SnowDens
223 INTEGER, DIMENSION(3), INTENT(in) :: DayofWeek_id
224
225 REAL(KIND(1D0)), INTENT(OUT) :: deltaLAI
226
227 REAL(KIND(1D0)) :: DecidCap_id
228 REAL(KIND(1D0)), INTENT(IN) :: DecidCap_id_prev
229 REAL(KIND(1D0)), INTENT(OUT) :: DecidCap_id_next
230 REAL(KIND(1D0)) :: albDecTr_id
231 REAL(KIND(1D0)), INTENT(IN) :: albDecTr_id_prev
232 REAL(KIND(1D0)), INTENT(OUT) :: albDecTr_id_next
233 REAL(KIND(1D0)) :: albEveTr_id
234 REAL(KIND(1D0)), INTENT(IN) :: albEveTr_id_prev
235 REAL(KIND(1D0)), INTENT(OUT) :: albEveTr_id_next
236 REAL(KIND(1D0)) :: albGrass_id
237 REAL(KIND(1D0)), INTENT(IN) :: albGrass_id_prev
238 REAL(KIND(1D0)), INTENT(OUT) :: albGrass_id_next
239 REAL(KIND(1D0)) :: porosity_id
240 REAL(KIND(1D0)), INTENT(INOUT) :: porosity_id_prev
241 REAL(KIND(1D0)), INTENT(INOUT) :: porosity_id_next
242 REAL(KIND(1D0)), DIMENSION(6, nsurf) :: StoreDrainPrm
243 REAL(KIND(1D0)), DIMENSION(6, nsurf), INTENT(in) :: StoreDrainPrm_prev
244 REAL(KIND(1D0)), DIMENSION(6, nsurf), INTENT(out) :: StoreDrainPrm_next
245
246 LOGICAL :: first_tstep_Q ! if this is the first tstep of a day
247 LOGICAL :: last_tstep_Q ! if this is the last tstep of a day
248 TYPE(datetime) :: time_now, time_prev, time_next
249
250 ! transfer values
251 lai_id = lai_id_prev
252 gdd_id = gdd_id_prev
253 sdd_id = sdd_id_prev
254 tmin_id = tmin_id_prev
255 tmax_id = tmax_id_prev
256 lenday_id = lenday_id_prev
257 storedrainprm = storedrainprm_prev
258 decidcap_id = decidcap_id_prev
259 albdectr_id = albdectr_id_prev
260 albevetr_id = albevetr_id_prev
261 albgrass_id = albgrass_id_prev
262 porosity_id = porosity_id_prev
263 hdd_id = hdd_id_prev
264
265 ! get timestamps
266 time_now = datetime(year=iy) + timedelta(days=id - 1, hours=it, minutes=imin, seconds=isec)
267 time_prev = time_now - timedelta(seconds=tstep_prev)
268 time_next = time_now + timedelta(seconds=tstep)
269
270 ! test if time at now is the first/last tstep of today
271 first_tstep_q = time_now%getDay() /= time_prev%getDay()
272 last_tstep_q = time_now%getDay() /= time_next%getDay()
273
274 ! --------------------------------------------------------------------------------
275 ! On first timestep of each day, define whether the day each a workday or weekend
276 IF (first_tstep_q) THEN
278 it, imin, & !input
279 hdd_id) !inout
280
281 ! reset certain GDD columns
282 tmin_id = temp_c !Daily min T in column 3
283 tmax_id = temp_c !Daily max T in column 4
284 lenday_id = 0 !Cumulate daytime hours
285 END IF
286
287 ! --------------------------------------------------------------------------------
288 ! regular update at all timesteps of a day
290 basetmethod, &
291 dayofweek_id, &
292 avkdn, & !input
293 temp_c, &
294 precip, &
295 baset_hc, &
296 baset_heating, baset_cooling, &
297 nsh_real, &
298 tmin_id, tmax_id, lenday_id, & !inout
299 hdd_id) !inout
300
301 ! Update snow density, albedo surface fraction
302 ! IF (SnowUse == 1) CALL SnowUpdate( &
303 ! nsurf, tstep, Temp_C, tau_a, tau_f, tau_r, &!input
304 ! SnowDensMax, SnowDensMin, SnowAlbMax, SnowAlbMin, SnowPack, &
305 ! SnowAlb, SnowDens)!inout
306
307 ! --------------------------------------------------------------------------------
308 ! On last timestep, perform the daily calculations -------------------------------
309 ! Daily values not correct until end of each day,
310 ! so main program should use values from the previous day
311 IF (last_tstep_q) THEN
313 id, it, imin, tstep, dt_since_start, & !input
314 tmin_id, tmax_id, lenday_id, &
315 laitype, ie_end, ie_start, laicalcyes, &
316 waterusemethod, dayofweek_id, &
317 albmax_dectr, albmax_evetr, albmax_grass, albmin_dectr, albmin_evetr, albmin_grass, &
318 baset, basete, capmax_dec, capmin_dec, daywat, daywatper, faut, gddfull, &
319 ie_a, ie_m, laimax, laimin, laipower, lat, pormax_dec, pormin_dec, sddfull, lai_obs, &
320 state_id, soilstore_id, soilstorecap, h_maintain, & !input
321 gdd_id, sdd_id, & !inout
322 hdd_id, &
323 lai_id, &
324 decidcap_id, &
325 albdectr_id, &
326 albevetr_id, &
327 albgrass_id, &
328 porosity_id, &
329 storedrainprm, &
330 wuday_id, deltalai) !output
331 END IF !End of section done only at the end of each day (i.e. only once per day)
332
333 ! translate values back
334 lai_id_next = lai_id
335 gdd_id_next = gdd_id
336 sdd_id_next = sdd_id
337 tmin_id_next = tmin_id
338 tmax_id_next = tmax_id
339 lenday_id_next = lenday_id
340 storedrainprm_next = storedrainprm
341 decidcap_id_next = decidcap_id
342 albdectr_id_next = albdectr_id
343 albevetr_id_next = albevetr_id
344 albgrass_id_next = albgrass_id
345 porosity_id_next = porosity_id
346 hdd_id_next = hdd_id
347 ! PRINT*, 'after_DailyState', iy,id,it,imin
348 ! PRINT*, 'HDD(id)', HDD(id,:)
349 ! PRINT*, 'HDD_id', HDD_id
350
351 ! RETURN
352
353 END SUBROUTINE suews_cal_dailystate
354
356 id, it, imin, tstep, dt_since_start, & !input
357 Tmin_id, Tmax_id, lenDay_id, &
358 LAIType, Ie_end, Ie_start, LAICalcYes, &
359 WaterUseMethod, DayofWeek_id, &
360 AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, &
361 BaseT, BaseTe, CapMax_dec, CapMin_dec, DayWat, DayWatPer, Faut, GDDFull, &
362 Ie_a, Ie_m, LAIMax, LAIMin, LAIPower, lat, PorMax_dec, PorMin_dec, SDDFull, LAI_obs, &
363 state_id, soilstore_id, SoilStoreCap, H_maintain, & !input
364 GDD_id, SDD_id, & !inout
365 HDD_id, &
366 LAI_id, &
367 DecidCap_id, &
368 albDecTr_id, &
369 albEveTr_id, &
370 albGrass_id, &
371 porosity_id, &
372 StoreDrainPrm, &
373 WUDay_id, deltaLAI) !output
374 IMPLICIT NONE
375
376 INTEGER, INTENT(IN) :: id
377 INTEGER, INTENT(IN) :: it
378 INTEGER, INTENT(IN) :: imin
379 INTEGER, INTENT(IN) :: tstep
380 INTEGER, INTENT(IN) :: dt_since_start
381 INTEGER, INTENT(IN) :: LAIType(nvegsurf)
382 INTEGER, INTENT(IN) :: Ie_end
383 INTEGER, INTENT(IN) :: Ie_start
384 INTEGER, INTENT(IN) :: LAICalcYes
385 INTEGER, INTENT(IN) :: WaterUseMethod
386 INTEGER, INTENT(in) :: DayofWeek_id(3)
387
388 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_DecTr
389 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_EveTr
390 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_Grass
391 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_DecTr
392 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_EveTr
393 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_Grass
394 REAL(KIND(1D0)), INTENT(IN) :: BaseT(nvegsurf)
395 REAL(KIND(1D0)), INTENT(IN) :: BaseTe(nvegsurf)
396 REAL(KIND(1D0)), INTENT(IN) :: CapMax_dec
397 REAL(KIND(1D0)), INTENT(IN) :: CapMin_dec
398 REAL(KIND(1D0)), INTENT(IN) :: DayWat(7)
399 REAL(KIND(1D0)), INTENT(IN) :: DayWatPer(7)
400 REAL(KIND(1D0)), INTENT(IN) :: Faut
401 REAL(KIND(1D0)), INTENT(IN) :: GDDFull(nvegsurf)
402 REAL(KIND(1D0)), INTENT(IN) :: Ie_a(3)
403 REAL(KIND(1D0)), INTENT(IN) :: Ie_m(3)
404 REAL(KIND(1D0)), INTENT(IN) :: LAIMax(nvegsurf)
405 REAL(KIND(1D0)), INTENT(IN) :: LAIMin(nvegsurf)
406 REAL(KIND(1D0)), INTENT(IN) :: LAIPower(4, nvegsurf)
407 REAL(KIND(1D0)), INTENT(IN) :: lat
408 REAL(KIND(1D0)), INTENT(IN) :: PorMax_dec
409 REAL(KIND(1D0)), INTENT(IN) :: PorMin_dec
410 REAL(KIND(1D0)), INTENT(IN) :: SDDFull(nvegsurf)
411 REAL(KIND(1D0)), INTENT(IN) :: LAI_obs
412 REAL(KIND(1D0)), INTENT(IN) :: Tmin_id
413 REAL(KIND(1D0)), INTENT(IN) :: Tmax_id
414 REAL(KIND(1D0)), INTENT(IN) :: lenDay_id
415 REAL(KIND(1D0)), INTENT(IN) :: H_maintain ! ponding water depth to maintain
416 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: state_id ! surface wetness [mm]
417 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: soilstore_id ! soil water store [mm]
418 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap !Capacity of soil store for each surface [mm]
419
420 REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: GDD_id ! Growing Degree Days (see SUEWS_DailyState.f95)
421 REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: SDD_id ! Senescence Degree Days (see SUEWS_DailyState.f95)
422 REAL(KIND(1D0)), DIMENSION(12), INTENT(INOUT) :: HDD_id
423 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(INOUT) :: LAI_id ! LAI for each veg surface [m2 m-2]
424
425 ! REAL(KIND(1d0)),DIMENSION(6),INTENT(INOUT)::HDD_id_use ! HDD of previous day
426 REAL(KIND(1D0)), DIMENSION(nvegsurf) :: LAI_id_in ! LAI of previous day
427
428 REAL(KIND(1D0)), DIMENSION(9), INTENT(OUT) :: WUDay_id
429 REAL(KIND(1D0)), INTENT(OUT) :: deltaLAI
430
431 REAL(KIND(1D0)), INTENT(INOUT) :: DecidCap_id
432 REAL(KIND(1D0)), INTENT(INOUT) :: albDecTr_id
433 REAL(KIND(1D0)), INTENT(INOUT) :: albEveTr_id
434 REAL(KIND(1D0)), INTENT(INOUT) :: albGrass_id
435 REAL(KIND(1D0)), INTENT(INOUT) :: porosity_id
436
437 REAL(KIND(1D0)), DIMENSION(6, nsurf), INTENT(inout) :: StoreDrainPrm
438
439 ! Calculate heating degree days ------------------------------------------
440 CALL update_hdd( &
441 dt_since_start, it, imin, tstep, & !input
442 hdd_id) !inout
443
444 ! Calculate modelled daily water use ------------------------------------------
445 CALL update_wateruse( &
446 id, waterusemethod, dayofweek_id, lat, faut, hdd_id, & !input
447 state_id, soilstore_id, soilstorecap, h_maintain, & !input
448 ie_a, ie_m, ie_start, ie_end, daywatper, daywat, &
449 wuday_id) !output
450
451 ! PRINT*, ''
452 ! PRINT*, 'WUDay(id)',WUDay(id,:)
453 ! PRINT*, 'WUDay_id',WUDay_id
454
455 !------------------------------------------------------------------------------
456 ! Calculation of LAI from growing degree days
457 ! This was revised and checked on 16 Feb 2014 by LJ
458 !------------------------------------------------------------------------------
459 ! save initial LAI_id
460 lai_id_in = lai_id
461
462 CALL update_gddlai( &
463 id, laicalcyes, & !input
464 lat, lai_obs, &
465 tmin_id, tmax_id, lenday_id, &
466 baset, basete, &
467 gddfull, sddfull, &
468 laimin, laimax, laipower, laitype, &
469 lai_id_in, &
470 gdd_id, sdd_id, & !inout
471 lai_id) !output
472
473 CALL update_veg( &
474 laimax, laimin, & !input
475 albmax_dectr, albmax_evetr, albmax_grass, &
476 albmin_dectr, albmin_evetr, albmin_grass, &
477 capmax_dec, capmin_dec, &
478 pormax_dec, pormin_dec, &
479 lai_id, lai_id_in, &
480 decidcap_id, & !inout
481 albdectr_id, &
482 albevetr_id, &
483 albgrass_id, &
484 porosity_id, &
485 storedrainprm, &
486 deltalai) !output
487
488 ! PRINT*, 'DecidCap',DecidCap(id),DecidCap_id
489 ! PRINT*, 'albDecTr',albDecTr(id),albDecTr_id
490 ! PRINT*, 'albEveTr',albEveTr(id),albEveTr_id
491 ! PRINT*, 'albGrass',albGrass(id),albGrass_id
492 ! PRINT*, 'porosity',porosity(id),porosity_id
493
494 END SUBROUTINE update_dailystate_end
495
497 BaseTMethod, &
498 DayofWeek_id, &
499 avkdn, & !input
500 Temp_C, &
501 Precip, &
502 BaseT_HC, &
503 BaseT_Heating, BaseT_Cooling, &
504 nsh_real, &
505 Tmin_id, Tmax_id, lenDay_id, & !inout
506 HDD_id) !inout
507 ! use time, only: id, id_prev_t
508 IMPLICIT NONE
509
510 INTEGER, INTENT(IN) :: BaseTMethod
511 INTEGER, DIMENSION(3), INTENT(in) :: DayofWeek_id
512
513 REAL(KIND(1D0)), INTENT(IN) :: avkdn
514 REAL(KIND(1D0)), INTENT(IN) :: Temp_C
515 REAL(KIND(1D0)), INTENT(IN) :: Precip
516 REAL(KIND(1D0)), INTENT(IN) :: BaseT_HC
517 REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) :: BaseT_Heating
518 REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) :: BaseT_Cooling
519 REAL(KIND(1D0)), INTENT(IN) :: nsh_real
520 REAL(KIND(1D0)), INTENT(INOUT) :: Tmin_id
521 REAL(KIND(1D0)), INTENT(INOUT) :: Tmax_id
522 REAL(KIND(1D0)), INTENT(INOUT) :: lenDay_id
523 ! REAL(KIND(1d0)), INTENT(out)::Tmin_id_next
524 ! REAL(KIND(1d0)), INTENT(out)::Tmax_id_next
525 ! REAL(KIND(1d0)), INTENT(out)::lenDay_id_next
526
527 ! REAL(KIND(1d0))::tstepcount
528 ! REAL(KIND(1d0)),DIMENSION(-4:366,6),INTENT(INOUT):: HDD
529 ! REAL(KIND(1d0)), DIMENSION(5), INTENT(INOUT):: GDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
530 REAL(KIND(1D0)), DIMENSION(12), INTENT(INOUT) :: HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
531 ! REAL(KIND(1d0)),DIMENSION(5),INTENT(OUT):: GDD_id_prev !Growing Degree Days (see SUEWS_DailyState.f95)
532 INTEGER :: iu ! flag for weekday/weekend
533
534 REAL(KIND(1D0)) :: dT_heating
535 REAL(KIND(1D0)) :: dT_cooling
536
537 REAL(KIND(1D0)) :: BaseT_Heating_use
538 REAL(KIND(1D0)) :: BaseT_Cooling_use
539
540 ! Set weekday/weekend counter
541 iu = 1 !Set to 1=weekday
542 IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) iu = 2 !Set to 2=weekend
543
544 SELECT CASE (basetmethod)
545 CASE (1)
546 baset_heating_use = baset_hc
547 baset_cooling_use = baset_hc
548 CASE (2)
549 baset_heating_use = baset_heating(iu)
550 baset_cooling_use = baset_cooling(iu)
551
552 CASE default
553 CALL errorhint(75, "RunControl.nml", -999, -999, -999)
554
555 END SELECT
556
557 ! Daily min and max temp (these get updated through the day) ---------------------
558 tmin_id = min(temp_c, tmin_id) !Daily min T in column 3
559 tmax_id = max(temp_c, tmax_id) !Daily max T in column 4
560 IF (avkdn > 10) THEN
561 lenday_id = lenday_id + 1/nsh_real !Cumulate daytime hours !Divide by nsh (HCW 01 Dec 2014)
562 END IF
563
564 ! Calculations related to heating and cooling degree days (HDD) ------------------
565 ! See Sailor & Vasireddy (2006) EMS Eq 1,2 (theirs is hourly timestep)
566 dt_heating = baset_heating_use - temp_c
567 dt_cooling = temp_c - baset_cooling_use
568
569 hdd_id(1) = hdd_id(1) + merge(dt_heating, 0d0, dt_heating >= 0) !Heating
570 hdd_id(2) = hdd_id(2) + merge(dt_cooling, 0d0, dt_cooling >= 0) !Cooling
571 hdd_id(3) = hdd_id(3) + temp_c !Will become daily average temperature
572 ! 4 ------------------------------------! !5-day running mean
573 hdd_id(5) = hdd_id(5) + precip !Daily precip total
574 ! 6 ------------------------------------! !Days since rain
575
576 END SUBROUTINE update_dailystate_day
577
578 SUBROUTINE update_veg( &
579 LAImax, LAIMin, & !input
580 AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, &
581 AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, &
582 CapMax_dec, CapMin_dec, &
583 PorMax_dec, PorMin_dec, &
584 LAI_id, LAI_id_prev, &
585 DecidCap_id, & !inout
586 albDecTr_id, &
587 albEveTr_id, &
588 albGrass_id, &
589 porosity_id, &
590 StoreDrainPrm, &
591 deltaLAI) !output
592
593 IMPLICIT NONE
594
595 ! INTEGER,INTENT(IN)::id
596 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAImax
597 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMin
598
599 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_DecTr
600 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_EveTr
601 REAL(KIND(1D0)), INTENT(IN) :: AlbMax_Grass
602 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_DecTr
603 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_EveTr
604 REAL(KIND(1D0)), INTENT(IN) :: AlbMin_Grass
605 REAL(KIND(1D0)), INTENT(IN) :: CapMax_dec
606 REAL(KIND(1D0)), INTENT(IN) :: CapMin_dec
607 REAL(KIND(1D0)), INTENT(IN) :: PorMax_dec
608 REAL(KIND(1D0)), INTENT(IN) :: PorMin_dec
609 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAI_id, LAI_id_prev
610
611 REAL(KIND(1D0)), INTENT(INOUT) :: DecidCap_id
612 REAL(KIND(1D0)), INTENT(INOUT) :: albDecTr_id
613 REAL(KIND(1D0)), INTENT(INOUT) :: albEveTr_id
614 REAL(KIND(1D0)), INTENT(INOUT) :: albGrass_id
615 REAL(KIND(1D0)), INTENT(INOUT) :: porosity_id
616
617 REAL(KIND(1D0)), DIMENSION(6, nsurf), INTENT(inout) :: StoreDrainPrm
618
619 REAL(KIND(1D0)), INTENT(OUT) :: deltaLAI
620
621 INTEGER :: iv
622
623 REAL(KIND(1D0)) :: albChangeDecTr
624 REAL(KIND(1D0)) :: albChangeEveTr
625 REAL(KIND(1D0)) :: albChangeGrass
626 REAL(KIND(1D0)) :: CapChange
627
628 REAL(KIND(1D0)) :: deltaLAIEveTr
629 REAL(KIND(1D0)) :: deltaLAIGrass
630 REAL(KIND(1D0)) :: porChange
631 !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
632 ! Calculate the development of vegetation cover
633 ! Albedo changes with LAI for each vegetation type
634 ! Storage capacity and porosity are updated based on DecTr LAI only (seasonal variation in Grass and EveTr assumed small)
635 ! If only LUMPS is used, set deciduous capacities to 0
636 ! QUESTION: Assume porosity Change based on GO99- Heisler?
637 deltalai = 0
638 deltalaievetr = 0
639 deltalaigrass = 0
640 capchange = 0
641 porchange = 0
642 albchangedectr = 0
643 albchangeevetr = 0
644 albchangegrass = 0
645
646 iv = ivdecid
647 IF ((lai_id(iv) - lai_id_prev(iv)) /= 0) THEN
648 deltalai = (lai_id(iv) - lai_id_prev(iv))/(laimax(iv) - laimin(iv))
649 albchangedectr = (albmax_dectr - albmin_dectr)*deltalai
650 capchange = (capmin_dec - capmax_dec)*deltalai
651 porchange = (pormin_dec - pormax_dec)*deltalai
652 END IF
653
654 iv = ivconif
655 IF ((lai_id(iv) - lai_id_prev(iv)) /= 0) THEN
656 deltalaievetr = (lai_id(iv) - lai_id_prev(iv))/(laimax(iv) - laimin(iv))
657 albchangeevetr = (albmax_evetr - albmin_evetr)*deltalaievetr
658 END IF
659
660 iv = ivgrass
661 IF ((lai_id(iv) - lai_id_prev(iv)) /= 0) THEN
662 deltalaigrass = (lai_id(iv) - lai_id_prev(iv))/(laimax(iv) - laimin(iv))
663 albchangegrass = (albmax_grass - albmin_grass)*deltalaigrass
664 END IF
665
666 iv = ivdecid
667
668 !write(*,*) deltaLAI, deltaLAIEveTr, deltaLAIGrass
669
670 decidcap_id = decidcap_id - capchange
671 storedrainprm(6, decidsurf) = decidcap_id !Change current storage capacity of deciduous trees
672 porosity_id = porosity_id + porchange !- changed to + by HCW 20 Aug 2015 (porosity greatest when LAI smallest)
673
674 ! update albedo values while limiting these to valid ranges
675 !albDecTr_id = min(max(albDecTr_id + albChangeDecTr, AlbMin_DecTr), AlbMax_DecTr)
676 !albEveTr_id = min(max(albEveTr_id + albChangeEveTr, AlbMin_EveTr), AlbMax_EveTr)
677 !albGrass_id = min(max(albGrass_id + albChangeGrass, AlbMin_Grass), AlbMax_Grass)
678 albdectr_id = albdectr_id + albchangedectr
679 albevetr_id = albevetr_id + albchangeevetr
680 albgrass_id = albgrass_id + albchangegrass
681
682 END SUBROUTINE update_veg
683
684 SUBROUTINE update_gddlai( &
685 id, LAICalcYes, & !input
686 lat, LAI_obs, &
687 Tmin_id_prev, Tmax_id_prev, lenDay_id_prev, &
688 BaseT, BaseTe, &
689 GDDFull, SDDFull, &
690 LAIMin, LAIMax, LAIPower, LAIType, &
691 LAI_id_prev, &
692 GDD_id, SDD_id, & !inout
693 LAI_id_next) !output
694 IMPLICIT NONE
695
696 !------------------------------------------------------------------------------
697 ! Calculation of LAI from growing degree days
698 ! This was revised and checked on 16 Feb 2014 by LJ
699 !------------------------------------------------------------------------------
700
701 INTEGER, INTENT(IN) :: id
702 INTEGER, INTENT(IN) :: LAICalcYes
703
704 REAL(KIND(1D0)), INTENT(IN) :: lat
705 REAL(KIND(1D0)), INTENT(IN) :: LAI_obs
706 REAL(KIND(1D0)), INTENT(IN) :: Tmin_id_prev
707 REAL(KIND(1D0)), INTENT(IN) :: Tmax_id_prev
708 REAL(KIND(1D0)), INTENT(IN) :: lenDay_id_prev
709
710 ! --- Vegetation phenology ---------------------------------------------------------------------
711 ! Parameters provided in input information for each vegetation surface (SUEWS_Veg.txt)
712 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: BaseT !Base temperature for growing degree days [degC]
713 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: BaseTe !Base temperature for senescence degree days [degC]
714 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: GDDFull !Growing degree days needed for full capacity [degC]
715 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: SDDFull !Senescence degree days needed to initiate leaf off [degC]
716 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMin !Min LAI [m2 m-2]
717 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMax !Max LAI [m2 m-2]
718 REAL(KIND(1D0)), DIMENSION(4, nvegsurf), INTENT(IN) :: LAIPower !Coeffs for LAI equation: 1,2 - leaf growth; 3,4 - leaf off
719 !! N.B. currently DecTr only, although input provided for all veg types
720 INTEGER, DIMENSION(nvegsurf), INTENT(IN) :: LAIType !LAI equation to use: original (0) or new (1)
721
722 REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: GDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
723 REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: SDD_id !Senescence Degree Days (see SUEWS_DailyState.f95)
724 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(OUT) :: LAI_id_next !LAI for each veg surface [m2 m-2]
725 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAI_id_prev ! LAI of previous day
726
727 REAL(KIND(1D0)) :: delta_SDD !Switches and checks for GDD
728 REAL(KIND(1D0)) :: delta_GDD !Switches and checks for GDD
729 REAL(KIND(1D0)) :: indHelp !Switches and checks for GDD
730 REAL(KIND(1D0)), DIMENSION(3) :: GDD_id_prev ! GDD of previous day
731 REAL(KIND(1D0)), DIMENSION(3) :: SDD_id_prev ! SDD of previous day
732
733 INTEGER :: critDays
734 INTEGER :: iv
735
736 ! translate values of previous day to local variables
737 gdd_id_prev = gdd_id
738 sdd_id_prev = sdd_id
739 ! LAI_id_prev = LAI_id_next
740
741 critdays = 50 !Critical limit for GDD when GDD or SDD is set to zero
742
743 ! Loop through vegetation types (iv)
744 DO iv = 1, nvegsurf
745 ! Calculate GDD for each day from the minimum and maximum air temperature
746 delta_gdd = ((tmin_id_prev + tmax_id_prev)/2 - baset(iv)) !Leaf on
747 delta_sdd = ((tmin_id_prev + tmax_id_prev)/2 - basete(iv)) !Leaf off
748
749 indhelp = 0 !Help switch to allow GDD to go to zero in sprint-time !! QUESTION: What does this mean? HCW
750
751 IF (delta_gdd < 0) THEN !GDD cannot be negative
752 indhelp = delta_gdd !Amount of negative GDD
753 delta_gdd = 0
754 END IF
755
756 IF (delta_sdd > 0) delta_sdd = 0 !SDD cannot be positive
757
758 ! Calculate cumulative growing and senescence degree days
759 gdd_id(iv) = gdd_id_prev(iv) + delta_gdd
760 sdd_id(iv) = sdd_id_prev(iv) + delta_sdd
761
762 ! Possibility for cold spring
763 IF (sdd_id(iv) <= sddfull(iv) .AND. indhelp < 0) THEN
764 gdd_id(iv) = 0
765 END IF
766
767 IF (gdd_id(iv) >= gddfull(iv)) THEN !Start senescence
768 gdd_id(iv) = gddfull(iv) !Leaves should not grow so delete yes from earlier
769 IF (sdd_id(iv) < -critdays) gdd_id(iv) = 0
770 END IF
771
772 IF (sdd_id(iv) <= sddfull(iv)) THEN !After senescence now start growing leaves
773 sdd_id(iv) = sddfull(iv) !Leaves off so add back earlier
774 IF (gdd_id(iv) > critdays) sdd_id(iv) = 0
775 END IF
776
777 ! With these limits SDD, GDD is set to zero
778 IF (sdd_id(iv) < -critdays .AND. sdd_id(iv) > sddfull(iv)) gdd_id(iv) = 0
779 IF (gdd_id(iv) > critdays .AND. gdd_id(iv) < gddfull(iv)) sdd_id(iv) = 0
780
781 ! Now calculate LAI itself
782 IF (lat >= 0) THEN !Northern hemispere
783 !If SDD is not zero by mid May, this is forced
784 IF (id == 140 .AND. sdd_id(iv) /= 0) sdd_id(iv) = 0
785 ! Set SDD to zero in summer time
786 IF (gdd_id(iv) > critdays .AND. id < 170) sdd_id(iv) = 0
787 ! Set GDD zero in winter time
788 IF (sdd_id(iv) < -critdays .AND. id > 170) gdd_id(iv) = 0
789
790 IF (laitype(iv) < 0.5) THEN !Original LAI type
791 IF (gdd_id(iv) > 0 .AND. gdd_id(iv) < gddfull(iv)) THEN !Leaves can still grow
792 lai_id_next(iv) = (lai_id_prev(iv)**laipower(1, iv)*gdd_id(iv)*laipower(2, iv)) + lai_id_prev(iv)
793 ELSEIF (sdd_id(iv) < 0 .AND. sdd_id(iv) > sddfull(iv)) THEN !Start senescence
794 lai_id_next(iv) = (lai_id_prev(iv)**laipower(3, iv)*sdd_id(iv)*laipower(4, iv)) + lai_id_prev(iv)
795 ELSE
796 lai_id_next(iv) = lai_id_prev(iv)
797 END IF
798 ELSEIF (laitype(iv) >= 0.5) THEN
799 IF (gdd_id(iv) > 0 .AND. gdd_id(iv) < gddfull(iv)) THEN !Leaves can still grow
800 lai_id_next(iv) = (lai_id_prev(iv)**laipower(1, iv)*gdd_id(iv)*laipower(2, iv)) + lai_id_prev(iv)
801 !! Use day length to start senescence at high latitudes (N hemisphere)
802 ELSEIF (lenday_id_prev <= 12 .AND. sdd_id(iv) > sddfull(iv)) THEN !Start senescence
803 lai_id_next(iv) = (lai_id_prev(iv)*laipower(3, iv)*(1 - sdd_id(iv))*laipower(4, iv)) + lai_id_prev(iv)
804 ELSE
805 lai_id_next(iv) = lai_id_prev(iv)
806 END IF
807 END IF
808
809 ELSEIF (lat < 0) THEN !Southern hemisphere !! N.B. not identical to N hemisphere - return to later
810 !If SDD is not zero by late Oct, this is forced
811 IF (id == 300 .AND. sdd_id(iv) /= 0) sdd_id(iv) = 0
812 ! Set SDD to zero in summer time
813 IF (gdd_id(iv) > critdays .AND. id > 250) sdd_id(iv) = 0
814 ! Set GDD zero in winter time
815 IF (sdd_id(iv) < -critdays .AND. id < 250) gdd_id(iv) = 0
816
817 IF (laitype(iv) < 0.5) THEN !Original LAI type
818 IF (gdd_id(iv) > 0 .AND. gdd_id(iv) < gddfull(iv)) THEN
819 lai_id_next(iv) = (lai_id_prev(iv)**laipower(1, iv)*gdd_id(iv)*laipower(2, iv)) + lai_id_prev(iv)
820 ELSEIF (sdd_id(iv) < 0 .AND. sdd_id(iv) > sddfull(iv)) THEN
821 lai_id_next(iv) = (lai_id_prev(iv)**laipower(3, iv)*sdd_id(iv)*laipower(4, iv)) + lai_id_prev(iv)
822 ELSE
823 lai_id_next(iv) = lai_id_prev(iv)
824 END IF
825 ELSE
826 IF (gdd_id(iv) > 0 .AND. gdd_id(iv) < gddfull(iv)) THEN
827 lai_id_next(iv) = (lai_id_prev(iv)**laipower(1, iv)*gdd_id(iv)*laipower(2, iv)) + lai_id_prev(iv)
828 !! Day length not used to start senescence in S hemisphere (not much land)
829 ELSEIF (sdd_id(iv) < 0 .AND. sdd_id(iv) > sddfull(iv)) THEN
830 lai_id_next(iv) = (lai_id_prev(iv)*laipower(3, iv)*(1 - sdd_id(iv))*laipower(4, iv)) + lai_id_prev(iv)
831 ELSE
832 lai_id_next(iv) = lai_id_prev(iv)
833 END IF
834 END IF
835 END IF !N or S hemisphere
836
837 ! Check LAI within limits; if not set to limiting value
838 IF (lai_id_next(iv) > laimax(iv)) THEN
839 lai_id_next(iv) = laimax(iv)
840 ELSEIF (lai_id_next(iv) < laimin(iv)) THEN
841 lai_id_next(iv) = laimin(iv)
842 END IF
843
844 END DO !End of loop over veg surfaces
845
846 IF (laicalcyes == 0) THEN ! moved to SUEWS_cal_DailyState, TS 18 Sep 2017
847 ! LAI(id-1,:)=LAI_obs ! check -- this is going to be a problem as it is not for each vegetation class
848 lai_id_next = lai_obs
849 END IF
850 !------------------------------------------------------------------------------
851
852 END SUBROUTINE update_gddlai
853
854 SUBROUTINE update_wateruse( &
855 id, WaterUseMethod, DayofWeek_id, lat, FrIrriAuto, HDD_id, & !input
856 state_id, soilstore_id, SoilStoreCap, H_maintain, & !input
857 Ie_a, Ie_m, Ie_start, Ie_end, DayWatPer, DayWat, &
858 WUDay_id) !output
859
860 IMPLICIT NONE
861
862 INTEGER, INTENT(IN) :: id
863 INTEGER, INTENT(IN) :: WaterUseMethod
864 INTEGER, INTENT(IN) :: Ie_start !Starting time of water use (DOY)
865 INTEGER, INTENT(IN) :: Ie_end !Ending time of water use (DOY)
866 INTEGER, DIMENSION(3), INTENT(IN) :: DayofWeek_id
867
868 REAL(KIND(1D0)), INTENT(IN) :: lat
869 REAL(KIND(1D0)), INTENT(IN) :: FrIrriAuto !Fraction of irrigated area using automatic irrigation
870
871 REAL(KIND(1D0)), DIMENSION(12), INTENT(IN) :: HDD_id
872 REAL(KIND(1D0)), DIMENSION(NVegSurf), INTENT(IN) :: Ie_a !Coefficients for automatic irrigation models
873 REAL(KIND(1D0)), DIMENSION(NVegSurf), INTENT(IN) :: Ie_m !Coefficients for manual irrigation models
874 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: DayWatPer !% of houses following daily water
875 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: DayWat !Days of watering allowed
876
877 ! ponding control related
878 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: state_id ! surface wetness [mm]
879 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(IN) :: soilstore_id ! soil water store [mm]
880 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap !Capacity of soil store for each surface [mm]
881 REAL(KIND(1D0)), INTENT(IN) :: H_maintain ! ponding water depth to maintain [mm]
882
883 REAL(KIND(1D0)), DIMENSION(9), INTENT(OUT) :: WUDay_id !Daily water use for EveTr, DecTr, Grass [mm] (see SUEWS_DailyState.f95)
884
885 REAL(KIND(1D0)), DIMENSION(3) :: h_need !water level to maintain: surface+soil [mm]
886 REAL(KIND(1D0)), DIMENSION(3) :: store_total !current water level: surface+soil [mm]
887 REAL(KIND(1D0)), DIMENSION(3) :: WUDay_P !water used to maintain ponding level [mm]
888 REAL(KIND(1D0)), DIMENSION(3) :: WUDay_A !automatic irrigation [mm]
889 REAL(KIND(1D0)), DIMENSION(3) :: WUDay_M !manual irrigation [mm]
890 REAL(KIND(1D0)), DIMENSION(3) :: WUDay_total !Coefficients for manual irrigation models
891
892 INTEGER :: wd !Water use calculation is done when calc = 1
893 INTEGER :: calc !Water use calculation is done when calc = 1
894 INTEGER :: i
895
896 REAL(KIND(1D0)) :: temp_avg
897 REAL(KIND(1D0)) :: days_since_rain
898
899 ! transfer HDD values
900 temp_avg = hdd_id(9)
901 days_since_rain = hdd_id(12)
902
903 ! initialise WUDay_id
904 wuday_id = 0
905
906 IF (waterusemethod == 0) THEN !If water use is to be modelled (rather than observed)
907
908 wd = dayofweek_id(1)
909
910 IF (daywat(wd) == 1.0) THEN !1 indicates watering permitted on this day
911 calc = 0
912 IF (lat >= 0) THEN !Northern Hemisphere
913 IF (id >= ie_start - 1 .AND. id <= ie_end + 1) calc = 1 !Day between irrigation period
914 ELSE !Southern Hemisphere
915 calc = 1
916 IF (id >= ie_end .AND. id <= ie_start) calc = 0 !Day between irrigation period
917 END IF
918
919 IF (calc == 1) THEN
920 ! Model daily water use based on days_since_rain (days since rain) and temp_avg (average temp)
921 ! WUDay is the amount of water [mm] per day, applied to each of the irrigated areas
922 ! N.B. These are the same for each vegetation type at the moment
923
924 ! ---- irrigation amount to maintain a certain water availability----
925 ! NB: H_maintain can be either positive or negative
926 h_need = soilstorecap(3:5) + h_maintain
927 store_total = state_id(3:5) + soilstore_id(3:5)
928 wuday_p = h_need - store_total
929 wuday_p = merge(wuday_p, 0d0, wuday_p > 0)
930
931 ! ---- automatic irrigation ----
932 wuday_a = frirriauto*(ie_a(1) + ie_a(2)*temp_avg + ie_a(3)*days_since_rain)*daywatper(wd)
933 wuday_a = merge(wuday_a, 0d0, wuday_a > 0)
934 ! add ponding-demand to auto-irrigation
935 wuday_a = wuday_a + wuday_p
936
937 ! ---- Manual irrigation----
938 wuday_m = (1 - frirriauto)*(ie_m(1) + ie_m(2)*temp_avg + ie_m(3)*days_since_rain)*daywatper(wd)
939 wuday_m = merge(wuday_m, 0d0, wuday_m > 0)
940
941 ! ---- total irrigation
942 wuday_total = wuday_p + wuday_a + wuday_m
943
944 ! transfer values to WUDay_id
945 wuday_id([((i - 1)*3 + 1, i=1, 3)]) = wuday_total
946 wuday_id([((i - 1)*3 + 2, i=1, 3)]) = wuday_a
947 wuday_id([((i - 1)*3 + 3, i=1, 3)]) = wuday_m
948
949 ELSE !If no irrigation on this day
950 wuday_id = 0
951 END IF
952 END IF
953 END IF
954
955 END SUBROUTINE update_wateruse
956
957 SUBROUTINE update_hdd( &
958 dt_since_start, it, imin, tstep, & !input
959 HDD_id) !inout
960 IMPLICIT NONE
961 INTEGER, INTENT(IN) :: dt_since_start, it, imin, tstep
962
963 REAL(KIND(1D0)), DIMENSION(12), INTENT(INOUT) :: HDD_id
964 ! REAL(KIND(1d0)),DIMENSION(6),INTENT(OUT):: HDD_id_use
965
966 INTEGER :: days_prev
967 REAL(KIND(1D0)) :: tstepcount
968
969 ! count of timesteps performed during day
970 tstepcount = (it*60 + imin)*60/tstep*1.
971 ! Heating degree days (HDD) -------------
972 hdd_id(1) = hdd_id(1)/tstepcount !Heating
973 hdd_id(2) = hdd_id(2)/tstepcount !Cooling
974 hdd_id(3) = hdd_id(3)/tstepcount !Average temp
975
976 ! Calculate a quasi-5-day-running-mean temp
977 days_prev = min(4, & ! dt_since_start >= 4 days
978 floor(dt_since_start/(24*60*60)*1.)) ! dt_since_start < 4 days
979 hdd_id(4) = (hdd_id(4)*days_prev + hdd_id(3))/(days_prev + 1)
980
981 ! Calculate number of days since rain
982 IF (hdd_id(5) > 0) THEN !Rain occurred
983 hdd_id(6) = 0
984 ELSE
985 hdd_id(6) = hdd_id(6) + 1 !Days since rain
986 END IF
987
988 ! save updated HDD_id(1:6) values to the last-half part (i.e., HDD_id(7:12))
989 hdd_id(6 + 1:6 + 6) = hdd_id(1:6)
990
991 END SUBROUTINE update_hdd
992
994 it, imin, & !input
995 HDD_id) !output
996 IMPLICIT NONE
997 INTEGER, INTENT(IN) :: it
998 INTEGER, INTENT(IN) :: imin
999
1000 REAL(KIND(1D0)), DIMENSION(6), INTENT(INOUT) :: HDD_id
1001 REAL(KIND(1D0)) :: HDD_id_mav, HDD_id_daysSR
1002
1003 ! reset HDD_id to ZERO except for:
1004 ! 5-day moving average
1005 hdd_id_mav = hdd_id(4)
1006 ! Days Since Rain
1007 hdd_id_dayssr = hdd_id(6)
1008 IF (it == 0 .AND. imin == 0) THEN
1009 hdd_id = 0
1010 hdd_id(4) = hdd_id_mav
1011 hdd_id(6) = hdd_id_dayssr
1012 END IF
1013
1014 END SUBROUTINE update_dailystate_start
1015
1017 id, datetimeline, & !input
1018 Gridiv, NumberOfGrids, &
1019 DailyStateLine, &
1020 dataOutDailyState) !inout
1021
1022 IMPLICIT NONE
1023
1024 ! INTEGER,INTENT(IN) ::iy
1025 INTEGER, INTENT(IN) :: id
1026 ! INTEGER,INTENT(IN) ::it
1027 ! INTEGER,INTENT(IN) ::imin
1028
1029 REAL(KIND(1D0)), DIMENSION(5), INTENT(IN) :: datetimeline
1030
1031 INTEGER, INTENT(IN) :: Gridiv
1032 INTEGER, INTENT(IN) :: NumberOfGrids
1033 REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutDailyState - 5), INTENT(IN) :: DailyStateLine
1034 REAL(KIND(1D0)), DIMENSION(ndays, ncolumnsDataOutDailyState, NumberOfGrids), INTENT(INOUT) :: dataOutDailyState
1035
1036 ! write out to dataOutDailyState
1037 dataoutdailystate(id, 1:5, gridiv) = datetimeline
1038 ! DailyStateLine will be -999 unless realistic values are calculated at the last timestep of each day
1039 dataoutdailystate(id, 6:ncolumnsdataoutdailystate, gridiv) = dailystateline
1040
1041 END SUBROUTINE suews_update_dailystate
1042
1043 ! transfer results to a one-line output for SUEWS_cal_DailyState
1045 it, imin, nsh_real, & !input
1046 GDD_id, HDD_id, LAI_id, &
1047 SDD_id, &
1048 Tmin_id, Tmax_id, lenday_id, &
1049 DecidCap_id, &
1050 albDecTr_id, &
1051 albEveTr_id, &
1052 albGrass_id, &
1053 porosity_id, &
1054 WUDay_id, &
1055 deltaLAI, VegPhenLumps, &
1056 SnowAlb, SnowDens, &
1057 a1, a2, a3, &
1058 DailyStateLine) !out
1059
1060 IMPLICIT NONE
1061
1062 ! INTEGER,INTENT(IN) ::iy
1063 ! INTEGER,INTENT(IN) ::id
1064 INTEGER, INTENT(IN) :: it
1065 INTEGER, INTENT(IN) :: imin
1066 REAL(KIND(1D0)), INTENT(IN) :: nsh_real
1067
1068 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: GDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
1069 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: SDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
1070 REAL(KIND(1D0)), DIMENSION(6), INTENT(IN) :: HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
1071 REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAI_id !LAI for each veg surface [m2 m-2]
1072
1073 REAL(KIND(1D0)), INTENT(IN) :: DecidCap_id
1074 REAL(KIND(1D0)), INTENT(IN) :: albDecTr_id
1075 REAL(KIND(1D0)), INTENT(IN) :: albEveTr_id
1076 REAL(KIND(1D0)), INTENT(IN) :: albGrass_id
1077 REAL(KIND(1D0)), INTENT(IN) :: porosity_id
1078 REAL(KIND(1D0)), INTENT(IN) :: Tmin_id
1079 REAL(KIND(1D0)), INTENT(IN) :: Tmax_id
1080 REAL(KIND(1D0)), INTENT(IN) :: lenday_id
1081 REAL(KIND(1D0)), DIMENSION(9), INTENT(IN) :: WUDay_id !Daily water use for EveTr, DecTr, Grass [mm] (see SUEWS_DailyState.f95)
1082
1083 REAL(KIND(1D0)), INTENT(IN) :: deltaLAI
1084 REAL(KIND(1D0)), INTENT(IN) :: VegPhenLumps
1085 REAL(KIND(1D0)), INTENT(IN) :: SnowAlb
1086 REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) :: SnowDens
1087 REAL(KIND(1D0)), INTENT(IN) :: a1
1088 REAL(KIND(1D0)), INTENT(IN) :: a2
1089 REAL(KIND(1D0)), INTENT(IN) :: a3
1090
1091 REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutDailyState - 5), INTENT(OUT) :: DailyStateLine
1092
1093 ! initialise DailyStateLine
1094 dailystateline = -999
1095 IF (it == 23 .AND. imin == (nsh_real - 1)/nsh_real*60) THEN
1096 ! Write actual data only at the last timesstep of each day
1097 ! DailyStateLine(1:2) = [iy,id]
1098 ! DailyStateLine(1:6) = HDD_id
1099 ! DailyStateLine(6 + 1:6 + 5) = GDD_id
1100 ! DailyStateLine(11 + 1:11 + 3) = LAI_id
1101 ! DailyStateLine(14 + 1:14 + 5) = [DecidCap_id, Porosity_id, AlbEveTr_id, AlbDecTr_id, AlbGrass_id]
1102 ! DailyStateLine(19 + 1:19 + 9) = WUDay_id(1:9)
1103 ! DailyStateLine(28 + 1) = deltaLAI
1104 ! DailyStateLine(29 + 1) = VegPhenLumps
1105 ! DailyStateLine(30 + 1:30 + 8) = [SnowAlb, SnowDens(1:7)]
1106 ! DailyStateLine(38 + 1:38 + 3) = [a1, a2, a3]
1107 dailystateline = [hdd_id, gdd_id, sdd_id, tmin_id, tmax_id, lenday_id, lai_id, decidcap_id, porosity_id, &
1108 albevetr_id, albdectr_id, albgrass_id, wuday_id, deltalai, vegphenlumps, snowalb, snowdens, &
1109 a1, a2, a3]
1110
1111 END IF
1112
1113 END SUBROUTINE update_dailystateline
1114
1115END MODULE dailystate_module
integer, parameter ivgrass
integer, parameter nvegsurf
integer, parameter ivdecid
integer, parameter ncolumnsdataoutdailystate
integer, parameter nsurf
integer, parameter ndays
integer, parameter ivconif
integer, parameter decidsurf
subroutine update_gddlai(id, LAICalcYes, lat, LAI_obs, Tmin_id_prev, Tmax_id_prev, lenDay_id_prev, BaseT, BaseTe, GDDFull, SDDFull, LAIMin, LAIMax, LAIPower, LAIType, LAI_id_prev, GDD_id, SDD_id, LAI_id_next)
subroutine suews_cal_dailystate(iy, id, it, imin, isec, tstep, tstep_prev, dt_since_start, DayofWeek_id, Tmin_id_prev, Tmax_id_prev, lenDay_id_prev, BaseTMethod, WaterUseMethod, Ie_start, Ie_end, LAICalcYes, LAIType, nsh_real, avkdn, Temp_C, Precip, BaseT_HC, BaseT_Heating, BaseT_Cooling, lat, Faut, LAI_obs, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, CapMax_dec, CapMin_dec, PorMax_dec, PorMin_dec, Ie_a, Ie_m, DayWatPer, DayWat, BaseT, BaseTe, GDDFull, SDDFull, LAIMin, LAIMax, LAIPower, DecidCap_id_prev, StoreDrainPrm_prev, LAI_id_prev, GDD_id_prev, SDD_id_prev, albDecTr_id_prev, albEveTr_id_prev, albGrass_id_prev, porosity_id_prev, HDD_id_prev, state_id, soilstore_id, SoilStoreCap, H_maintain, HDD_id_next, Tmin_id_next, Tmax_id_next, lenDay_id_next, albDecTr_id_next, albEveTr_id_next, albGrass_id_next, porosity_id_next, DecidCap_id_next, StoreDrainPrm_next, LAI_id_next, GDD_id_next, SDD_id_next, deltaLAI, WUDay_id)
subroutine update_hdd(dt_since_start, it, imin, tstep, HDD_id)
subroutine update_dailystate_day(BaseTMethod, DayofWeek_id, avkdn, Temp_C, Precip, BaseT_HC, BaseT_Heating, BaseT_Cooling, nsh_real, Tmin_id, Tmax_id, lenDay_id, HDD_id)
subroutine update_dailystateline(it, imin, nsh_real, GDD_id, HDD_id, LAI_id, SDD_id, Tmin_id, Tmax_id, lenday_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, WUDay_id, deltaLAI, VegPhenLumps, SnowAlb, SnowDens, a1, a2, a3, DailyStateLine)
subroutine update_wateruse(id, WaterUseMethod, DayofWeek_id, lat, FrIrriAuto, HDD_id, state_id, soilstore_id, SoilStoreCap, H_maintain, Ie_a, Ie_m, Ie_start, Ie_end, DayWatPer, DayWat, WUDay_id)
subroutine update_dailystate_start(it, imin, HDD_id)
subroutine update_veg(LAImax, LAIMin, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, CapMax_dec, CapMin_dec, PorMax_dec, PorMin_dec, LAI_id, LAI_id_prev, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, StoreDrainPrm, deltaLAI)
subroutine update_dailystate_end(id, it, imin, tstep, dt_since_start, Tmin_id, Tmax_id, lenDay_id, LAIType, Ie_end, Ie_start, LAICalcYes, WaterUseMethod, DayofWeek_id, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, BaseT, BaseTe, CapMax_dec, CapMin_dec, DayWat, DayWatPer, Faut, GDDFull, Ie_a, Ie_m, LAIMax, LAIMin, LAIPower, lat, PorMax_dec, PorMin_dec, SDDFull, LAI_obs, state_id, soilstore_id, SoilStoreCap, H_maintain, GDD_id, SDD_id, HDD_id, LAI_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, StoreDrainPrm, WUDay_id, deltaLAI)
subroutine suews_update_dailystate(id, datetimeline, Gridiv, NumberOfGrids, DailyStateLine, dataOutDailyState)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)