SUEWS API Site
Documentation of SUEWS source code
suews_phys_resist.f95
Go to the documentation of this file.
2  implicit none
3 
4 contains
5  SUBROUTINE aerodynamicresistance( &
6  ZZD, &! input:
7  z0m, &
8  AVU1, &
9  L_mod, &
10  UStar, &
11  VegFraction, &
12  AerodynamicResistanceMethod, &
13  StabilityMethod, &
14  RoughLenHeatMethod, &
15  RA)! output:
16 
17  ! Returns Aerodynamic resistance (RA) to the main program SUEWS_Calculations
18  ! All RA equations reported in Thom & Oliver (1977)
19  ! Modified by TS 08 Aug 2017 - interface modified
20  ! Modified by LJ
21  ! -Removal of tabs and cleaning the code
22  ! Modified by HCW 03 Dec 2015 - changed lower limit on RA from 2 s m-1 to 10 s m-1 (to avoid unrealistically high evaporation rates)
23  ! Modified by LJ in 12 April to not to be used with global variables
24  ! To Do:
25  ! - Check whether the thresholds 2-200 s m-1 are suitable over a range of z0!! HCW 04 Mar 2015
26  ! OUTPUT: RA - Aerodynamic resistance [s m^-1]
27  ! INPUT: AerodynamicResistanceMethod = Method to calculate RA
28  ! StabilityMethod = defines the method to calculate atmospheric stability
29  ! RoughLenHeatMethod = Method to calculate heat roughness length
30  ! *Measurement height minus* Displacement height (m) (was incorrectly labelled, corrected HCW 25 May 2016
31  ! z0m = Aerodynamic roughness length (m)
32  ! k2 = Power of Van Karman's constant (= 0.16 = 0.4^2)
33  ! AVU1 = Mean wind speed
34  ! L_mod = Obukhov length (m)
35  ! UStar = Friction velocity (m s-1)
36  ! VegFraction = Fraction of vegetation
37  ! (changed from veg_fr which also includes water surface by HCW 05 Nov 2015)
38 
40  use sues_data, only: psih
41 
42  IMPLICIT NONE
43 
44  REAL(KIND(1d0)), INTENT(in)::ZZD!Active measurement height (meas. height-displac. height)
45  REAL(KIND(1d0)), INTENT(in)::z0m!Aerodynamic roughness length
46  REAL(KIND(1d0)), INTENT(in)::AVU1!Average wind speed
47  REAL(KIND(1d0)), INTENT(in)::L_mod!Monin-Obukhov length (either measured or modelled)
48  REAL(KIND(1d0)), INTENT(in)::UStar!Friction velocity
49  REAL(KIND(1d0)), INTENT(in)::VegFraction!Fraction of vegetation
50 
51  INTEGER, INTENT(in)::AerodynamicResistanceMethod
52  INTEGER, INTENT(in)::StabilityMethod
53  INTEGER, INTENT(in)::RoughLenHeatMethod
54 
55  REAL(KIND(1d0)), INTENT(out)::RA !Aerodynamic resistance [s m^-1]
56 
57  INTEGER, PARAMETER :: notUsedI = -55
58 
59  REAL(KIND(1d0)), PARAMETER :: &
60  notUsed = -55.5, &
61  k2 = 0.16, & !Power of Van Karman's constant (= 0.16 = 0.4^2)
62  muu = 1.46e-5 !molecular viscosity
63  REAL(KIND(1d0)):: psim
64  ! REAL(KIND(1d0)):: psih
65  REAL(KIND(1d0)):: z0V
66 
67  !1)Monteith (1965)-neutral stability
68  IF (aerodynamicresistancemethod == 1) THEN
69  ra = (log(zzd/z0m)**2)/(k2*avu1)
70 
71  !2) Non-neutral stability
72  ! PSIM - stability function for momentum
73  ! psih - stability function for heat
74  ! assuming stability functions the same for heat and water
75  ELSEIF (aerodynamicresistancemethod == 2) THEN !Dyer (1974)
76 
77  psim = stab_psi_mom(stabilitymethod, zzd/l_mod)
78  psih = stab_psi_heat(stabilitymethod, zzd/l_mod)
79 
80  !Z0V roughness length for vapour
81  z0v = cal_z0v(roughlenheatmethod, z0m, vegfraction, ustar)
82 
83  IF (zzd/l_mod == 0 .OR. ustar == 0) THEN
84  ra = (log(zzd/z0m)*log(zzd/z0v))/(k2*avu1) !Use neutral equation
85  ELSE
86  ra = ((log(zzd/z0m) - psim)*(log(zzd/z0v) - psih))/(k2*avu1)
87  ENDIF
88 
89  !3) Thom and Oliver (1977)
90  ELSEIF (aerodynamicresistancemethod == 3) THEN
91  ra = (4.72*log(zzd/z0m)**2)/(1 + 0.54*avu1)
92  ENDIF
93 
94  !If RA outside permitted range, adjust extreme values !!Check whether these thresholds are suitable over a range of z0
95  IF (ra > 200) THEN !was 175
96  CALL errorhint(7, 'In AerodynamicResistance.f95, calculated RA > 200 s m-1; RA set to 200 s m-1', ra, notused, notusedi)
97  ra = 200
98  ELSEIF (ra < 10) THEN !found By Shiho - fix Dec 2012 !Threshold changed from 2 to 10 s m-1 (HCW 03 Dec 2015)
99  CALL errorhint(7, 'In AerodynamicResistance.f95, calculated RA < 10 s m-1; RA set to 10 s m-1', ra, notused, notusedi)
100  ra = 10
101  ! RA=(log(ZZD/z0m))**2/(k2*AVU1)
102  IF (avu1 < 0) WRITE (*, *) avu1, ra
103  ENDIF
104 
105  RETURN
106  END SUBROUTINE aerodynamicresistance
107 
108  SUBROUTINE surfaceresistance( &
109  id, it, &! input:
110  SMDMethod, SnowFrac, sfr, avkdn, Temp_C, dq, xsmd, vsmd, MaxConductance, &
111  LAIMax, LAI_id, gsModel, Kmax, &
112  G1, G2, G3, G4, G5, G6, TH, TL, S1, S2, &
113  gfunc, gsc, ResistSurf)! output:
114  ! Calculates bulk surface resistance (ResistSurf [s m-1]) based on Jarvis 1976 approach
115  ! Last modified -----------------------------------------------------
116  ! MH 01 Feb 2019: gsModel choices to model with air temperature or 2 meter temperature. Added gfunc for photosynthesis calculations
117  ! HCW 21 Jul 2016: If no veg surfaces, vsmd = NaN so QE & QH = NaN; if water surfaces only, smd = NaN so QE & QH = NaN.
118  ! Add checks here so that gs (soil part) = 0 in either of these situations.
119  ! This shouldn't change results but handles NaN error.
120  ! HCW 01 Mar 2016: SM dependence is now on modelled smd for vegetated surfaces only (vsmd) (Note: obs smd still not operational!)
121  ! HCW 18 Jun 2015: Alternative gs parameterisation added using different functional forms and new coefficients
122  ! HCW 31 Jul 2014: Modified condition on g6 part to select meas/mod smd
123  ! LJ 24 Apr 2013: Added impact of snow fraction in LAI and in soil moisture deficit
124  ! -------------------------------------------------------------------
125 
126  ! USE allocateArray
127  ! USE data_in
128  ! USE defaultNotUsed
129  ! USE gis_data
130  ! USE moist
131  ! USE resist
132  ! USE sues_data
133 
134  IMPLICIT NONE
135  ! INTEGER,PARAMETER::BldgSurf=2
136  INTEGER, PARAMETER::ConifSurf = 3
137  INTEGER, PARAMETER::DecidSurf = 4
138  INTEGER, PARAMETER::GrassSurf = 5
139  ! INTEGER,PARAMETER::ivConif=1
140  ! INTEGER,PARAMETER::ivGrass=3
141  ! INTEGER,PARAMETER::MaxNumberOfGrids=2000
142  ! INTEGER,PARAMETER::ndays=366
143  INTEGER, PARAMETER::nsurf = 7
144  ! INTEGER,PARAMETER::NVegSurf=3
145  ! INTEGER,PARAMETER::PavSurf=1
146  INTEGER, PARAMETER::WaterSurf = 7
147 
148  INTEGER, INTENT(in)::id
149  INTEGER, INTENT(in)::it ! time: day of year and hour
150  INTEGER, INTENT(in)::gsModel!Choice of gs parameterisation (1 = Ja11, 2 = Wa16)
151  INTEGER, INTENT(in)::SMDMethod!Method of measured soil moisture
152  ! INTEGER,INTENT(in)::ConifSurf!= 3, surface code
153  ! INTEGER,INTENT(in)::DecidSurf!= 4, surface code
154  ! INTEGER,INTENT(in)::GrassSurf!= 5, surface code
155  ! INTEGER,INTENT(in)::WaterSurf!= 7, surface code
156  ! INTEGER,INTENT(in)::nsurf!= 7, Total number of surfaces
157 
158  REAL(KIND(1d0)), INTENT(in)::avkdn!Average downwelling shortwave radiation
159  REAL(KIND(1d0)), INTENT(in)::Temp_C!Air temperature
160  REAL(KIND(1d0)), INTENT(in)::Kmax!Annual maximum hourly solar radiation
161  REAL(KIND(1d0)), INTENT(in)::G1!Fitted parameters related to surface res. calculations
162  REAL(KIND(1d0)), INTENT(in)::G2!Fitted parameters related to surface res. calculations
163  REAL(KIND(1d0)), INTENT(in)::G3!Fitted parameters related to surface res. calculations
164  REAL(KIND(1d0)), INTENT(in)::G4!Fitted parameters related to surface res. calculations
165  REAL(KIND(1d0)), INTENT(in)::G5!Fitted parameters related to surface res. calculations
166  REAL(KIND(1d0)), INTENT(in)::G6!Fitted parameters related to surface res. calculations
167  REAL(KIND(1d0)), INTENT(in)::S1!Fitted parameters related to surface res. calculations
168  REAL(KIND(1d0)), INTENT(in)::S2!Fitted parameters related to surface res. calculations
169  REAL(KIND(1d0)), INTENT(in)::TH!Maximum temperature limit
170  REAL(KIND(1d0)), INTENT(in)::TL!Minimum temperature limit
171  REAL(KIND(1d0)), INTENT(in)::dq!Specific humidity deficit
172  REAL(KIND(1d0)), INTENT(in)::xsmd!Measured soil moisture deficit
173  REAL(KIND(1d0)), INTENT(in)::vsmd!QUESTION: Soil moisture deficit for vegetated surfaces only (what about BSoil?)
174 
175  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::MaxConductance!Max conductance [mm s-1]
176  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::LAIMax !Max LAI [m2 m-2]
177  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::LAI_id !=LAI(id-1,:), LAI for each veg surface [m2 m-2]
178  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac !Surface fraction of snow cover
179  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr !Surface fractions [-]
180 
181  REAL(KIND(1d0)), INTENT(out)::gfunc!gdq*gtemp*gs*gq for photosynthesis calculations
182  REAL(KIND(1d0)), INTENT(out)::gsc!Surface Layer Conductance
183  REAL(KIND(1d0)), INTENT(out)::ResistSurf!Surface resistance
184 
185  REAL(KIND(1d0)):: &
186  gl, &!G(LAI)
187  QNM, &!QMAX/(QMAX+G2)
188  gq, &!G(Q*)
189  gdq, &!G(dq)
190  TC, &!Temperature parameter 1
191  TC2, &!Temperature parameter 2
192  gtemp, &!G(T)
193  sdp, &!S1/G6+S2
194  gs!G(Soil moisture deficit)
195 
196  INTEGER:: iv
197  REAL(KIND(1d0)):: id_real
198 
199  REAL(KIND(1d0)), PARAMETER :: notUsed = -55
200  ! REAL(KIND(1d0)),PARAMETER :: notUsedi=-55.5
201 
202  ! initialisation
203  gdq = 0.5
204  gtemp = 0.5
205  gs = 0.5
206  gq = 0.5
207 
208  id_real = REAL(id) !Day of year in real number
209 
210  !gsModel = 1 - original parameterisation (Jarvi et al. 2011)
211  !gsModel = 2 - new parameterisation (Ward et al. 2016)
212  !gsModel = 3 - original parameterisation (Jarvi et al. 2011) with 2 m temperature
213  !gsModel = 4 - new parameterisation (Ward et al. 2016) with 2 m temperature
214 
215  IF (gsmodel == 1 .OR. gsmodel == 3) THEN
216  ! kdown ----
217  qnm = kmax/(kmax + g2)
218  !gq=(qn1/(g2+qn1))/qnm !With net all-wave radiation
219  gq = (avkdn/(g2 + avkdn))/qnm !With Kdown
220  ! specific humidity deficit ----
221  IF (dq < g4) THEN
222  gdq = 1 - g3*dq
223  ELSE
224  gdq = 1 - g3*g4
225  ENDIF
226  ! air temperature ----
227  tc = (th - g5)/(g5 - tl)
228  tc2 = (g5 - tl)*(th - g5)**tc
229  !If air temperature below TL or above TH, fit it to TL+0.1/TH-0.1
230  IF (temp_c <= tl) THEN
231  gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
232  !Call error only if no snow on ground
233  ! IF (MIN(SnowFrac(1),SnowFrac(2),SnowFrac(3),SnowFrac(4),SnowFrac(5),SnowFrac(6))/=1) THEN
234  IF (minval(snowfrac(1:6)) /= 1) THEN
235  CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TL=0.1,Temp_c,id,it', &
236  REAL(Temp_c, KIND(1d0)), id_real, it)
237  ENDIF
238  ELSEIF (temp_c >= th) THEN
239  gtemp = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
240  CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TH=39.9,Temp_c,id,it', &
241  REAL(Temp_c, KIND(1d0)), id_real, it)
242  ELSE
243  gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
244  ENDIF
245 
246  ! soil moisture deficit ----
247  sdp = s1/g6 + s2
248  IF (smdmethod > 0) THEN !Modified from ==1 to > 0 by HCW 31/07/2014
249  gs = 1 - exp(g6*(xsmd - sdp)) !Measured soil moisture deficit is used
250  ELSE
251  gs = 1 - exp(g6*(vsmd - sdp)) !Modelled is used
252  IF (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) == 0 .OR. sfr(watersurf) == 1) THEN
253  gs = 0 !If no veg so no vsmd, or all water so no smd, set gs=0 (HCW 21 Jul 2016)
254  ENDIF
255  ENDIF
256  !gs = gs*(1 - SUM(SnowFrac(1:6))/6)
257  IF (gs < 0) THEN
258  CALL errorhint(65, &
259  'subroutine SurfaceResistance.f95 (gsModel=1): g(smd) < 0 calculated, setting to 0.0001', &
260  gs, id_real, it)
261  gs = 0.0001
262  ENDIF
263 
264  !LAI
265  !Original way
266  !gl=((LAI(id,2)*areaunir/lm)+areair)/(areair+areaunir)
267  !New way
268  gl = 0 !First initialize
269  ! vegetated surfaces
270  ! check basis for values koe - maximum surface conductance
271  ! print*,id,it,sfr
272  ! DO iv=ivConif,ivGrass
273  DO iv = 1, 3
274  ! gl=gl+(sfr(iv+2)*(1-SnowFrac(iv+2)))*LAI(id-1,iv)/LAIMax(iv)*MaxConductance(iv)
275  gl = gl + (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
276  ENDDO
277 
278  IF (avkdn <= 0) THEN !At nighttime set gsc at arbitrary low value: gsc=0.1 mm/s (Shuttleworth, 1988b)
279  gsc = 0.1
280  ELSE
281  ! Multiply parts together
282  gsc = (g1*gq*gdq*gtemp*gs*gl)
283  ENDIF
284 
285  IF (gsc <= 0) THEN
286  CALL errorhint(65, 'subroutine SurfaceResistance.f95 (gsModel=1): gs <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
287  gsc = 0.1
288  ENDIF
289 
290  ELSEIF (gsmodel == 2 .OR. gsmodel == 4) THEN
291 
292  ! ---- g(kdown)----
293  qnm = kmax/(kmax + g2)
294  gq = (avkdn/(avkdn + g2))/qnm
295  IF (avkdn >= kmax) THEN !! Add proper error handling later - HCW!!
296  WRITE (*, *) 'Kmax exceeds Kdn setting to g(Kdn) to 1'
297  gq = 1
298  ENDIF
299 
300  ! ---- g(delq) ----
301  gdq = g3 + (1 - g3)*(g4**dq) !Ogink-Hendriks (1995) Eq 12 (using G3 as Kshd and G4 as r)
302 
303  ! ---- g(Tair) ----
304  tc = (th - g5)/(g5 - tl)
305  tc2 = (g5 - tl)*(th - g5)**tc
306  ! If air temperature below TL or above TH, then use value for TL+0.1 or TH-0.1
307  IF (temp_c <= tl) THEN
308  gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
309  ! Call error only if no snow on ground
310  IF (min(snowfrac(1), snowfrac(2), snowfrac(3), snowfrac(4), snowfrac(5), snowfrac(6)) /= 1) THEN
311  CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TL+0.1,Temp_C,id,it', &
312  REAL(Temp_c, KIND(1d0)), id_real, it)
313  ENDIF
314  ELSEIF (temp_c >= th) THEN
315  gtemp = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
316  CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TH-0.1,Temp_C,id,it', &
317  REAL(Temp_c, KIND(1d0)), id_real, it)
318  ELSE
319  gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
320  ENDIF
321  ! ---- g(smd) ----
322  sdp = s1/g6 + s2
323  IF (smdmethod > 0) THEN !Modified from ==1 to > 0 by HCW 31/07/2014
324  gs = (1 - exp(g6*(xsmd - sdp)))/(1 - exp(g6*(-sdp))) !Use measured smd
325  ELSE
326  !gs=1-EXP(g6*(vsmd-sdp)) !Use modelled smd
327  gs = (1 - exp(g6*(vsmd - sdp)))/(1 - exp(g6*(-sdp)))
328  IF (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) == 0 .OR. sfr(watersurf) == 1) THEN
329  gs = 0 !If no veg so no vsmd, or all water so no smd, set gs=0 HCW 21 Jul 2016
330  ENDIF
331  ENDIF
332 
333  !gs = gs*(1 - SUM(SnowFrac(1:6))/6)
334 
335  IF (gs < 0) THEN
336  CALL errorhint(65, &
337  'subroutine SurfaceResistance.f95 (gsModel=2): gs < 0 calculated, setting to 0.0001', &
338  gs, id_real, it)
339  gs = 0.0001
340  ENDIF
341 
342  ! ---- g(LAI) ----
343  gl = 0 !Initialise
344  ! DO iv=ivConif,ivGrass !For vegetated surfaces
345  DO iv = 1, 3 !For vegetated surfaces
346  ! gl=gl+(sfr(iv+2)*(1-SnowFrac(iv+2)))*LAI(id-1,iv)/LAIMax(iv)*MaxConductance(iv)
347  gl = gl + (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
348  ENDDO
349 
350  IF (avkdn <= 0) THEN !At nighttime set gsc at arbitrary low value: gsc=0.1 mm/s (Shuttleworth, 1988b)
351  gsc = 0.1
352  ELSE
353  ! Multiply parts together
354  gsc = (g1*gq*gdq*gtemp*gs*gl)
355  ENDIF
356 
357  IF (gsc <= 0) THEN
358  CALL errorhint(65, 'subroutine SurfaceResistance.f95 (gsModel=2): gsc <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
359  gsc = 0.1
360  ENDIF
361 
362  ELSEIF (gsmodel < 1 .OR. gsmodel > 4) THEN
363  CALL errorhint(71, 'Value of gsModel not recognised.', notused, notused, gsmodel)
364  ENDIF
365 
366  resistsurf = 1/(gsc/1000) ![s m-1]
367  gfunc = gdq*gtemp*gs*gq
368 
369  RETURN
370  END SUBROUTINE surfaceresistance
371 
372  SUBROUTINE boundarylayerresistance( &
373  zzd, & ! input: !Active measurement height (meas. height-displac. height)
374  z0m, & !Aerodynamic roughness length
375  avu1, & !Average wind speed
376  ustar, &! input/output:
377  rb)! output:
379  IMPLICIT NONE
380 
381  REAL(KIND(1d0)), INTENT(in)::zzd !Active measurement height (meas. height-displac. height)
382  REAL(KIND(1d0)), INTENT(in)::z0m !Aerodynamic roughness length
383  REAL(KIND(1d0)), INTENT(in)::avU1 !Average wind speed
384 
385  REAL(KIND(1d0)), INTENT(inout)::UStar!Friction velocity
386 
387  REAL(KIND(1d0)), INTENT(out)::rb !boundary layer resistance shuttleworth
388 
389  REAL(KIND(1d0)), PARAMETER :: k = 0.4
390 
391  IF (ustar < 0.01) THEN
392  ustar = avu1/log(zzd/z0m)*k
393  END IF
394 
395  rb = (1.1/ustar) + (5.6*(ustar**0.333333))!rb - boundary layer resistance shuttleworth
396 
397  RETURN
398  END SUBROUTINE boundarylayerresistance
399 
400  SUBROUTINE suews_cal_roughnessparameters( &
401  RoughLenMomMethod, &! input:
402  sfr, &! surface fractions
403  bldgH, EveTreeH, DecTreeH, &
404  porosity_id, FAIBldg, FAIEveTree, FAIDecTree, &
405  z0m_in, zdm_in, Z, &
406  planF, &! output:
407  Zh, z0m, zdm, ZZD)
408  ! Get surface covers and frontal area fractions (LJ 11/2010)
409  ! Last modified:
410  ! TS 18 Sep 2017 - added explicit interface
411  ! HCW 08 Feb 2017 - fixed bug in Zh between grids, added default z0m, zdm
412  ! HCW 03 Mar 2015
413  ! sg feb 2012 - made separate subroutine
414  !--------------------------------------------------------------------------------
415 
416  IMPLICIT NONE
417 
418  INTEGER, PARAMETER:: nsurf = 7 ! number of surface types
419  INTEGER, PARAMETER:: PavSurf = 1 !When all surfaces considered together (1-7)
420  INTEGER, PARAMETER:: BldgSurf = 2
421  INTEGER, PARAMETER:: ConifSurf = 3
422  INTEGER, PARAMETER:: DecidSurf = 4
423  INTEGER, PARAMETER:: GrassSurf = 5 !New surface classes: Grass = 5th/7 surfaces
424  INTEGER, PARAMETER:: BSoilSurf = 6 !New surface classes: Bare soil = 6th/7 surfaces
425  INTEGER, PARAMETER:: WaterSurf = 7
426 
427  INTEGER, INTENT(in) ::RoughLenMomMethod
428 
429  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::sfr! surface fractions
430 
431  REAL(KIND(1d0)), INTENT(in) ::bldgH
432  REAL(KIND(1d0)), INTENT(in) ::EveTreeH
433  REAL(KIND(1d0)), INTENT(in) ::DecTreeH
434  REAL(KIND(1d0)), INTENT(in) ::porosity_id
435  REAL(KIND(1d0)), INTENT(in) ::FAIBldg
436  REAL(KIND(1d0)), INTENT(in) ::FAIEveTree
437  REAL(KIND(1d0)), INTENT(in) ::FAIDecTree
438  REAL(KIND(1d0)), INTENT(in) ::z0m_in ! z0m set in SiteSelect
439  REAL(KIND(1d0)), INTENT(in) ::zdm_in ! zdm set in SiteSelect
440  REAL(KIND(1d0)), INTENT(in) ::Z
441 
442  REAL(KIND(1d0)), INTENT(out) ::planF
443  REAL(KIND(1d0)), INTENT(out) ::Zh ! effective height of bluff bodies
444  REAL(KIND(1d0)), INTENT(out) ::z0m ! aerodynamic roughness length
445  REAL(KIND(1d0)), INTENT(out) ::zdm ! zero-plance displacement
446  REAL(KIND(1d0)), INTENT(out) ::ZZD ! z-zdm
447 
448  REAL(KIND(1d0)) ::areaZh
449  INTEGER, PARAMETER :: notUsedI = -55
450  REAL(KIND(1d0)), PARAMETER:: notUsed = -55.5
451  REAL(KIND(1D0)):: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water !Default values for roughness lengths [m]
452 
453  areazh = (sfr(bldgsurf) + sfr(conifsurf) + sfr(decidsurf)) !Total area of buildings and trees
454 
455  ! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
456  z0m4paved = 0.003 !estimate
457  z0m4grass = 0.02
458  z0m4bsoil = 0.002
459  z0m4water = 0.0005
460 
461  !------------------------------------------------------------------------------
462  !If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
463  IF (areazh /= 0) THEN
464  zh = dot_product([bldgh, evetreeh, dectreeh*(1 - porosity_id)], sfr([bldgsurf, conifsurf, decidsurf]))/areazh
465  planf = dot_product([faibldg, faievetree, faidectree*(1 - porosity_id)], sfr([bldgsurf, conifsurf, decidsurf]))/areazh
466  ELSE
467  zh = 0 !Set Zh to zero if areaZh = 0
468  planf = 1e-5
469  ENDIF
470 
471  IF (zh /= 0) THEN
472  !Calculate z0m and zdm depending on the Z0 method
473  IF (roughlenmommethod == 2) THEN !Rule of thumb (G&O 1999)
474  z0m = 0.1*zh
475  zdm = 0.7*zh
476  ELSEIF (roughlenmommethod == 3) THEN !MacDonald 1998
477  IF (areazh /= 0) THEN !Plan area fraction
478  !planF=FAIBldg*sfr(BldgSurf)/areaZh+FAItree*sfr(ConifSurf)/areaZh+FAItree*(1-porosity_id)*sfr(DecidSurf)/areaZh
479  ! planF = FAIBldg*sfr(BldgSurf)/areaZh + FAIEveTree*sfr(ConifSurf)/areaZh + FAIDecTree*(1 - porosity_id)*sfr(DecidSurf)/areaZh
480  ELSE
481  ! planF = 0.00001
482  zh = 1
483  ENDIF
484  zdm = (1 + 4.43**(-sfr(bldgsurf))*(sfr(bldgsurf) - 1))*zh
485  z0m = ((1 - zdm/zh)*exp(-(0.5*1.0*1.2/0.4**2*(1 - zdm/zh)*planf)**(-0.5)))*zh
486  ENDIF
487  ELSEIF (zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
488  IF (areazh /= 0) CALL errorhint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, areazh, notusedi)
489  !Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
490  IF (areazh /= 1) THEN
491  z0m = (z0m4paved*sfr(pavsurf) &
492  + z0m4grass*sfr(grasssurf) &
493  + z0m4bsoil*sfr(bsoilsurf) &
494  + z0m4water*sfr(watersurf))/(1 - areazh)
495  zdm = 0
496  CALL errorhint(15, 'Setting z0m and zdm using default values', z0m, zdm, notusedi)
497  ELSEIF (areazh == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
498  z0m = 1
499  zdm = 7
500  CALL errorhint(15, 'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
501  ENDIF
502  ENDIF
503 
504  IF (roughlenmommethod == 1) THEN !use values set in SiteSelect
505  z0m = z0m_in
506  zdm = zdm_in
507  ENDIF
508 
509  zzd = z - zdm
510 
511  ! Error messages if aerodynamic parameters negative
512  IF (z0m < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, z0 < 0 m.', z0m, notused, notusedi)
513  IF (zdm < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, zd < 0 m.', zdm, notused, notusedi)
514  IF (zzd < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, (z-zd) < 0 m.', zzd, notused, notusedi)
515  END SUBROUTINE suews_cal_roughnessparameters
516 
517  FUNCTION cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar) RESULT(z0V)
518  ! TS 31 Jul 2018: make this a separate funciton for reuse
519  !Z0V roughness length for vapour
520  IMPLICIT NONE
521  INTEGER, INTENT(in)::RoughLenHeatMethod
522  REAL(KIND(1d0)), INTENT(in)::z0m!Aerodynamic roughness length
523  REAL(KIND(1d0)), INTENT(in)::VegFraction!Fraction of vegetation
524  REAL(KIND(1d0)), INTENT(in)::UStar!Friction velocity
525 
526  REAL(KIND(1d0))::z0V!roughness length for vapor/heat
527 
528  REAL(KIND(1d0)), PARAMETER:: muu = 1.46e-5 !molecular viscosity
529 
530  z0v = 0.01 ! initialise as 0.01
531 
532  !Z0V roughness length for vapour
533  IF (roughlenheatmethod == 1) THEN !Brutasert (1982) Z0v=z0/10(see Grimmond & Oke, 1986)
534  z0v = z0m/10
535  ELSEIF (roughlenheatmethod == 2) THEN ! Kawai et al. (2007)
536  !z0V=z0m*exp(2-(1.2-0.9*veg_fr**0.29)*(UStar*z0m/muu)**0.25)
537  ! Changed by HCW 05 Nov 2015 (veg_fr includes water; VegFraction = veg + bare soil)
538  z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
539  ELSEIF (roughlenheatmethod == 3) THEN
540  z0v = z0m*exp(-20.) ! Voogt and Grimmond, JAM, 2000
541  ELSEIF (roughlenheatmethod == 4) THEN
542  z0v = z0m*exp(2 - 1.29*(ustar*z0m/muu)**0.25) !See !Kanda and Moriwaki (2007),Loridan et al. (2010)
543  ELSEIF (roughlenheatmethod == 5) THEN
544  ! an adaptive way to determine z0v; TS 06 Feb 2020
545  if (vegfraction > .999) then
546  ! fully pervious surface
547  z0v = z0m/10
548  else
549  ! impervious-pervious mixed surface
550  z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
551  endif
552  ENDIF
553 
554  END FUNCTION cal_z0v
555 
556 end module resist_module
subroutine boundarylayerresistance(zzd, z0m, avU1, UStar, rb)
subroutine surfaceresistance(id, it, SMDMethod, SnowFrac, sfr, avkdn, Temp_C, dq, xsmd, vsmd, MaxConductance, LAIMax, LAI_id, gsModel, Kmax, G1, G2, G3, G4, G5, G6, TH, TL, S1, S2, gfunc, gsc, ResistSurf)
real(kind(1d0)) function cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar)
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) psih
subroutine suews_cal_roughnessparameters(RoughLenMomMethod, sfr, bldgH, EveTreeH, DecTreeH, porosity_id, FAIBldg, FAIEveTree, FAIDecTree, z0m_in, zdm_in, Z, planF, Zh, z0m, zdm, ZZD)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
subroutine aerodynamicresistance(ZZD, z0m, AVU1, L_mod, UStar, VegFraction, AerodynamicResistanceMethod, StabilityMethod, RoughLenHeatMethod, RA)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)