12 AerodynamicResistanceMethod, &
44 REAL(KIND(1D0)),
INTENT(in) :: ZZD
45 REAL(KIND(1D0)),
INTENT(in) :: z0m
46 REAL(KIND(1D0)),
INTENT(in) :: AVU1
47 REAL(KIND(1D0)),
INTENT(in) :: L_mod
48 REAL(KIND(1D0)),
INTENT(in) :: UStar
49 REAL(KIND(1D0)),
INTENT(in) :: VegFraction
51 INTEGER,
INTENT(in) :: AerodynamicResistanceMethod
52 INTEGER,
INTENT(in) :: StabilityMethod
53 INTEGER,
INTENT(in) :: RoughLenHeatMethod
55 REAL(KIND(1D0)),
INTENT(out) :: RA_h
56 REAL(KIND(1D0)),
INTENT(out) :: z0V
58 INTEGER,
PARAMETER :: notUsedI = -55
60 REAL(KIND(1D0)),
PARAMETER :: &
64 REAL(KIND(1D0)) :: psim
68 z0v =
cal_z0v(roughlenheatmethod, z0m, vegfraction, ustar)
71 IF (aerodynamicresistancemethod == 1)
THEN
72 ra_h = (log(zzd/z0m)**2)/(k2*avu1)
78 ELSEIF (aerodynamicresistancemethod == 2)
THEN
85 IF (zzd/l_mod == 0 .OR. ustar == 0)
THEN
86 ra_h = (log(zzd/z0m)*log(zzd/z0v))/(k2*avu1)
88 ra_h = ((log(zzd/z0m) - psim)*(log(zzd/z0v) -
psih))/(k2*avu1)
93 ELSEIF (aerodynamicresistancemethod == 3)
THEN
94 ra_h = (4.72*log(zzd/z0m)**2)/(1 + 0.54*avu1)
99 CALL errorhint(7,
'In AerodynamicResistance.f95, calculated RA > 200 s m-1; RA set to 200 s m-1', ra_h, notused, notusedi)
101 ELSEIF (ra_h < 10)
THEN
102 CALL errorhint(7,
'In AerodynamicResistance.f95, calculated RA < 10 s m-1; RA set to 10 s m-1', ra_h, notused, notusedi)
105 IF (avu1 < 0)
WRITE (*, *) avu1, ra_h
113 SMDMethod, SnowFrac, sfr_surf, avkdn, Temp_C, dq, xsmd, vsmd, MaxConductance, &
114 LAIMax, LAI_id, gsModel, Kmax, &
115 G1, G2, G3, G4, G5, G6, TH, TL, S1, S2, &
116 gfunc, gsc, ResistSurf)
139 INTEGER,
PARAMETER :: ConifSurf = 3
140 INTEGER,
PARAMETER :: DecidSurf = 4
141 INTEGER,
PARAMETER :: GrassSurf = 5
146 INTEGER,
PARAMETER :: nsurf = 7
149 INTEGER,
PARAMETER :: WaterSurf = 7
151 INTEGER,
INTENT(in) :: id
152 INTEGER,
INTENT(in) :: it
153 INTEGER,
INTENT(in) :: gsModel
154 INTEGER,
INTENT(in) :: SMDMethod
161 REAL(KIND(1D0)),
INTENT(in) :: avkdn
162 REAL(KIND(1D0)),
INTENT(in) :: Temp_C
163 REAL(KIND(1D0)),
INTENT(in) :: Kmax
164 REAL(KIND(1D0)),
INTENT(in) :: G1
165 REAL(KIND(1D0)),
INTENT(in) :: G2
166 REAL(KIND(1D0)),
INTENT(in) :: G3
167 REAL(KIND(1D0)),
INTENT(in) :: G4
168 REAL(KIND(1D0)),
INTENT(in) :: G5
169 REAL(KIND(1D0)),
INTENT(in) :: G6
170 REAL(KIND(1D0)),
INTENT(in) :: S1
171 REAL(KIND(1D0)),
INTENT(in) :: S2
172 REAL(KIND(1D0)),
INTENT(in) :: TH
173 REAL(KIND(1D0)),
INTENT(in) :: TL
174 REAL(KIND(1D0)),
INTENT(in) :: dq
175 REAL(KIND(1D0)),
INTENT(in) :: xsmd
176 REAL(KIND(1D0)),
INTENT(in) :: vsmd
178 REAL(KIND(1D0)),
DIMENSION(3),
INTENT(in) :: MaxConductance
179 REAL(KIND(1D0)),
DIMENSION(3),
INTENT(in) :: LAIMax
180 REAL(KIND(1D0)),
DIMENSION(3),
INTENT(in) :: LAI_id
181 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: SnowFrac
182 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: sfr_surf
184 REAL(KIND(1D0)),
INTENT(out) :: gfunc
185 REAL(KIND(1D0)),
INTENT(out) :: gsc
186 REAL(KIND(1D0)),
INTENT(out) :: ResistSurf
200 REAL(KIND(1D0)) :: id_real
202 REAL(KIND(1D0)),
PARAMETER :: notUsed = -55
218 IF (gsmodel == 1 .OR. gsmodel == 3)
THEN
220 qnm = kmax/(kmax + g2)
222 gq = (avkdn/(g2 + avkdn))/qnm
230 tc = (th - g5)/(g5 - tl)
231 tc2 = (g5 - tl)*(th - g5)**tc
233 IF (temp_c <= tl)
THEN
234 gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
237 IF (minval(snowfrac(1:6)) /= 1)
THEN
238 CALL errorhint(29,
'subroutine SurfaceResistance.f95: T changed to fit limits TL=0.1,Temp_c,id,it', &
239 REAL(Temp_c, KIND(1D0)), id_real, it)
241 ELSEIF (temp_c >= th)
THEN
242 gtemp = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
243 CALL errorhint(29,
'subroutine SurfaceResistance.f95: T changed to fit limits TH=39.9,Temp_c,id,it', &
244 REAL(Temp_c, KIND(1D0)), id_real, it)
246 gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
251 IF (smdmethod > 0)
THEN
252 gs = 1 - exp(g6*(xsmd - sdp))
254 gs = 1 - exp(g6*(vsmd - sdp))
255 IF (sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf) == 0 .OR. sfr_surf(watersurf) == 1)
THEN
262 'subroutine SurfaceResistance.f95 (gsModel=1): g(smd) < 0 calculated, setting to 0.0001', &
278 gl = gl + (sfr_surf(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
285 gsc = (g1*gq*gdq*gtemp*gs*gl)
289 CALL errorhint(65,
'subroutine SurfaceResistance.f95 (gsModel=1): gs <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
293 ELSEIF (gsmodel == 2 .OR. gsmodel == 4)
THEN
296 qnm = kmax/(kmax + g2)
297 gq = (avkdn/(avkdn + g2))/qnm
298 IF (avkdn >= kmax)
THEN
299 WRITE (*, *)
'Kmax exceeds Kdn setting to g(Kdn) to 1'
304 gdq = g3 + (1 - g3)*(g4**dq)
307 tc = (th - g5)/(g5 - tl)
308 tc2 = (g5 - tl)*(th - g5)**tc
310 IF (temp_c <= tl)
THEN
311 gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
313 IF (min(snowfrac(1), snowfrac(2), snowfrac(3), snowfrac(4), snowfrac(5), snowfrac(6)) /= 1)
THEN
314 CALL errorhint(29,
'subroutine SurfaceResistance.f95: T changed to fit limits TL+0.1,Temp_C,id,it', &
315 REAL(Temp_c, KIND(1D0)), id_real, it)
317 ELSEIF (temp_c >= th)
THEN
318 gtemp = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
319 CALL errorhint(29,
'subroutine SurfaceResistance.f95: T changed to fit limits TH-0.1,Temp_C,id,it', &
320 REAL(Temp_c, KIND(1D0)), id_real, it)
322 gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
326 IF (smdmethod > 0)
THEN
327 gs = (1 - exp(g6*(xsmd - sdp)))/(1 - exp(g6*(-sdp)))
330 gs = (1 - exp(g6*(vsmd - sdp)))/(1 - exp(g6*(-sdp)))
331 IF (sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf) == 0 .OR. sfr_surf(watersurf) == 1)
THEN
340 'subroutine SurfaceResistance.f95 (gsModel=2): gs < 0 calculated, setting to 0.0001', &
350 gl = gl + (sfr_surf(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
357 gsc = (g1*gq*gdq*gtemp*gs*gl)
361 CALL errorhint(65,
'subroutine SurfaceResistance.f95 (gsModel=2): gsc <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
365 ELSEIF (gsmodel < 1 .OR. gsmodel > 4)
THEN
366 CALL errorhint(71,
'Value of gsModel not recognised.', notused, notused, gsmodel)
369 resistsurf = 1/(gsc/1000)
370 gfunc = gdq*gtemp*gs*gq
376 zzd, & ! input: !Active measurement height (meas. height-displac. height)
384 REAL(KIND(1D0)),
INTENT(in) :: zzd
385 REAL(KIND(1D0)),
INTENT(in) :: z0m
386 REAL(KIND(1D0)),
INTENT(in) :: avU1
388 REAL(KIND(1D0)),
INTENT(inout) :: UStar
390 REAL(KIND(1D0)),
INTENT(out) :: rb
392 REAL(KIND(1D0)),
PARAMETER :: k = 0.4
394 IF (ustar < 0.01)
THEN
395 ustar = avu1/log(zzd/z0m)*k
398 rb = (1.1/ustar) + (5.6*(ustar**0.333333))
404 RoughLenMomMethod, & ! input:
405 sfr_surf, & ! surface fractions
406 bldgH, EveTreeH, DecTreeH, &
407 porosity_dectr, FAIBldg, FAIEveTree, FAIDecTree, &
409 FAI, PAI, & ! output:
421 INTEGER,
PARAMETER :: nsurf = 7
422 INTEGER,
PARAMETER :: PavSurf = 1
423 INTEGER,
PARAMETER :: BldgSurf = 2
424 INTEGER,
PARAMETER :: ConifSurf = 3
425 INTEGER,
PARAMETER :: DecidSurf = 4
426 INTEGER,
PARAMETER :: GrassSurf = 5
427 INTEGER,
PARAMETER :: BSoilSurf = 6
428 INTEGER,
PARAMETER :: WaterSurf = 7
429 REAL(KIND(1D0)),
PARAMETER :: porosity_evetr = 0.32
431 INTEGER,
INTENT(in) :: RoughLenMomMethod
433 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: sfr_surf
435 REAL(KIND(1D0)),
INTENT(in) :: bldgH
436 REAL(KIND(1D0)),
INTENT(in) :: EveTreeH
437 REAL(KIND(1D0)),
INTENT(in) :: DecTreeH
438 REAL(KIND(1D0)),
INTENT(in) :: porosity_dectr
439 REAL(KIND(1D0)),
INTENT(in) :: FAIBldg
440 REAL(KIND(1D0)),
INTENT(in) :: FAIEveTree
441 REAL(KIND(1D0)),
INTENT(in) :: FAIDecTree
442 REAL(KIND(1D0)),
INTENT(in) :: z0m_in
443 REAL(KIND(1D0)),
INTENT(in) :: zdm_in
444 REAL(KIND(1D0)),
INTENT(in) :: Z
446 REAL(KIND(1D0)),
INTENT(out) :: FAI
447 REAL(KIND(1D0)),
INTENT(out) :: PAI
448 REAL(KIND(1D0)),
INTENT(out) :: Zh
449 REAL(KIND(1D0)),
INTENT(out) :: z0m
450 REAL(KIND(1D0)),
INTENT(out) :: zdm
451 REAL(KIND(1D0)),
INTENT(out) :: ZZD
453 INTEGER,
PARAMETER :: notUsedI = -55
454 REAL(KIND(1D0)),
PARAMETER :: notUsed = -55.5
455 REAL(KIND(1D0)) :: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water
460 pai = dot_product(sfr_surf([bldgsurf, conifsurf, decidsurf]), [1d0, 1 - porosity_evetr, 1 - porosity_dectr])
472 [bldgh, evetreeh*(1 - porosity_evetr), dectreeh*(1 - porosity_dectr)], &
473 sfr_surf([bldgsurf, conifsurf, decidsurf]))/pai
474 fai = sum(merge([faibldg, faievetree*(1 - porosity_evetr), faidectree*(1 - porosity_dectr)], &
476 sfr_surf([bldgsurf, conifsurf, decidsurf]) > 0))
488 IF (roughlenmommethod == 2)
THEN
491 ELSEIF (roughlenmommethod == 3)
THEN
492 zdm = (1 + 4.43**(-sfr_surf(bldgsurf))*(sfr_surf(bldgsurf) - 1))*zh
493 z0m = ((1 - zdm/zh)*exp(-(0.5*1.0*1.2/0.4**2*(1 - zdm/zh)*fai)**(-0.5)))*zh
494 ELSEIF (roughlenmommethod == 4)
THEN
496 zdm = (-0.182 + 0.722*
sigmoid(-1.16 + 3.89*pai) + 0.493*
sigmoid(-5.17 + 32.7*pai))*zh
498 0.0165*min(pai, .7) + 2.52*min(pai, .7)**2 + &
499 3.21*min(pai, .7)**3 - 43.6*min(pai, .7)**4 + &
500 76.5*min(pai, .7)**5 - 40.*min(pai, .7)**6)*zh
502 ELSEIF (zh == 0)
THEN
503 IF (pai /= 0)
CALL errorhint(15,
'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, pai, notusedi)
506 z0m = (z0m4paved*sfr_surf(pavsurf) &
507 + z0m4grass*sfr_surf(grasssurf) &
508 + z0m4bsoil*sfr_surf(bsoilsurf) &
509 + z0m4water*sfr_surf(watersurf))/(1 - pai)
511 CALL errorhint(15,
'Setting z0m and zdm using default values', z0m, zdm, notusedi)
512 ELSEIF (pai == 1)
THEN
515 CALL errorhint(15,
'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
519 IF (roughlenmommethod == 1)
THEN
527 IF (z0m < 0)
CALL errorhint(14,
'In SUEWS_cal_RoughnessParameters, z0 < 0 m.', z0m, notused, notusedi)
528 IF (zdm < 0)
CALL errorhint(14,
'In SUEWS_cal_RoughnessParameters, zd < 0 m.', zdm, notused, notusedi)
529 IF (zzd < 0)
CALL errorhint(14,
'In SUEWS_cal_RoughnessParameters, (z-zd) < 0 m.', zzd, notused, notusedi)
532 FUNCTION cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar)
RESULT(z0V)
536 INTEGER,
INTENT(in) :: roughlenheatmethod
537 REAL(kind(1d0)),
INTENT(in) :: z0m
538 REAL(kind(1d0)),
INTENT(in) :: vegfraction
539 REAL(kind(1d0)),
INTENT(in) :: ustar
541 REAL(kind(1d0)) :: z0v
543 REAL(kind(1d0)),
PARAMETER :: muu = 1.46e-5
548 IF (roughlenheatmethod == 1)
THEN
550 ELSEIF (roughlenheatmethod == 2)
THEN
553 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
554 ELSEIF (roughlenheatmethod == 3)
THEN
556 ELSEIF (roughlenheatmethod == 4)
THEN
557 z0v = z0m*exp(2 - 1.29*(ustar*z0m/muu)**0.25)
558 ELSEIF (roughlenheatmethod == 5)
THEN
560 IF (vegfraction > .999)
THEN
565 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
573 REAL(kind(1d0)),
INTENT(in) :: x
574 REAL(kind(1d0)) :: res
576 res = 1/(1 + exp(-x))
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
subroutine surfaceresistance(id, it, SMDMethod, SnowFrac, sfr_surf, 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)
subroutine suews_cal_roughnessparameters(RoughLenMomMethod, sfr_surf, bldgH, EveTreeH, DecTreeH, porosity_dectr, FAIBldg, FAIEveTree, FAIDecTree, z0m_in, zdm_in, Z, FAI, PAI, Zh, z0m, zdm, ZZD)
subroutine aerodynamicresistance(ZZD, z0m, AVU1, L_mod, UStar, VegFraction, AerodynamicResistanceMethod, StabilityMethod, RoughLenHeatMethod, RA_h, z0V)
subroutine boundarylayerresistance(zzd, z0m, avU1, UStar, rb)
real(kind(1d0)) function sigmoid(x)
real(kind(1d0)) function cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)