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, G1, G2, G3, G4, G5, G6, TH, TL, S1, S2, gfunc, gsc, ResistSurf)
 
subroutine boundarylayerresistance (zzd, z0m, avU1, UStar, rb)
 
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)
 
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
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)
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 > 200) 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 = 200
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
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
integer stabilitymethod
integer aerodynamicresistancemethod
real(kind(1d0)) ustar
real(kind(1d0)) l_mod
real(kind(1d0)) psih
integer roughlenheatmethod
real(kind(1d0)) psim
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)

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().

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 375 of file suews_phys_resist.f95.

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

Referenced by suews_driver::suews_cal_resistance().

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 532 of file suews_phys_resist.f95.

533 ! TS 31 Jul 2018: make this a separate funciton for reuse
534 !Z0V roughness length for vapour
535 IMPLICIT NONE
536 INTEGER, INTENT(in) :: RoughLenHeatMethod
537 REAL(KIND(1D0)), INTENT(in) :: z0m !Aerodynamic roughness length
538 REAL(KIND(1D0)), INTENT(in) :: VegFraction !Fraction of vegetation
539 REAL(KIND(1D0)), INTENT(in) :: UStar !Friction velocity
540
541 REAL(KIND(1D0)) :: z0V !roughness length for vapor/heat
542
543 REAL(KIND(1D0)), PARAMETER :: muu = 1.46e-5 !molecular viscosity
544
545 z0v = 0.01 ! initialise as 0.01
546
547 !Z0V roughness length for vapour
548 IF (roughlenheatmethod == 1) THEN !Brutsaert (1982) Z0v=z0/10(see Grimmond & Oke, 1986)
549 z0v = z0m/10
550 ELSEIF (roughlenheatmethod == 2) THEN ! Kawai et al. (2007)
551 !z0V=z0m*exp(2-(1.2-0.9*veg_fr**0.29)*(UStar*z0m/muu)**0.25)
552 ! Changed by HCW 05 Nov 2015 (veg_fr includes water; VegFraction = veg + bare soil)
553 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
554 ELSEIF (roughlenheatmethod == 3) THEN
555 z0v = z0m*exp(-20.) ! Voogt and Grimmond, JAM, 2000
556 ELSEIF (roughlenheatmethod == 4) THEN
557 z0v = z0m*exp(2 - 1.29*(ustar*z0m/muu)**0.25) !See !Kanda and Moriwaki (2007),Loridan et al. (2010)
558 ELSEIF (roughlenheatmethod == 5) THEN
559 ! an adaptive way to determine z0v; TS 06 Feb 2020
560 IF (vegfraction > .999) THEN
561 ! fully pervious surface
562 z0v = z0m/10
563 ELSE
564 ! impervious-pervious mixed surface
565 z0v = z0m*exp(2 - (1.2 - 0.9*vegfraction**0.29)*(ustar*z0m/muu)**0.25)
566 END IF
567 END IF
568

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 571 of file suews_phys_resist.f95.

572 IMPLICIT NONE
573 REAL(KIND(1D0)), INTENT(in) :: x
574 REAL(KIND(1D0)) :: res
575
576 res = 1/(1 + exp(-x))
577

Referenced by suews_cal_roughnessparameters().

Here is the caller graph for this function:

◆ suews_cal_roughnessparameters()

subroutine resist_module::suews_cal_roughnessparameters ( integer, intent(in)  RoughLenMomMethod,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
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 403 of file suews_phys_resist.f95.

411 ! Get surface covers and frontal area fractions (LJ 11/2010)
412 ! Last modified:
413 ! TS 18 Sep 2017 - added explicit interface
414 ! HCW 08 Feb 2017 - fixed bug in Zh between grids, added default z0m, zdm
415 ! HCW 03 Mar 2015
416 ! sg feb 2012 - made separate subroutine
417 !--------------------------------------------------------------------------------
418
419 IMPLICIT NONE
420
421 INTEGER, PARAMETER :: nsurf = 7 ! number of surface types
422 INTEGER, PARAMETER :: PavSurf = 1 !When all surfaces considered together (1-7)
423 INTEGER, PARAMETER :: BldgSurf = 2
424 INTEGER, PARAMETER :: ConifSurf = 3
425 INTEGER, PARAMETER :: DecidSurf = 4
426 INTEGER, PARAMETER :: GrassSurf = 5 !New surface classes: Grass = 5th/7 surfaces
427 INTEGER, PARAMETER :: BSoilSurf = 6 !New surface classes: Bare soil = 6th/7 surfaces
428 INTEGER, PARAMETER :: WaterSurf = 7
429 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
430
431 INTEGER, INTENT(in) :: RoughLenMomMethod
432
433 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! surface fractions
434
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 ! z0m set in SiteSelect
443 REAL(KIND(1D0)), INTENT(in) :: zdm_in ! zdm set in SiteSelect
444 REAL(KIND(1D0)), INTENT(in) :: Z
445
446 REAL(KIND(1D0)), INTENT(out) :: FAI
447 REAL(KIND(1D0)), INTENT(out) :: PAI
448 REAL(KIND(1D0)), INTENT(out) :: Zh ! effective height of bluff bodies
449 REAL(KIND(1D0)), INTENT(out) :: z0m ! aerodynamic roughness length
450 REAL(KIND(1D0)), INTENT(out) :: zdm ! zero-plance displacement
451 REAL(KIND(1D0)), INTENT(out) :: ZZD ! z-zdm
452
453 INTEGER, PARAMETER :: notUsedI = -55
454 REAL(KIND(1D0)), PARAMETER :: notUsed = -55.5
455 REAL(KIND(1D0)) :: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water !Default values for roughness lengths [m]
456
457 !Total area of buildings and trees
458 ! areaZh = (sfr_surf(BldgSurf) + sfr_surf(ConifSurf) + sfr_surf(DecidSurf))
459 ! TS 19 Jun 2022: take porosity of trees into account; to be consistent with PAI calculation in RSL
460 pai = dot_product(sfr_surf([bldgsurf, conifsurf, decidsurf]), [1d0, 1 - porosity_evetr, 1 - porosity_dectr])
461
462 ! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
463 z0m4paved = 0.003 !estimate
464 z0m4grass = 0.02
465 z0m4bsoil = 0.002
466 z0m4water = 0.0005
467
468 !------------------------------------------------------------------------------
469 !If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
470 IF (pai /= 0) THEN
471 zh = dot_product( &
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)], &
475 [0d0, 0d0, 0d0], &
476 sfr_surf([bldgsurf, conifsurf, decidsurf]) > 0))
477
478 ! `1e-5` set to avoid numerical difficulty
479 fai = max(fai, 1e-5)
480
481 ELSE
482 zh = 0 !Set Zh to zero if areaZh = 0
483 fai = 1e-5
484 END IF
485
486 IF (zh /= 0) THEN
487 !Calculate z0m and zdm depending on the Z0 method
488 IF (roughlenmommethod == 2) THEN !Rule of thumb (G&O 1999)
489 z0m = 0.1*zh
490 zdm = 0.7*zh
491 ELSEIF (roughlenmommethod == 3) THEN !MacDonald 1998
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 ! lambdaP dependent as in Fig.1a of G&O (1999)
495 ! these are derived using digitalised points
496 zdm = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*pai) + 0.493*sigmoid(-5.17 + 32.7*pai))*zh
497 z0m = (0.00208 + &
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
501 END IF
502 ELSEIF (zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
503 IF (pai /= 0) CALL errorhint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, pai, notusedi)
504 !Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
505 IF (pai /= 1) THEN
506 z0m = (z0m4paved*sfr_surf(pavsurf) &
507 + z0m4grass*sfr_surf(grasssurf) &
508 + z0m4bsoil*sfr_surf(bsoilsurf) &
509 + z0m4water*sfr_surf(watersurf))/(1 - pai)
510 zdm = 0
511 CALL errorhint(15, 'Setting z0m and zdm using default values', z0m, zdm, notusedi)
512 ELSEIF (pai == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
513 z0m = 1
514 zdm = 7
515 CALL errorhint(15, 'Assuming mean height = 10 m, Setting z0m and zdm to default value', z0m, zdm, notusedi)
516 END IF
517 END IF
518
519 IF (roughlenmommethod == 1) THEN !use values set in SiteSelect
520 z0m = z0m_in
521 zdm = zdm_in
522 END IF
523
524 zzd = z - zdm
525
526 ! Error messages if aerodynamic parameters negative
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)

References errorhint(), and sigmoid().

Referenced by initialstate(), and suews_driver::suews_cal_main().

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)  G1,
real(kind(1d0)), intent(in)  G2,
real(kind(1d0)), intent(in)  G3,
real(kind(1d0)), intent(in)  G4,
real(kind(1d0)), intent(in)  G5,
real(kind(1d0)), intent(in)  G6,
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)  gfunc,
real(kind(1d0)), intent(out)  gsc,
real(kind(1d0)), intent(out)  ResistSurf 
)

Definition at line 111 of file suews_phys_resist.f95.

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

References errorhint().

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

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