SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
resist_module Module Reference

Functions/Subroutines

subroutine aerodynamicresistance (zzd, z0m, avu1, l_mod, ustar, vegfraction, aerodynamicresistancemethod, stabilitymethod, roughlenheatmethod, ra_h, z0v)
 
subroutine surfaceresistance (id, it, smdmethod, snowfrac, sfr_surf, avkdn, temp_c, dq, xsmd, vsmd, maxconductance, laimax, lai_id, gsmodel, kmax, g_max, g_k, g_q_base, g_q_shape, g_t, g_sm, th, tl, s1, s2, g_kdown, g_dq, g_ta, g_smd, g_lai, gfunc, gsc, rs)
 
subroutine boundarylayerresistance (zzd, z0m, avu1, ustar, rb)
 
subroutine suews_cal_roughnessparameters (roughlenmommethod, faimethod, sfr_surf, surfacearea, bldgh, evetreeh, dectreeh, porosity_dectr, faibldg, faievetree, faidectree, z0m_in, zdm_in, z, fai, pai, zh, z0m, zdm, zzd)
 
subroutine suews_cal_roughnessparameters_dts (roughlenmommethod, faimethod, sfr_paved, sfr_bldg, sfr_evetr, sfr_dectr, sfr_grass, sfr_bsoil, sfr_water, surfacearea, bldgh, evetreeh, dectreeh, porosity_dectr, faibldg, faievetree, faidectree, z0m_in, zdm_in, z, faibldg_use, faievetree_use, faidectree_use, fai, pai, zh, z0m, zdm, zzd)
 
real(kind(1d0)) function cal_z0v (roughlenheatmethod, z0m, vegfraction, ustar)
 
real(kind(1d0)) function sigmoid (x)
 

Function/Subroutine Documentation

◆ aerodynamicresistance()

subroutine resist_module::aerodynamicresistance ( real(kind(1d0)), intent(in) zzd,
real(kind(1d0)), intent(in) z0m,
real(kind(1d0)), intent(in) avu1,
real(kind(1d0)), intent(in) l_mod,
real(kind(1d0)), intent(in) ustar,
real(kind(1d0)), intent(in) vegfraction,
integer, intent(in) aerodynamicresistancemethod,
integer, intent(in) stabilitymethod,
integer, intent(in) roughlenheatmethod,
real(kind(1d0)), intent(out) ra_h,
real(kind(1d0)), intent(out) z0v )

Definition at line 5 of file suews_phys_resist.f95.

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
39 USE atmmoiststab_module, ONLY: stab_psi_heat, stab_psi_mom
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_h !Aerodynamic resistance for heat/vapour [s m^-1]
56 REAL(KIND(1D0)), INTENT(out) :: z0V
57
58 INTEGER, PARAMETER :: notUsedI = -55
59
60 REAL(KIND(1D0)), PARAMETER :: &
61 notused = -55.5, &
62 k2 = 0.16, & !Power of Van Karman's constant (= 0.16 = 0.4^2)
63 muu = 1.46e-5 !molecular viscosity
64 REAL(KIND(1D0)) :: psim
65 ! REAL(KIND(1d0)):: psih
66
67 !Z0V roughness length for vapour
68 z0v = cal_z0v(roughlenheatmethod, z0m, vegfraction, ustar)
69
70 !1)Monteith (1965)-neutral stability
71 IF (aerodynamicresistancemethod == 1) THEN
72 ra_h = (log(zzd/z0m)**2)/(k2*avu1)
73
74 !2) Non-neutral stability
75 ! PSIM - stability function for momentum
76 ! psih - stability function for heat
77 ! assuming stability functions the same for heat and water
78 ELSEIF (aerodynamicresistancemethod == 2) THEN !Dyer (1974)
79
80 ! psim = stab_psi_mom(StabilityMethod, zzd/L_mod)
81 ! psih = stab_psi_heat(StabilityMethod, ZZD/L_mod)
82 psim = stab_psi_mom(stabilitymethod, zzd/l_mod) - stab_psi_mom(stabilitymethod, z0m/l_mod)
83 psih = stab_psi_heat(stabilitymethod, zzd/l_mod) - stab_psi_heat(stabilitymethod, z0v/l_mod)
84
85 IF (zzd/l_mod == 0 .OR. ustar == 0) THEN
86 ra_h = (log(zzd/z0m)*log(zzd/z0v))/(k2*avu1) !Use neutral equation
87 ELSE
88 ra_h = ((log(zzd/z0m) - psim)*(log(zzd/z0v) - psih))/(k2*avu1)
89 ! RA = AVU1/UStar**2
90 END IF
91
92 !3) Thom and Oliver (1977)
93 ELSEIF (aerodynamicresistancemethod == 3) THEN
94 ra_h = (4.72*log(zzd/z0m)**2)/(1 + 0.54*avu1)
95 END IF
96
97 !If RA outside permitted range, adjust extreme values !!Check whether these thresholds are suitable over a range of z0
98 IF (ra_h > 120) THEN !was 175
99 CALL errorhint(7, 'In AerodynamicResistance.f95, calculated RA > 200 s m-1; RA set to 200 s m-1', ra_h, notused, notusedi)
100 ra_h = 120
101 ELSEIF (ra_h < 10) THEN !found By Shiho - fix Dec 2012 !Threshold changed from 2 to 10 s m-1 (HCW 03 Dec 2015)
102 CALL errorhint(7, 'In AerodynamicResistance.f95, calculated RA < 10 s m-1; RA set to 10 s m-1', ra_h, notused, notusedi)
103 ra_h = 10
104 ! RA=(log(ZZD/z0m))**2/(k2*AVU1)
105 IF (avu1 < 0) WRITE (*, *) avu1, ra_h
106 END IF
107
108 RETURN
atmmoiststab_module
Definition suews_phys_atmmoiststab.f95:1
atmmoiststab_module::stab_psi_mom
real(kind(1d0)) function stab_psi_mom(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:272
atmmoiststab_module::stab_psi_heat
real(kind(1d0)) function stab_psi_heat(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:293
sues_data
Definition suews_ctrl_const.f95:1526
sues_data::stabilitymethod
integer stabilitymethod
Definition suews_ctrl_const.f95:1544
sues_data::aerodynamicresistancemethod
integer aerodynamicresistancemethod
Definition suews_ctrl_const.f95:1551
sues_data::ustar
real(kind(1d0)) ustar
Definition suews_ctrl_const.f95:1625
sues_data::l_mod
real(kind(1d0)) l_mod
Definition suews_ctrl_const.f95:1625
sues_data::psih
real(kind(1d0)) psih
Definition suews_ctrl_const.f95:1625
sues_data::roughlenheatmethod
integer roughlenheatmethod
Definition suews_ctrl_const.f95:1544
sues_data::psim
real(kind(1d0)) psim
Definition suews_ctrl_const.f95:1625
errorhint
subroutine errorhint(errh, problemfile, value, value2, valuei)
Definition suews_ctrl_error.f95:2

References cal_z0v(), errorhint(), sues_data::psih, atmmoiststab_module::stab_psi_heat(), and atmmoiststab_module::stab_psi_mom().

Referenced by suews_driver::suews_cal_resistance(), and suews_driver::suews_cal_resistance_dts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ boundarylayerresistance()

subroutine resist_module::boundarylayerresistance ( real(kind(1d0)), intent(in) zzd,
real(kind(1d0)), intent(in) z0m,
real(kind(1d0)), intent(in) avu1,
real(kind(1d0)), intent(inout) ustar,
real(kind(1d0)), intent(out) rb )

Definition at line 381 of file suews_phys_resist.f95.

387
388 IMPLICIT NONE
389
390 REAL(KIND(1D0)), INTENT(in) :: zzd !Active measurement height (meas. height-displac. height)
391 REAL(KIND(1D0)), INTENT(in) :: z0m !Aerodynamic roughness length
392 REAL(KIND(1D0)), INTENT(in) :: avU1 !Average wind speed
393
394 REAL(KIND(1D0)), INTENT(inout) :: UStar !Friction velocity
395
396 REAL(KIND(1D0)), INTENT(out) :: rb !boundary layer resistance shuttleworth
397
398 REAL(KIND(1D0)), PARAMETER :: k = 0.4
399
400 IF (ustar < 0.01) THEN
401 ustar = avu1/log(zzd/z0m)*k
402 END IF
403
404 rb = (1.1/ustar) + (5.6*(ustar**0.333333)) !rb - boundary layer resistance shuttleworth
405
406 RETURN

Referenced by suews_driver::suews_cal_resistance(), and suews_driver::suews_cal_resistance_dts().

Here is the caller graph for this function:

◆ cal_z0v()

real(kind(1d0)) function resist_module::cal_z0v ( integer, intent(in) roughlenheatmethod,
real(kind(1d0)), intent(in) z0m,
real(kind(1d0)), intent(in) vegfraction,
real(kind(1d0)), intent(in) ustar )

Definition at line 731 of file suews_phys_resist.f95.

732 ! TS 31 Jul 2018: make this a separate funciton for reuse
733 !Z0V roughness length for vapour
734 IMPLICIT NONE
735 INTEGER, INTENT(in) :: RoughLenHeatMethod
736 REAL(KIND(1D0)), INTENT(in) :: z0m !Aerodynamic roughness length
737 REAL(KIND(1D0)), INTENT(in) :: VegFraction !Fraction of vegetation
738 REAL(KIND(1D0)), INTENT(in) :: UStar !Friction velocity
739
740 REAL(KIND(1D0)) :: z0V !roughness length for vapor/heat
741
742 REAL(KIND(1D0)), PARAMETER :: muu = 1.46e-5 !molecular viscosity
743
744 z0v = 0.01 ! initialise as 0.01
745
746 !Z0V roughness length for vapour
747 IF (roughlenheatmethod == 1) THEN !Brutsaert (1982) Z0v=z0/10(see Grimmond & Oke, 1986)
748 z0v = z0m/10
749 ELSEIF (roughlenheatmethod == 2) THEN ! Kawai et al. (2007)
750 !z0V=z0m*exp(2-(1.2-0.9*veg_fr**0.29)*(UStar*z0m/muu)**0.25)
751 ! Changed by HCW 05 Nov 2015 (veg_fr includes water; VegFraction = veg + bare soil)
752 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
753 ELSEIF (roughlenheatmethod == 3) THEN
754 z0v = z0m*exp(-20.) ! Voogt and Grimmond, JAM, 2000
755 ELSEIF (roughlenheatmethod == 4) THEN
756 z0v = z0m*exp(2 - 1.29*(ustar*z0m/muu)**0.25) !See !Kanda and Moriwaki (2007),Loridan et al. (2010)
757 ELSEIF (roughlenheatmethod == 5) THEN
758 ! an adaptive way to determine z0v; TS 06 Feb 2020
759 IF (vegfraction > .999) THEN
760 ! fully pervious surface
761 z0v = z0m/10
762 ELSE
763 ! impervious-pervious mixed surface
764 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
765 END IF
766 END IF
767

Referenced by aerodynamicresistance().

Here is the caller graph for this function:

◆ sigmoid()

real(kind(1d0)) function resist_module::sigmoid ( real(kind(1d0)), intent(in) x)

Definition at line 770 of file suews_phys_resist.f95.

771 IMPLICIT NONE
772 REAL(KIND(1D0)), INTENT(in) :: x
773 REAL(KIND(1D0)) :: res
774
775 res = 1/(1 + exp(-x))
776

Referenced by suews_cal_roughnessparameters(), and suews_cal_roughnessparameters_dts().

Here is the caller graph for this function:

◆ suews_cal_roughnessparameters()

subroutine resist_module::suews_cal_roughnessparameters ( integer, intent(in) roughlenmommethod,
integer, intent(in) faimethod,
real(kind(1d0)), dimension(nsurf), intent(in) sfr_surf,
real(kind(1d0)), intent(in) surfacearea,
real(kind(1d0)), intent(in) bldgh,
real(kind(1d0)), intent(in) evetreeh,
real(kind(1d0)), intent(in) dectreeh,
real(kind(1d0)), intent(in) porosity_dectr,
real(kind(1d0)), intent(in) faibldg,
real(kind(1d0)), intent(in) faievetree,
real(kind(1d0)), intent(in) faidectree,
real(kind(1d0)), intent(in) z0m_in,
real(kind(1d0)), intent(in) zdm_in,
real(kind(1d0)), intent(in) z,
real(kind(1d0)), intent(out) fai,
real(kind(1d0)), intent(out) pai,
real(kind(1d0)), intent(out) zh,
real(kind(1d0)), intent(out) z0m,
real(kind(1d0)), intent(out) zdm,
real(kind(1d0)), intent(out) zzd )

Definition at line 409 of file suews_phys_resist.f95.

418 ! Get surface covers and frontal area fractions (LJ 11/2010)
419 ! Last modified:
420 ! TS 18 Sep 2017 - added explicit interface
421 ! HCW 08 Feb 2017 - fixed bug in Zh between grids, added default z0m, zdm
422 ! HCW 03 Mar 2015
423 ! sg feb 2012 - made separate subroutine
424 !--------------------------------------------------------------------------------
425
426 IMPLICIT NONE
427
428 INTEGER, PARAMETER :: nsurf = 7 ! number of surface types
429 INTEGER, PARAMETER :: PavSurf = 1 !When all surfaces considered together (1-7)
430 INTEGER, PARAMETER :: BldgSurf = 2
431 INTEGER, PARAMETER :: ConifSurf = 3
432 INTEGER, PARAMETER :: DecidSurf = 4
433 INTEGER, PARAMETER :: GrassSurf = 5 !New surface classes: Grass = 5th/7 surfaces
434 INTEGER, PARAMETER :: BSoilSurf = 6 !New surface classes: Bare soil = 6th/7 surfaces
435 INTEGER, PARAMETER :: WaterSurf = 7
436 REAL(KIND(1D0)), PARAMETER :: porosity_evetr = 0.32 ! assumed porosity of evergreen trees, ref: Lai et al. (2022), http://dx.doi.org/10.2139/ssrn.4058842
437
438 INTEGER, INTENT(in) :: RoughLenMomMethod
439 INTEGER, INTENT(in) :: FAImethod ! 0 = use FAI provided, 1 = use the simple scheme
440
441 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! surface fractions
442
443 REAL(KIND(1D0)), INTENT(in) :: surfaceArea ! surface area of whole grid cell
444 REAL(KIND(1D0)), INTENT(in) :: bldgH
445 REAL(KIND(1D0)), INTENT(in) :: EveTreeH
446 REAL(KIND(1D0)), INTENT(in) :: DecTreeH
447 REAL(KIND(1D0)), INTENT(in) :: porosity_dectr
448 REAL(KIND(1D0)), INTENT(in) :: FAIBldg
449 REAL(KIND(1D0)), INTENT(in) :: FAIEveTree
450 REAL(KIND(1D0)), INTENT(in) :: FAIDecTree
451 REAL(KIND(1D0)), INTENT(in) :: z0m_in ! z0m set in SiteSelect
452 REAL(KIND(1D0)), INTENT(in) :: zdm_in ! zdm set in SiteSelect
453 REAL(KIND(1D0)), INTENT(in) :: Z
454
455 REAL(KIND(1D0)), INTENT(out) :: FAI
456 REAL(KIND(1D0)), INTENT(out) :: PAI
457 REAL(KIND(1D0)), INTENT(out) :: Zh ! effective height of bluff bodies
458 REAL(KIND(1D0)), INTENT(out) :: z0m ! aerodynamic roughness length
459 REAL(KIND(1D0)), INTENT(out) :: zdm ! zero-plance displacement
460 REAL(KIND(1D0)), INTENT(out) :: ZZD ! z-zdm
461
462 INTEGER, PARAMETER :: notUsedI = -55
463 REAL(KIND(1D0)), PARAMETER :: notUsed = -55.5
464 REAL(KIND(1D0)) :: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water !Default values for roughness lengths [m]
465
466 ! calculated values of FAI
467 REAL(KIND(1D0)) :: FAIBldg_use
468 REAL(KIND(1D0)) :: FAIEveTree_use
469 REAL(KIND(1D0)) :: FAIDecTree_use
470
471 !Total area of buildings and trees
472 ! areaZh = (sfr_surf(BldgSurf) + sfr_surf(ConifSurf) + sfr_surf(DecidSurf))
473 ! TS 19 Jun 2022: take porosity of trees into account; to be consistent with PAI calculation in RSL
474 pai = dot_product(sfr_surf([bldgsurf, conifsurf, decidsurf]), [1d0, 1 - porosity_evetr, 1 - porosity_dectr])
475
476 ! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
477 z0m4paved = 0.003 !estimate
478 z0m4grass = 0.02
479 z0m4bsoil = 0.002
480 z0m4water = 0.0005
481
482 !------------------------------------------------------------------------------
483 !If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
484 IF (pai > 0.02) THEN ! 0.02 is a threshold value, below which we assume no buildings and trees, TS 06 Jun 2023
485 zh = dot_product( &
486 [bldgh, evetreeh*(1 - porosity_evetr), dectreeh*(1 - porosity_dectr)], &
487 sfr_surf([bldgsurf, conifsurf, decidsurf]))/pai
488 IF (faimethod == 0) THEN
489 ! use FAI provided
490 faibldg_use = faibldg
491 faievetree_use = faievetree*(1 - porosity_evetr)
492 faidectree_use = faidectree*(1 - porosity_dectr)
493
494 ELSEIF (faimethod == 1) THEN
495 ! use the simple scheme, details in #192
496
497 ! buildings
498 faibldg_use = sqrt(sfr_surf(bldgsurf)/surfacearea)*bldgh
499
500 ! evergreen trees
501 faievetree_use = 1.07*sfr_surf(conifsurf)
502
503 ! deciduous trees
504 faidectree_use = 1.66*(1 - porosity_dectr)*sfr_surf(decidsurf)
505
506 END IF
507 fai = sum(merge([faibldg_use, faievetree_use, faidectree_use], &
508 [0d0, 0d0, 0d0], &
509 sfr_surf([bldgsurf, conifsurf, decidsurf]) > 0))
510 ! `1e-5` set to avoid numerical difficulty
511 fai = max(fai, 1e-5)
512
513 ELSE
514 zh = 0 !Set Zh to zero if areaZh = 0
515 fai = 1e-5
516 END IF
517
518 IF (zh /= 0) THEN
519 !Calculate z0m and zdm depending on the Z0 method
520 IF (roughlenmommethod == 2) THEN !Rule of thumb (G&O 1999)
521 z0m = 0.1*zh
522 zdm = 0.7*zh
523 ELSEIF (roughlenmommethod == 3) THEN !MacDonald 1998
524 ! zdm = (1 + 4.43**(-sfr_surf(BldgSurf))*(sfr_surf(BldgSurf) - 1))*Zh
525 zdm = (1 + 4.43**(-pai)*(pai - 1))*zh ! changed sfr_surf(BldgSurf) to PAI to be consistent with z0m calculation, TS 04 Jun 2023
526 z0m = ((1 - zdm/zh)*exp(-(0.5*1.0*1.2/0.4**2*(1 - zdm/zh)*fai)**(-0.5)))*zh
527 ELSEIF (roughlenmommethod == 4) THEN ! lambdaP dependent as in Fig.1a of G&O (1999)
528 ! these are derived using digitalised points
529 zdm = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*pai) + 0.493*sigmoid(-5.17 + 32.7*pai))*zh
530 z0m = (0.00208 + &
531 0.0165*min(pai, .7) + 2.52*min(pai, .7)**2 + &
532 3.21*min(pai, .7)**3 - 43.6*min(pai, .7)**4 + &
533 76.5*min(pai, .7)**5 - 40.*min(pai, .7)**6)*zh
534 END IF
535 ELSEIF (zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
536 IF (pai /= 0) CALL errorhint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, pai, notusedi)
537 !Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
538 IF (pai /= 1) THEN
539 z0m = (z0m4paved*sfr_surf(pavsurf) &
540 + z0m4grass*sfr_surf(grasssurf) &
541 + z0m4bsoil*sfr_surf(bsoilsurf) &
542 + z0m4water*sfr_surf(watersurf))/(1 - pai)
543 zdm = 0
544 CALL errorhint(15, 'Setting z0m and zdm using default values', z0m, zdm, notusedi)
545 ELSEIF (pai == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
546 z0m = 1
547 zdm = 7
548 CALL errorhint(15, 'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
549 END IF
550 END IF
551
552 IF (roughlenmommethod == 1) THEN !use values set in SiteSelect
553 z0m = z0m_in
554 zdm = zdm_in
555 END IF
556
557 zzd = z - zdm
558
559 ! Error messages if aerodynamic parameters negative
560 IF (z0m < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, z0 < 0 m.', z0m, notused, notusedi)
561 IF (zdm < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, zd < 0 m.', zdm, notused, notusedi)
562 IF (zzd < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, (z-zd) < 0 m.', zzd, notused, notusedi)

References errorhint(), and sigmoid().

Referenced by suews_driver::suews_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ suews_cal_roughnessparameters_dts()

subroutine resist_module::suews_cal_roughnessparameters_dts ( integer, intent(in) roughlenmommethod,
integer, intent(in) faimethod,
real(kind(1d0)), intent(in) sfr_paved,
real(kind(1d0)), intent(in) sfr_bldg,
real(kind(1d0)), intent(in) sfr_evetr,
real(kind(1d0)), intent(in) sfr_dectr,
real(kind(1d0)), intent(in) sfr_grass,
real(kind(1d0)), intent(in) sfr_bsoil,
real(kind(1d0)), intent(in) sfr_water,
real(kind(1d0)), intent(in) surfacearea,
real(kind(1d0)), intent(in) bldgh,
real(kind(1d0)), intent(in) evetreeh,
real(kind(1d0)), intent(in) dectreeh,
real(kind(1d0)), intent(in) porosity_dectr,
real(kind(1d0)), intent(in) faibldg,
real(kind(1d0)), intent(in) faievetree,
real(kind(1d0)), intent(in) faidectree,
real(kind(1d0)), intent(in) z0m_in,
real(kind(1d0)), intent(in) zdm_in,
real(kind(1d0)), intent(in) z,
real(kind(1d0)), intent(out) faibldg_use,
real(kind(1d0)), intent(out) faievetree_use,
real(kind(1d0)), intent(out) faidectree_use,
real(kind(1d0)), intent(out) fai,
real(kind(1d0)), intent(out) pai,
real(kind(1d0)), intent(out) zh,
real(kind(1d0)), intent(out) z0m,
real(kind(1d0)), intent(out) zdm,
real(kind(1d0)), intent(out) zzd )

Definition at line 565 of file suews_phys_resist.f95.

575 ! Get surface covers and frontal area fractions (LJ 11/2010)
576 ! Last modified:
577 ! TS 18 Sep 2017 - added explicit interface
578 ! HCW 08 Feb 2017 - fixed bug in Zh between grids, added default z0m, zdm
579 ! HCW 03 Mar 2015
580 ! sg feb 2012 - made separate subroutine
581 !--------------------------------------------------------------------------------
582
583 IMPLICIT NONE
584
585 INTEGER, PARAMETER :: nsurf = 7 ! number of surface types
586 INTEGER, PARAMETER :: PavSurf = 1 !When all surfaces considered together (1-7)
587 INTEGER, PARAMETER :: BldgSurf = 2
588 INTEGER, PARAMETER :: ConifSurf = 3
589 INTEGER, PARAMETER :: DecidSurf = 4
590 INTEGER, PARAMETER :: GrassSurf = 5 !New surface classes: Grass = 5th/7 surfaces
591 INTEGER, PARAMETER :: BSoilSurf = 6 !New surface classes: Bare soil = 6th/7 surfaces
592 INTEGER, PARAMETER :: WaterSurf = 7
593 REAL(KIND(1D0)), PARAMETER :: porosity_evetr = 0.32 ! assumed porosity of evergreen trees, ref: Lai et al. (2022), http://dx.doi.org/10.2139/ssrn.4058842
594
595 INTEGER, INTENT(in) :: RoughLenMomMethod
596 INTEGER, INTENT(in) :: FAImethod ! 0 = use FAI provided, 1 = use the simple scheme
597
598 REAL(KIND(1D0)), INTENT(IN) :: sfr_paved
599 REAL(KIND(1D0)), INTENT(IN) :: sfr_bldg
600 REAL(KIND(1D0)), INTENT(IN) :: sfr_evetr
601 REAL(KIND(1D0)), INTENT(IN) :: sfr_dectr
602 REAL(KIND(1D0)), INTENT(IN) :: sfr_grass
603 REAL(KIND(1D0)), INTENT(IN) :: sfr_bsoil
604 REAL(KIND(1D0)), INTENT(IN) :: sfr_water
605 REAL(KIND(1D0)), DIMENSION(nsurf) :: sfr_surf ! surface fractions
606
607 REAL(KIND(1D0)), INTENT(in) :: surfaceArea ! surface area of whole grid cell
608 REAL(KIND(1D0)), INTENT(in) :: bldgH
609 REAL(KIND(1D0)), INTENT(in) :: EveTreeH
610 REAL(KIND(1D0)), INTENT(in) :: DecTreeH
611 REAL(KIND(1D0)), INTENT(in) :: porosity_dectr
612 REAL(KIND(1D0)), INTENT(in) :: FAIBldg
613 REAL(KIND(1D0)), INTENT(in) :: FAIEveTree
614 REAL(KIND(1D0)), INTENT(in) :: FAIDecTree
615 REAL(KIND(1D0)), INTENT(in) :: z0m_in ! z0m set in SiteSelect
616 REAL(KIND(1D0)), INTENT(in) :: zdm_in ! zdm set in SiteSelect
617 REAL(KIND(1D0)), INTENT(in) :: Z
618
619 REAL(KIND(1D0)), INTENT(out) :: FAI
620 REAL(KIND(1D0)), INTENT(out) :: PAI
621 REAL(KIND(1D0)), INTENT(out) :: Zh ! effective height of bluff bodies
622 REAL(KIND(1D0)), INTENT(out) :: z0m ! aerodynamic roughness length
623 REAL(KIND(1D0)), INTENT(out) :: zdm ! zero-plance displacement
624 REAL(KIND(1D0)), INTENT(out) :: ZZD ! z-zdm
625
626 INTEGER, PARAMETER :: notUsedI = -55
627 REAL(KIND(1D0)), PARAMETER :: notUsed = -55.5
628 REAL(KIND(1D0)) :: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water !Default values for roughness lengths [m]
629
630 ! calculated values of FAI
631 REAL(KIND(1D0)), INTENT(out) :: FAIBldg_use
632 REAL(KIND(1D0)), INTENT(out) :: FAIEveTree_use
633 REAL(KIND(1D0)), INTENT(out) :: FAIDecTree_use
634
635 sfr_surf = [sfr_paved, sfr_bldg, sfr_evetr, sfr_dectr, sfr_grass, sfr_bsoil, sfr_water]
636
637 !Total area of buildings and trees
638 ! areaZh = (sfr_surf(BldgSurf) + sfr_surf(ConifSurf) + sfr_surf(DecidSurf))
639 ! TS 19 Jun 2022: take porosity of trees into account; to be consistent with PAI calculation in RSL
640 pai = dot_product(sfr_surf([bldgsurf, conifsurf, decidsurf]), [1d0, 1 - porosity_evetr, 1 - porosity_dectr])
641
642 ! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
643 z0m4paved = 0.003 !estimate
644 z0m4grass = 0.02
645 z0m4bsoil = 0.002
646 z0m4water = 0.0005
647
648 !------------------------------------------------------------------------------
649 !If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
650 IF (pai /= 0) THEN
651 zh = dot_product( &
652 [bldgh, evetreeh*(1 - porosity_evetr), dectreeh*(1 - porosity_dectr)], &
653 sfr_surf([bldgsurf, conifsurf, decidsurf]))/pai
654 IF (faimethod == 0) THEN
655 ! use FAI provided
656 faibldg_use = faibldg
657 faievetree_use = faievetree*(1 - porosity_evetr)
658 faidectree_use = faidectree*(1 - porosity_dectr)
659
660 ELSEIF (faimethod == 1) THEN
661 ! use the simple scheme, details in #192
662
663 ! buildings
664 faibldg_use = sqrt(sfr_surf(bldgsurf)/surfacearea)*bldgh
665
666 ! evergreen trees
667 faievetree_use = 1.07*sfr_surf(conifsurf)
668
669 ! deciduous trees
670 faidectree_use = 1.66*(1 - porosity_dectr)*sfr_surf(decidsurf)
671
672 END IF
673 fai = sum(merge([faibldg_use, faievetree_use, faidectree_use], &
674 [0d0, 0d0, 0d0], &
675 sfr_surf([bldgsurf, conifsurf, decidsurf]) > 0))
676
677 ! `1e-5` set to avoid numerical difficulty
678 fai = max(fai, 1e-5)
679
680 ELSE
681 zh = 0 !Set Zh to zero if areaZh = 0
682 fai = 1e-5
683 END IF
684
685 IF (zh /= 0) THEN
686 !Calculate z0m and zdm depending on the Z0 method
687 IF (roughlenmommethod == 2) THEN !Rule of thumb (G&O 1999)
688 z0m = 0.1*zh
689 zdm = 0.7*zh
690 ELSEIF (roughlenmommethod == 3) THEN !MacDonald 1998
691 zdm = (1 + 4.43**(-sfr_surf(bldgsurf))*(sfr_surf(bldgsurf) - 1))*zh
692 z0m = ((1 - zdm/zh)*exp(-(0.5*1.0*1.2/0.4**2*(1 - zdm/zh)*fai)**(-0.5)))*zh
693 ELSEIF (roughlenmommethod == 4) THEN ! lambdaP dependent as in Fig.1a of G&O (1999)
694 ! these are derived using digitalised points
695 zdm = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*pai) + 0.493*sigmoid(-5.17 + 32.7*pai))*zh
696 z0m = (0.00208 + &
697 0.0165*min(pai, .7) + 2.52*min(pai, .7)**2 + &
698 3.21*min(pai, .7)**3 - 43.6*min(pai, .7)**4 + &
699 76.5*min(pai, .7)**5 - 40.*min(pai, .7)**6)*zh
700 END IF
701 ELSEIF (zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
702 IF (pai /= 0) CALL errorhint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, pai, notusedi)
703 !Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
704 IF (pai /= 1) THEN
705 z0m = (z0m4paved*sfr_surf(pavsurf) &
706 + z0m4grass*sfr_surf(grasssurf) &
707 + z0m4bsoil*sfr_surf(bsoilsurf) &
708 + z0m4water*sfr_surf(watersurf))/(1 - pai)
709 zdm = 0
710 CALL errorhint(15, 'Setting z0m and zdm using default values', z0m, zdm, notusedi)
711 ELSEIF (pai == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
712 z0m = 1
713 zdm = 7
714 CALL errorhint(15, 'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
715 END IF
716 END IF
717
718 IF (roughlenmommethod == 1) THEN !use values set in SiteSelect
719 z0m = z0m_in
720 zdm = zdm_in
721 END IF
722
723 zzd = z - zdm
724
725 ! Error messages if aerodynamic parameters negative
726 IF (z0m < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, z0 < 0 m.', z0m, notused, notusedi)
727 IF (zdm < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, zd < 0 m.', zdm, notused, notusedi)
728 IF (zzd < 0) CALL errorhint(14, 'In SUEWS_cal_RoughnessParameters, (z-zd) < 0 m.', zzd, notused, notusedi)

References errorhint(), and sigmoid().

Referenced by suews_driver::suews_cal_main_dts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ surfaceresistance()

subroutine resist_module::surfaceresistance ( integer, intent(in) id,
integer, intent(in) it,
integer, intent(in) smdmethod,
real(kind(1d0)), dimension(nsurf), intent(in) snowfrac,
real(kind(1d0)), dimension(nsurf), intent(in) sfr_surf,
real(kind(1d0)), intent(in) avkdn,
real(kind(1d0)), intent(in) temp_c,
real(kind(1d0)), intent(in) dq,
real(kind(1d0)), intent(in) xsmd,
real(kind(1d0)), intent(in) vsmd,
real(kind(1d0)), dimension(3), intent(in) maxconductance,
real(kind(1d0)), dimension(3), intent(in) laimax,
real(kind(1d0)), dimension(3), intent(in) lai_id,
integer, intent(in) gsmodel,
real(kind(1d0)), intent(in) kmax,
real(kind(1d0)), intent(in) g_max,
real(kind(1d0)), intent(in) g_k,
real(kind(1d0)), intent(in) g_q_base,
real(kind(1d0)), intent(in) g_q_shape,
real(kind(1d0)), intent(in) g_t,
real(kind(1d0)), intent(in) g_sm,
real(kind(1d0)), intent(in) th,
real(kind(1d0)), intent(in) tl,
real(kind(1d0)), intent(in) s1,
real(kind(1d0)), intent(in) s2,
real(kind(1d0)), intent(out) g_kdown,
real(kind(1d0)), intent(out) g_dq,
real(kind(1d0)), intent(out) g_ta,
real(kind(1d0)), intent(out) g_smd,
real(kind(1d0)), intent(out) g_lai,
real(kind(1d0)), intent(out) gfunc,
real(kind(1d0)), intent(out) gsc,
real(kind(1d0)), intent(out) rs )

Definition at line 111 of file suews_phys_resist.f95.

118 ! Calculates bulk surface resistance (ResistSurf [s m-1]) based on Jarvis 1976 approach
119 ! Last modified -----------------------------------------------------
120 ! MH 01 Feb 2019: gsModel choices to model with air temperature or 2 meter temperature. Added gfunc for photosynthesis calculations
121 ! HCW 21 Jul 2016: If no veg surfaces, vsmd = NaN so QE & QH = NaN; if water surfaces only, smd = NaN so QE & QH = NaN.
122 ! Add checks here so that gs (soil part) = 0 in either of these situations.
123 ! This shouldn't change results but handles NaN error.
124 ! HCW 01 Mar 2016: SM dependence is now on modelled smd for vegetated surfaces only (vsmd) (Note: obs smd still not operational!)
125 ! HCW 18 Jun 2015: Alternative gs parameterisation added using different functional forms and new coefficients
126 ! HCW 31 Jul 2014: Modified condition on g6 part to select meas/mod smd
127 ! LJ 24 Apr 2013: Added impact of snow fraction in LAI and in soil moisture deficit
128 ! -------------------------------------------------------------------
129
130 ! USE allocateArray
131 ! USE data_in
132 ! USE defaultNotUsed
133 ! USE gis_data
134 ! USE moist
135 ! USE resist
136 ! USE sues_data
137
138 IMPLICIT NONE
139 ! INTEGER,PARAMETER::BldgSurf=2
140 INTEGER, PARAMETER :: ConifSurf = 3
141 INTEGER, PARAMETER :: DecidSurf = 4
142 INTEGER, PARAMETER :: GrassSurf = 5
143 ! INTEGER,PARAMETER::ivConif=1
144 ! INTEGER,PARAMETER::ivGrass=3
145 ! INTEGER,PARAMETER::MaxNumberOfGrids=2000
146 ! INTEGER,PARAMETER::ndays=366
147 INTEGER, PARAMETER :: nsurf = 7
148 ! INTEGER,PARAMETER::NVegSurf=3
149 ! INTEGER,PARAMETER::PavSurf=1
150 INTEGER, PARAMETER :: WaterSurf = 7
151
152 INTEGER, INTENT(in) :: id
153 INTEGER, INTENT(in) :: it ! time: day of year and hour
154 INTEGER, INTENT(in) :: gsModel !Choice of gs parameterisation (1 = Ja11, 2 = Wa16)
155 INTEGER, INTENT(in) :: SMDMethod !Method of measured soil moisture
156 ! INTEGER,INTENT(in)::ConifSurf!= 3, surface code
157 ! INTEGER,INTENT(in)::DecidSurf!= 4, surface code
158 ! INTEGER,INTENT(in)::GrassSurf!= 5, surface code
159 ! INTEGER,INTENT(in)::WaterSurf!= 7, surface code
160 ! INTEGER,INTENT(in)::nsurf!= 7, Total number of surfaces
161
162 REAL(KIND(1D0)), INTENT(in) :: avkdn !Average downwelling shortwave radiation
163 REAL(KIND(1D0)), INTENT(in) :: Temp_C !Air temperature
164 REAL(KIND(1D0)), INTENT(in) :: Kmax !Annual maximum hourly solar radiation
165 REAL(KIND(1D0)), INTENT(in) :: G_max !Fitted parameters related to max surface resistance
166 REAL(KIND(1D0)), INTENT(in) :: G_k !Fitted parameters related to incoming solar radiation
167 REAL(KIND(1D0)), INTENT(in) :: g_q_base !Fitted parameters related to humidity deficit, base value
168 REAL(KIND(1D0)), INTENT(in) :: g_q_shape !Fitted parameters related to humidity deficit, shape parameter
169 REAL(KIND(1D0)), INTENT(in) :: G_t !Fitted parameters related to air temperature
170 REAL(KIND(1D0)), INTENT(in) :: G_sm !Fitted parameters related to soil moisture deficit
171 REAL(KIND(1D0)), INTENT(in) :: S1 !Fitted parameters related to wilting point
172 REAL(KIND(1D0)), INTENT(in) :: S2 !Fitted parameters related to wilting point
173 REAL(KIND(1D0)), INTENT(in) :: TH !Maximum temperature limit
174 REAL(KIND(1D0)), INTENT(in) :: TL !Minimum temperature limit
175 REAL(KIND(1D0)), INTENT(in) :: dq !Specific humidity deficit
176 REAL(KIND(1D0)), INTENT(in) :: xsmd !Measured soil moisture deficit
177 REAL(KIND(1D0)), INTENT(in) :: vsmd !QUESTION: Soil moisture deficit for vegetated surfaces only (what about BSoil?)
178
179 REAL(KIND(1D0)), DIMENSION(3), INTENT(in) :: MaxConductance !Max conductance [mm s-1]
180 REAL(KIND(1D0)), DIMENSION(3), INTENT(in) :: LAIMax !Max LAI [m2 m-2]
181 REAL(KIND(1D0)), DIMENSION(3), INTENT(in) :: LAI_id !=LAI(id-1,:), LAI for each veg surface [m2 m-2]
182 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowFrac !Surface fraction of snow cover
183 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf !Surface fractions [-]
184 ! g_kdown*g_dq*g_ta*g_smd*g_lai
185 REAL(KIND(1D0)), INTENT(out) :: g_kdown !gdq*gtemp*gs*gq for photosynthesis calculations
186 REAL(KIND(1D0)), INTENT(out) :: g_dq !gdq*gtemp*gs*gq for photosynthesis calculations
187 REAL(KIND(1D0)), INTENT(out) :: g_ta !gdq*gtemp*gs*gq for photosynthesis calculations
188 REAL(KIND(1D0)), INTENT(out) :: g_smd !gdq*gtemp*gs*gq for photosynthesis calculations
189 REAL(KIND(1D0)), INTENT(out) :: g_lai !gdq*gtemp*gs*gq for photosynthesis calculations
190 REAL(KIND(1D0)), INTENT(out) :: gfunc !gdq*gtemp*gs*gq for photosynthesis calculations
191 REAL(KIND(1D0)), INTENT(out) :: gsc !Surface Layer Conductance
192 REAL(KIND(1D0)), INTENT(out) :: RS !Surface resistance
193
194 REAL(KIND(1D0)) :: &
195 ! g_lai, & !G(LAI)
196 qnm, & !QMAX/(QMAX+G2)
197 ! g_kdown, & !G(Q*)
198 ! g_dq, & !G(dq)
199 tc, & !Temperature parameter 1
200 tc2, & !Temperature parameter 2
201 ! g_ta, & !G(T)
202 sdp !S1/G6+S2
203 ! g_smd !G(Soil moisture deficit)
204
205 INTEGER :: iv
206 REAL(KIND(1D0)) :: id_real
207
208 REAL(KIND(1D0)), PARAMETER :: notUsed = -55
209 ! REAL(KIND(1d0)),PARAMETER :: notUsedi=-55.5
210
211 ! initialisation
212 g_dq = 0.5
213 g_ta = 0.5
214 g_smd = 0.5
215 g_kdown = 0.5
216
217 id_real = real(id) !Day of year in real number
218
219 !gsModel = 1 - original parameterisation (Jarvi et al. 2011)
220 !gsModel = 2 - new parameterisation (Ward et al. 2016)
221 !gsModel = 3 - original parameterisation (Jarvi et al. 2011) with 2 m temperature
222 !gsModel = 4 - new parameterisation (Ward et al. 2016) with 2 m temperature
223
224 IF (gsmodel == 1 .OR. gsmodel == 3) THEN
225 ! kdown ----
226 qnm = kmax/(kmax + g_k)
227 !gq=(qn1/(g2+qn1))/qnm !With net all-wave radiation
228 g_kdown = (avkdn/(g_k + avkdn))/qnm !With Kdown
229 ! specific humidity deficit ----
230 IF (dq < g_q_shape) THEN
231 g_dq = 1 - g_q_base*dq
232 ELSE
233 g_dq = 1 - g_q_base*g_q_shape
234 END IF
235 ! air temperature ----
236 tc = (th - g_t)/(g_t - tl)
237 tc2 = (g_t - tl)*(th - g_t)**tc
238 !If air temperature below TL or above TH, fit it to TL+0.1/TH-0.1
239 IF (temp_c <= tl) THEN
240 g_ta = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
241 !Call error only if no snow on ground
242 ! IF (MIN(SnowFrac(1),SnowFrac(2),SnowFrac(3),SnowFrac(4),SnowFrac(5),SnowFrac(6))/=1) THEN
243 IF (minval(snowfrac(1:6)) /= 1) THEN
244 CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TL=0.1,Temp_c,id,it', &
245 REAL(Temp_c, KIND(1D0)), id_real, it)
246 END IF
247 ELSEIF (temp_c >= th) THEN
248 g_ta = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
249 CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TH=39.9,Temp_c,id,it', &
250 REAL(Temp_c, KIND(1D0)), id_real, it)
251 ELSE
252 g_ta = (temp_c - tl)*(th - temp_c)**tc/tc2
253 END IF
254
255 ! soil moisture deficit ----
256 sdp = s1/g_sm + s2
257 IF (smdmethod > 0) THEN !Modified from ==1 to > 0 by HCW 31/07/2014
258 g_smd = 1 - exp(g_sm*(xsmd - sdp)) !Measured soil moisture deficit is used
259 ELSE
260 g_smd = 1 - exp(g_sm*(vsmd - sdp)) !Modelled is used
261 IF (sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf) == 0 .OR. sfr_surf(watersurf) == 1) THEN
262 g_smd = 0 !If no veg so no vsmd, or all water so no smd, set gs=0 (HCW 21 Jul 2016)
263 END IF
264 END IF
265
266 IF (g_smd < 0) THEN
267 CALL errorhint(65, &
268 'subroutine SurfaceResistance.f95 (gsModel=1): g(smd) < 0 calculated, setting to 0.0001', &
269 g_smd, id_real, it)
270 g_smd = 0.0001
271 END IF
272
273 !LAI
274 !Original way
275 !gl=((LAI(id,2)*areaunir/lm)+areair)/(areair+areaunir)
276 !New way
277 g_lai = 0 !First initialize
278 ! vegetated surfaces
279 ! check basis for values koe - maximum surface conductance
280 ! print*,id,it,sfr_surf
281 ! DO iv=ivConif,ivGrass
282 DO iv = 1, 3
283 ! gl=gl+(sfr_surf(iv+2)*(1-SnowFrac(iv+2)))*LAI(id-1,iv)/LAIMax(iv)*MaxConductance(iv)
284 g_lai = g_lai + (sfr_surf(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
285 END DO
286
287 IF (avkdn <= 0) THEN !At nighttime set gsc at arbitrary low value: gsc=0.1 mm/s (Shuttleworth, 1988b)
288 gsc = 0.1
289 ELSE
290 ! Multiply parts together
291 gsc = (g_max*g_kdown*g_dq*g_ta*g_smd*g_lai)
292 END IF
293
294 IF (gsc <= 0) THEN
295 CALL errorhint(65, 'subroutine SurfaceResistance.f95 (gsModel=1): gs <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
296 gsc = 0.1
297 END IF
298
299 ELSEIF (gsmodel == 2 .OR. gsmodel == 4) THEN
300
301 ! ---- g(kdown)----
302 qnm = kmax/(kmax + g_k)
303 g_kdown = (avkdn/(avkdn + g_k))/qnm
304 IF (avkdn >= kmax) THEN !! Add proper error handling later - HCW!!
305 WRITE (*, *) 'Kmax exceeds Kdn setting to g(Kdn) to 1'
306 g_kdown = 1
307 END IF
308
309 ! ---- g(delq) ----
310 g_dq = g_q_base + (1 - g_q_base)*(g_q_shape**dq) !Ogink-Hendriks (1995) Eq 12 (using G3 as Kshd and G4 as r)
311
312 ! ---- g(Tair) ----
313 tc = (th - g_t)/(g_t - tl)
314 tc2 = (g_t - tl)*(th - g_t)**tc
315 ! If air temperature below TL or above TH, then use value for TL+0.1 or TH-0.1
316 IF (temp_c <= tl) THEN
317 g_ta = (tl + 0.1 - tl)*(th - (tl + 0.1))**tc/tc2
318 ! Call error only if no snow on ground
319 IF (min(snowfrac(1), snowfrac(2), snowfrac(3), snowfrac(4), snowfrac(5), snowfrac(6)) /= 1) THEN
320 CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TL+0.1,Temp_C,id,it', &
321 REAL(Temp_c, KIND(1D0)), id_real, it)
322 END IF
323 ELSEIF (temp_c >= th) THEN
324 g_ta = ((th - 0.1) - tl)*(th - (th - 0.1))**tc/tc2
325 CALL errorhint(29, 'subroutine SurfaceResistance.f95: T changed to fit limits TH-0.1,Temp_C,id,it', &
326 REAL(Temp_c, KIND(1D0)), id_real, it)
327 ELSE
328 g_ta = (temp_c - tl)*(th - temp_c)**tc/tc2
329 END IF
330 ! ---- g(smd) ----
331 sdp = s1/g_sm + s2
332 IF (smdmethod > 0) THEN !Modified from ==1 to > 0 by HCW 31/07/2014
333 g_smd = (1 - exp(g_sm*(xsmd - sdp)))/(1 - exp(g_sm*(-sdp))) !Use measured smd
334 ELSE
335 !gs=1-EXP(g6*(vsmd-sdp)) !Use modelled smd
336 g_smd = (1 - exp(g_sm*(vsmd - sdp)))/(1 - exp(g_sm*(-sdp)))
337 IF (sfr_surf(conifsurf) + sfr_surf(decidsurf) + sfr_surf(grasssurf) == 0 .OR. sfr_surf(watersurf) == 1) THEN
338 g_smd = 0 !If no veg so no vsmd, or all water so no smd, set gs=0 HCW 21 Jul 2016
339 END IF
340 END IF
341
342 !gs = gs*(1 - SUM(SnowFrac(1:6))/6)
343
344 IF (g_smd < 0) THEN
345 CALL errorhint(65, &
346 'subroutine SurfaceResistance.f95 (gsModel=2): gs < 0 calculated, setting to 0.0001', &
347 g_smd, id_real, it)
348 g_smd = 0.0001
349 END IF
350
351 ! ---- g(LAI) ----
352 g_lai = 0 !Initialise
353 ! DO iv=ivConif,ivGrass !For vegetated surfaces
354 DO iv = 1, 3 !For vegetated surfaces
355 ! gl=gl+(sfr_surf(iv+2)*(1-SnowFrac(iv+2)))*LAI(id-1,iv)/LAIMax(iv)*MaxConductance(iv)
356 g_lai = g_lai + (sfr_surf(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)/laimax(iv)*maxconductance(iv)
357 END DO
358
359 IF (avkdn <= 0) THEN !At nighttime set gsc at arbitrary low value: gsc=0.1 mm/s (Shuttleworth, 1988b)
360 gsc = 0.1
361 ELSE
362 ! Multiply parts together
363 gsc = (g_max*g_kdown*g_dq*g_ta*g_smd*g_lai)
364 END IF
365
366 IF (gsc <= 0) THEN
367 CALL errorhint(65, 'subroutine SurfaceResistance.f95 (gsModel=2): gsc <= 0, setting to 0.1 mm s-1', gsc, id_real, it)
368 gsc = 0.1
369 END IF
370
371 ELSEIF (gsmodel < 1 .OR. gsmodel > 4) THEN
372 CALL errorhint(71, 'Value of gsModel not recognised.', notused, notused, gsmodel)
373 END IF
374
375 rs = 1/(gsc/1000) ![s m-1]
376 gfunc = g_dq*g_ta*g_smd*g_kdown
377
378 RETURN

References errorhint().

Referenced by suews_driver::suews_cal_biogenco2(), suews_driver::suews_cal_biogenco2_dts(), suews_driver::suews_cal_resistance(), and suews_driver::suews_cal_resistance_dts().

Here is the call graph for this function:
Here is the caller graph for this function:
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0