SUEWS API Site
Documentation of SUEWS source code
suews_phys_evap.f95
Go to the documentation of this file.
1 module evap_module
2  implicit none
3 
4 contains
5  SUBROUTINE cal_evap( &
6  EvapMethod, state_is, WetThresh_is, capStore_is, &!input
7  vpd_hPa, avdens, avcp, qn_e, s_hPa, psyc_hPa, ResistSurf, RA, rb, tlv, &
8  rss, ev, qe) !output
9  !------------------------------------------------------------------------------
10  !-Calculates evaporation for each surface from modified Penman-Monteith eqn
11  !-State determines whether each surface type is dry or wet (wet/transition)
12  !-Wet surfaces below storage capacity are in transition
13  ! and QE depends on the state and storage capacity (i.e. varies with surface);
14  ! for wet or dry surfaces QE does not vary between surface types
15  !-See Sect 2.4 of Jarvi et al. (2011) Ja11
16  !
17  !Last modified:
18  ! HCW 06 Jul 2016
19  ! Moved rss declaration to LUMPS_Module_Constants so it can be written out
20  ! HCW 11 Jun 2015
21  ! Added WetThresh to distinguish wet/partially wet surfaces from the storage capacities used in SUEWS_drain
22  ! HCW 30 Jan 2015
23  ! Removed StorCap input because it is provided by module allocateArray
24  ! Tidied and commented code
25  ! LJ 10/2010
26  !------------------------------------------------------------------------------
27 
28  IMPLICIT NONE
29  INTEGER, INTENT(in) :: EvapMethod!Evaporation calculated according to Rutter (1) or Shuttleworth (2)
30 
31  REAL(KIND(1d0)), INTENT(in)::state_is ! wetness status
32  REAL(KIND(1d0)), INTENT(in)::WetThresh_is!When State > WetThresh, RS=0 limit in SUEWS_evap [mm] (specified in input files)
33  REAL(KIND(1d0)), INTENT(in)::capStore_is ! = StoreDrainPrm(6,is), current storage capacity [mm]
34 
35  REAL(KIND(1d0)), INTENT(in)::vpd_hPa ! vapour pressure deficit [hPa]
36  REAL(KIND(1d0)), INTENT(in)::avdens ! air density
37  REAL(KIND(1d0)), INTENT(in)::avcp ! air heat capacity
38  REAL(KIND(1d0)), INTENT(in)::qn_e !net available energy for evaporation
39  REAL(KIND(1d0)), INTENT(in)::s_hPa!Vapour pressure versus temperature slope in hPa
40  REAL(KIND(1d0)), INTENT(in)::psyc_hPa!Psychometric constant in hPa
41  REAL(KIND(1d0)), INTENT(in)::ResistSurf!Surface resistance
42  ! REAL(KIND(1d0)),INTENT(in)::sp!Term in calculation of E
43  REAL(KIND(1d0)), INTENT(in)::RA!Aerodynamic resistance
44  REAL(KIND(1d0)), INTENT(in)::rb!Boundary layer resistance
45  REAL(KIND(1d0)), INTENT(in)::tlv!Latent heat of vaporization per timestep [J kg-1 s-1], (tlv=lv_J_kg/tstep_real)
46 
47  REAL(KIND(1d0)), INTENT(out)::rss !Redefined surface resistance for wet
48  REAL(KIND(1d0)), INTENT(out)::ev ! evapotranspiration [mm]
49  REAL(KIND(1d0)), INTENT(out)::qe ! latent heat flux [W m-2]
50 
51  REAL(KIND(1d0))::numPM!numerator of P-M eqn
52  REAL(KIND(1d0))::rbsg !Boundary-layer resistance x (slope/psychrometric const + 1) [s m-1]
53  REAL(KIND(1d0))::rsrbsg !RS + rbsg [s m-1]
54  REAL(KIND(1d0))::rst
55  REAL(KIND(1d0))::W !Depends on the amount of water on the canopy [-]
56  REAL(KIND(1d0))::x
57  REAL(KIND(1d0))::r
58 
59  REAL(KIND(1d0)), PARAMETER:: NAN = -999
60 
61  ! Use Penman-Monteith eqn modified for urban areas (Eq6, Jarvi et al. 2011)
62  ! Calculation independent of surface characteristics
63  ! Uses value of RS for whole area (calculated based on LAI of veg surfaces in SUEWS_SurfaceResistance.f95)
64 
65  ! PRINT*, 'is',is,'SMOIS',state(is)
66  ! PRINT*, 'SMOIS',state_is,state_is<=0.001
67  ! PRINT*, 'EvapMethod',EvapMethod
68 
69  !numerator of P-M eqn, refer to Eq6, Jarvi et al. 2011
70  numpm = s_hpa*qn_e + vpd_hpa*avdens*avcp/ra !s_haPa - slope of svp vs t curve.
71 
72  ! Dry surface ---------------------------------------------------------------
73  IF (state_is <= 0.001) THEN
74  qe = numpm/(s_hpa + psyc_hpa*(1 + resistsurf/ra)) !QE [W m-2] (numPM = numerator of P-M eqn)
75  ev = qe/tlv !Ev [mm] (qe[W m-2]/tlv[J kg-1 s-1]*1/density_water[1000 kg m-3])
76  w = nan !W not needed for dry surfaces (set to -999)
77  rst = 1 !Set flag indicating dry surface(1)
78 
79  ! Wet surface ---------------------------------------------------------------
80  ELSE
81  rst = 0 !Set flag=0 indicating wet surface(0)
82 
83  ! Evaporation calculated according to Rutter(EvapMethod=1) or Shuttleworth(EvapMethod=2).
84  !Set in SUEWS_initial (so not an input to the model)
85  IF (evapmethod == 2) THEN !-- Shuttleworth (1978) --
86  rbsg = rb*(s_hpa/psyc_hpa + 1) !Boundary-layer resistance x (slope/psychro + 1)
87  rsrbsg = resistsurf + rbsg !RS + rsbg
88 
89  ! If surface is completely wet, set RS to zero -------------------
90  !if(state(is)>=StoreDrainPrm(6,is).or.ResistSurf<25) then !If at storage capacity or RS is small
91  IF (state_is >= wetthresh_is .OR. resistsurf < 25) THEN !If at storage capacity or RS is small
92  w = 1 !So that RS=0 (Eq7, Jarvi et al. 2011)
93  ! If surface is in transition, use rss ---------------------------
94  ELSE !if((state(is)<StorCap).and.(state(is)>0.001).or.(ResistSurf<50)) then
95  r = (resistsurf/ra)*(ra - rb)/rsrbsg
96  w = (r - 1)/(r - (wetthresh_is/state_is))
97  ENDIF
98 
99  ! PRINT*, 'r',r
100  ! PRINT*, 'W',W
101 
102  rss = (1/((w/rbsg) + ((1 - w)/rsrbsg))) - rbsg !Redefined surface resistance for wet
103  ! PRINT*, 'resistances:',rbsg,rsrbsg,rss
104  !surfaces (zero if W=1). Eq7, Jarvi et al. (2011)
105  qe = numpm/(s_hpa + psyc_hpa*(1 + rss/ra)) !QE [W m-2]
106  ev = qe/tlv !Ev [mm]
107  ! PRINT*, 'numPM',numPM
108  ! PRINT*, 'qe',qe
109 
110  ELSEIF (evapmethod == 1) THEN !-- Rutter --
111  qe = numpm/(s_hpa + psyc_hpa)
112  ev = qe/tlv
113 
114  x = merge(1d0, state_is/capstore_is, state_is > capstore_is)
115  ev = ev*x !QE [W m-2]
116  qe = ev*tlv !Ev [mm]
117  ENDIF !Rutter/Shuttleworth calculation
118  ENDIF !Wet/dry surface
119 
120  ! IF ( id>190 ) THEN
121  ! STOP "stop in Evap_SUEWS_new"
122  !
123  ! END IF
124  END SUBROUTINE cal_evap
125 
126 end module evap_module
subroutine cal_evap(EvapMethod, state_is, WetThresh_is, capStore_is, vpd_hPa, avdens, avcp, qn_e, s_hPa, psyc_hPa, ResistSurf, RA, rb, tlv, rss, ev, qe)