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
57 INTEGER,
PARAMETER :: notUsedI = -55
59 REAL(KIND(1d0)),
PARAMETER :: &
63 REAL(KIND(1d0)):: psim
68 IF (aerodynamicresistancemethod == 1)
THEN 69 ra = (log(zzd/z0m)**2)/(k2*avu1)
75 ELSEIF (aerodynamicresistancemethod == 2)
THEN 81 z0v =
cal_z0v(roughlenheatmethod, z0m, vegfraction, ustar)
83 IF (zzd/l_mod == 0 .OR. ustar == 0)
THEN 84 ra = (log(zzd/z0m)*log(zzd/z0v))/(k2*avu1)
86 ra = ((log(zzd/z0m) - psim)*(log(zzd/z0v) -
psih))/(k2*avu1)
90 ELSEIF (aerodynamicresistancemethod == 3)
THEN 91 ra = (4.72*log(zzd/z0m)**2)/(1 + 0.54*avu1)
96 CALL errorhint(7,
'In AerodynamicResistance.f95, calculated RA > 200 s m-1; RA set to 200 s m-1', ra, notused, notusedi)
99 CALL errorhint(7,
'In AerodynamicResistance.f95, calculated RA < 10 s m-1; RA set to 10 s m-1', ra, notused, notusedi)
102 IF (avu1 < 0)
WRITE (*, *) avu1, ra
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)
136 INTEGER,
PARAMETER::ConifSurf = 3
137 INTEGER,
PARAMETER::DecidSurf = 4
138 INTEGER,
PARAMETER::GrassSurf = 5
143 INTEGER,
PARAMETER::nsurf = 7
146 INTEGER,
PARAMETER::WaterSurf = 7
148 INTEGER,
INTENT(in)::id
149 INTEGER,
INTENT(in)::it
150 INTEGER,
INTENT(in)::gsModel
151 INTEGER,
INTENT(in)::SMDMethod
158 REAL(KIND(1d0)),
INTENT(in)::avkdn
159 REAL(KIND(1d0)),
INTENT(in)::Temp_C
160 REAL(KIND(1d0)),
INTENT(in)::Kmax
161 REAL(KIND(1d0)),
INTENT(in)::G1
162 REAL(KIND(1d0)),
INTENT(in)::G2
163 REAL(KIND(1d0)),
INTENT(in)::G3
164 REAL(KIND(1d0)),
INTENT(in)::G4
165 REAL(KIND(1d0)),
INTENT(in)::G5
166 REAL(KIND(1d0)),
INTENT(in)::G6
167 REAL(KIND(1d0)),
INTENT(in)::S1
168 REAL(KIND(1d0)),
INTENT(in)::S2
169 REAL(KIND(1d0)),
INTENT(in)::TH
170 REAL(KIND(1d0)),
INTENT(in)::TL
171 REAL(KIND(1d0)),
INTENT(in)::dq
172 REAL(KIND(1d0)),
INTENT(in)::xsmd
173 REAL(KIND(1d0)),
INTENT(in)::vsmd
175 REAL(KIND(1d0)),
DIMENSION(3),
INTENT(in) ::MaxConductance
176 REAL(KIND(1d0)),
DIMENSION(3),
INTENT(in) ::LAIMax
177 REAL(KIND(1d0)),
DIMENSION(3),
INTENT(in) ::LAI_id
178 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(in)::SnowFrac
179 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(in)::sfr
181 REAL(KIND(1d0)),
INTENT(out)::gfunc
182 REAL(KIND(1d0)),
INTENT(out)::gsc
183 REAL(KIND(1d0)),
INTENT(out)::ResistSurf
197 REAL(KIND(1d0)):: id_real
199 REAL(KIND(1d0)),
PARAMETER :: notUsed = -55
215 IF (gsmodel == 1 .OR. gsmodel == 3)
THEN 217 qnm = kmax/(kmax + g2)
219 gq = (avkdn/(g2 + avkdn))/qnm
227 tc = (th - g5)/(g5 - tl)
228 tc2 = (g5 - tl)*(th - g5)**tc
230 IF (temp_c <= tl)
THEN 231 gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
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)
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)
243 gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
248 IF (smdmethod > 0)
THEN 249 gs = 1 - exp(g6*(xsmd - sdp))
251 gs = 1 - exp(g6*(vsmd - sdp))
252 IF (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) == 0 .OR. sfr(watersurf) == 1)
THEN 259 'subroutine SurfaceResistance.f95 (gsModel=1): g(smd) < 0 calculated, setting to 0.0001', &
275 gl = gl + (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
282 gsc = (g1*gq*gdq*gtemp*gs*gl)
286 CALL errorhint(65,
'subroutine SurfaceResistance.f95 (gsModel=1): gs <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
290 ELSEIF (gsmodel == 2 .OR. gsmodel == 4)
THEN 293 qnm = kmax/(kmax + g2)
294 gq = (avkdn/(avkdn + g2))/qnm
295 IF (avkdn >= kmax)
THEN 296 WRITE (*, *)
'Kmax exceeds Kdn setting to g(Kdn) to 1' 301 gdq = g3 + (1 - g3)*(g4**dq)
304 tc = (th - g5)/(g5 - tl)
305 tc2 = (g5 - tl)*(th - g5)**tc
307 IF (temp_c <= tl)
THEN 308 gtemp = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
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)
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)
319 gtemp = (temp_c - tl)*(th - temp_c)**tc/tc2
323 IF (smdmethod > 0)
THEN 324 gs = (1 - exp(g6*(xsmd - sdp)))/(1 - exp(g6*(-sdp)))
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 337 'subroutine SurfaceResistance.f95 (gsModel=2): gs < 0 calculated, setting to 0.0001', &
347 gl = gl + (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
354 gsc = (g1*gq*gdq*gtemp*gs*gl)
358 CALL errorhint(65,
'subroutine SurfaceResistance.f95 (gsModel=2): gsc <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
362 ELSEIF (gsmodel < 1 .OR. gsmodel > 4)
THEN 363 CALL errorhint(71,
'Value of gsModel not recognised.', notused, notused, gsmodel)
366 resistsurf = 1/(gsc/1000)
367 gfunc = gdq*gtemp*gs*gq
373 zzd, & ! input: !Active measurement height (meas. height-displac. height)
381 REAL(KIND(1d0)),
INTENT(in)::zzd
382 REAL(KIND(1d0)),
INTENT(in)::z0m
383 REAL(KIND(1d0)),
INTENT(in)::avU1
385 REAL(KIND(1d0)),
INTENT(inout)::UStar
387 REAL(KIND(1d0)),
INTENT(out)::rb
389 REAL(KIND(1d0)),
PARAMETER :: k = 0.4
391 IF (ustar < 0.01)
THEN 392 ustar = avu1/log(zzd/z0m)*k
395 rb = (1.1/ustar) + (5.6*(ustar**0.333333))
401 RoughLenMomMethod, &! input:
402 sfr, &! surface fractions
403 bldgH, EveTreeH, DecTreeH, &
404 porosity_id, FAIBldg, FAIEveTree, FAIDecTree, &
418 INTEGER,
PARAMETER:: nsurf = 7
419 INTEGER,
PARAMETER:: PavSurf = 1
420 INTEGER,
PARAMETER:: BldgSurf = 2
421 INTEGER,
PARAMETER:: ConifSurf = 3
422 INTEGER,
PARAMETER:: DecidSurf = 4
423 INTEGER,
PARAMETER:: GrassSurf = 5
424 INTEGER,
PARAMETER:: BSoilSurf = 6
425 INTEGER,
PARAMETER:: WaterSurf = 7
427 INTEGER,
INTENT(in) ::RoughLenMomMethod
429 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(in) ::sfr
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
439 REAL(KIND(1d0)),
INTENT(in) ::zdm_in
440 REAL(KIND(1d0)),
INTENT(in) ::Z
442 REAL(KIND(1d0)),
INTENT(out) ::planF
443 REAL(KIND(1d0)),
INTENT(out) ::Zh
444 REAL(KIND(1d0)),
INTENT(out) ::z0m
445 REAL(KIND(1d0)),
INTENT(out) ::zdm
446 REAL(KIND(1d0)),
INTENT(out) ::ZZD
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
453 areazh = (sfr(bldgsurf) + sfr(conifsurf) + sfr(decidsurf))
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
473 IF (roughlenmommethod == 2)
THEN 476 ELSEIF (roughlenmommethod == 3)
THEN 477 IF (areazh /= 0)
THEN 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
487 ELSEIF (zh == 0)
THEN 488 IF (areazh /= 0)
CALL errorhint(15,
'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, areazh, notusedi)
490 IF (areazh /= 1)
THEN 491 z0m = (z0m4paved*sfr(pavsurf) &
492 + z0m4grass*sfr(grasssurf) &
493 + z0m4bsoil*sfr(bsoilsurf) &
494 + z0m4water*sfr(watersurf))/(1 - areazh)
496 CALL errorhint(15,
'Setting z0m and zdm using default values', z0m, zdm, notusedi)
497 ELSEIF (areazh == 1)
THEN 500 CALL errorhint(15,
'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
504 IF (roughlenmommethod == 1)
THEN 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)
517 FUNCTION cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar)
RESULT(z0V)
521 INTEGER,
INTENT(in)::RoughLenHeatMethod
522 REAL(KIND(1d0)),
INTENT(in)::z0m
523 REAL(KIND(1d0)),
INTENT(in)::VegFraction
524 REAL(KIND(1d0)),
INTENT(in)::UStar
528 REAL(KIND(1d0)),
PARAMETER:: muu = 1.46e-5
533 IF (roughlenheatmethod == 1)
THEN 535 ELSEIF (roughlenheatmethod == 2)
THEN 538 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
539 ELSEIF (roughlenheatmethod == 3)
THEN 541 ELSEIF (roughlenheatmethod == 4)
THEN 542 z0v = z0m*exp(2 - 1.29*(ustar*z0m/muu)**0.25)
543 ELSEIF (roughlenheatmethod == 5)
THEN 545 if (vegfraction > .999)
then 550 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
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)
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)