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)
 
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)
 
subroutine boundarylayerresistance (zzd, z0m, avU1, UStar, rb)
 
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 cal_z0v (RoughLenHeatMethod, z0m, VegFraction, UStar)
 

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 
)

Definition at line 16 of file suews_phys_resist.f95.

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

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
real(kind(1d0)) notused
real(kind(1d0)) k2
real(kind(1d0)) zzd
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) psih
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
real(kind(1e10)) z0v
real(kind(1d0)) z0m
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
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 378 of file suews_phys_resist.f95.

Referenced by suews_driver::suews_cal_resistance().

378 
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
real(kind(1d0)) k
real(kind(1d0)) zzd
real(kind(1d0)) z0m
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 518 of file suews_phys_resist.f95.

Referenced by aerodynamicresistance().

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 
real(kind(1e10)) z0v
real(kind(1d0)) z0m
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,
real(kind(1d0)), intent(in)  bldgH,
real(kind(1d0)), intent(in)  EveTreeH,
real(kind(1d0)), intent(in)  DecTreeH,
real(kind(1d0)), intent(in)  porosity_id,
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)  planF,
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 408 of file suews_phys_resist.f95.

References errorhint().

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

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)
real(kind(1d0)) z
real(kind(1d0)) notused
real(kind(1d0)) zdm
real(kind(1d0)) zdm_in
real(kind(1d0)) zzd
real(kind(1d0)) z0m_in
real(kind(1d0)) bldgh
real(kind(1d0)) z0m
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
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,
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 114 of file suews_phys_resist.f95.

References errorhint().

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

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
real(kind(1d0)) tc2
real(kind(1d0)) g1
real(kind(1d0)) s2
integer gsmodel
real(kind(1d0)) g5
real(kind(1d0)) kmax
real(kind(1d0)) th
real(kind(1d0)) notused
real(kind(1d0)) g2
real(kind(1d0)) tc
real(kind(1d0)) g3
integer id
real(kind(1d0)) dq
real(kind(1d0)) g4
real(kind(1d0)) g6
integer it
real(kind(1d0)) tl
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
real(kind(1d0)) s1
Here is the call graph for this function:
Here is the caller graph for this function: