SUEWS API Site
Documentation of SUEWS source code
suews_phys_rslprof.f95
Go to the documentation of this file.
1 module rsl_module
3  USE meteo, ONLY: rh2qa, qa2rh
4  USE allocatearray, ONLY: &
6  implicit none
7 
8  INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
9 
10 contains
11 
12  SUBROUTINE rslprofile( &
13  Zh, z0m, zdm, &
14  L_MOD, sfr, FAI, StabilityMethod, &
15  avcp, lv_J_kg, avdens, &
16  avU1, Temp_C, avRH, Press_hPa, zMeas, qh, qe, & ! input
17  T2_C, q2_gkg, U10_ms, RH2, &!output
18  dataoutLineRSL) ! output
19  !-----------------------------------------------------
20  ! calculates windprofiles using MOST with a RSL-correction
21  ! based on Harman & Finnigan 2007
22  !
23  ! last modified by:
24  ! NT 16 Mar 2019: initial version
25  ! TS 16 Oct 2019: improved consistency in parameters/varaibles within SUEWS
26  ! TODO how to improve the speed of this code
27  !
28  !-----------------------------------------------------
29 
30  IMPLICIT NONE
31 
32  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::sfr! surface fractions [-]
33  REAL(KIND(1d0)), INTENT(in):: zMeas ! height of atmospheric forcing [m]
34  REAL(KIND(1d0)), INTENT(in):: avU1 ! Wind speed at forcing height [m s-1]
35  REAL(KIND(1d0)), INTENT(in):: Temp_C ! Air temperature at forcing height [C]
36  REAL(KIND(1d0)), INTENT(in):: avRH ! relative humidity at forcing height [-]
37  REAL(KIND(1d0)), INTENT(in):: Press_hPa ! pressure at forcing height [hPa]
38  REAL(KIND(1d0)), INTENT(in):: L_MOD ! Obukhov length [m]
39  REAL(KIND(1d0)), INTENT(in):: avcp ! specific heat capacity [J kg-1 K-1]
40  REAL(KIND(1d0)), INTENT(in):: lv_J_kg ! Latent heat of vaporization in [J kg-1]
41  REAL(KIND(1d0)), INTENT(in):: avdens ! air density [kg m-3]
42  REAL(KIND(1d0)), INTENT(in):: qh ! sensible heat flux [W m-2]
43  REAL(KIND(1d0)), INTENT(in):: qe ! Latent heat flux [W m-2]
44  REAL(KIND(1d0)), INTENT(in):: Zh ! Mean building height [m]
45  REAL(KIND(1d0)), INTENT(in):: z0m ! Mean building height [m]
46  REAL(KIND(1d0)), INTENT(in):: zdm ! Mean building height [m]
47  REAL(KIND(1d0)), INTENT(in):: FAI ! Frontal area index [-]
48 
49  INTEGER, INTENT(in)::StabilityMethod
50 
51  REAL(KIND(1d0)), INTENT(out):: T2_C ! Air temperature at 2 m [C]
52  REAL(KIND(1d0)), INTENT(out):: q2_gkg ! Air specific humidity at 2 m [g kg-1]
53  REAL(KIND(1d0)), INTENT(out):: U10_ms ! wind speed at 10 m [m s-1]
54  REAL(KIND(1d0)), INTENT(out):: RH2 ! Air relative humidity [-]
55 
56  INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
57 
58  REAL(KIND(1d0)), PARAMETER:: cd_tree = 1.2, & ! drag coefficient tree canopy !!!!needs adjusting!!!
59  a_tree = 0.05, & ! the foliage area per unit volume !!!!needs adjusting!!!
60  kappa = 0.40, &! von karman constant
61  ! lv_J_kg = 2.5E6, &! latent heat for water vapor!!! make consistant with rest of code
62  beta_n = 0.40, & ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
63  pi = 4.*atan(1.0), r = 0.1, &
64  a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
65 
66  ! REAL(KIND(1d0)), INTENT(out), DIMENSION(nz*3):: zarrays ! Height array
67  REAL(KIND(1d0)), INTENT(out), DIMENSION(nz*4):: dataoutLineRSL ! Variables array (/U,T,q/)
68  REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z
69  REAL(KIND(1d0)), DIMENSION(nz):: psihath_z
70  REAL(KIND(1d0)), DIMENSION(nz):: dif
71  ! REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z, psihath_z
72  REAL(KIND(1d0)), DIMENSION(nz):: zarray
73  REAL(KIND(1d0)), DIMENSION(nz):: dataoutLineURSL ! wind speed array [m s-1]
74  REAL(KIND(1d0)), DIMENSION(nz):: dataoutLineTRSL ! Temperature array [C]
75  REAL(KIND(1d0)), DIMENSION(nz):: dataoutLineqRSL ! Specific humidity array [g kg-1]
76 
77  REAL(KIND(1d0))::z0 ! roughness length from H&F
78  REAL(KIND(1d0))::zd ! zero-plane displacement
79 
80  REAL(KIND(1d0))::Lc_build, Lc_tree, Lc ! canopy drag length scale
81  REAL(KIND(1d0))::Scc ! Schmidt number for temperature and humidity
82  REAL(KIND(1d0))::phim, psimz, psimZh, psimz0, psimza, phi_hatmZh, phi_hathZh, phimzp, phimz, phihzp, phihz, psihz, psihza ! stability function for momentum
83  REAL(KIND(1d0))::betaHF, betaNL, beta, betaN2 ! beta coefficient from Harman 2012
84  REAL(KIND(1d0))::elm ! mixing length
85  REAL(KIND(1d0))::xxm1, xxm1_2, xxh1, xxh1_2, err, z01, dphi, dphih ! dummy variables for stability functions
86  REAL(KIND(1d0))::f, cm, c2, ch, c2h ! H&F'07 and H&F'08 'constants'
87  REAL(KIND(1d0))::t_h, q_h ! H&F'08 canopy corrections
88  REAL(KIND(1d0))::TStar_RSL ! temperature scale
89  REAL(KIND(1d0))::UStar_RSL ! friction velocity used in RSL
90  REAL(KIND(1d0))::sfr_zh ! land cover fraction of bluff bodies: buildings and trees
91  REAL(KIND(1d0))::sfr_tr ! land cover fraction of trees
92  REAL(KIND(1d0))::L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
93  REAL(KIND(1d0))::zH_RSL ! mean canyon height used in RSL module with thresholds applied
94  REAL(KIND(1d0))::dz! initial height step
95  REAL(KIND(1d0)), parameter::Zh_min = 0.15! limit for minimum canyon height used in RSL module
96  REAL(KIND(1d0)), parameter::ratio_dz = 1.618! ratio between neighbouring height steps
97 
98  REAL(KIND(1d0))::qa_gkg, qStar_RSL ! specific humidity scale
99  INTEGER :: I, z, it, idx_can, idx_za, idx_2m, idx_10m
100  INTEGER :: nz_can ! number of heights in canyon
101  !
102  ! Step 1: Calculate grid-cell dependent constants
103  ! Step 2: Calculate Beta (crucial for H&F method)
104  ! Step 3: calculate the stability dependent H&F constants
105  ! Step 4: determine psihat at levels above the canopy
106  ! Step 5: Calculate z0 iteratively
107  ! Step 6: Calculate mean variables above canopy
108  ! Step 7: Calculate mean variables in canopy
109  !
110  ! ! Step 1
111  ! ! Start setting up the parameters
112 
113  call rsl_cal_prms( &
114  zh_min, &
115  z0m, zdm, &
116  stabilitymethod, zh, l_mod, sfr, fai, &!input
117  l_mod_rsl, zh_rsl, lc, beta, zd, z0, elm, scc, f)
118 
119  ! Define the height array with consideration of key heights
120  ! set number of heights within canopy
121  IF (zh_rsl <= 2) THEN
122  nz_can = 5
123  ELSE IF (zh_rsl <= 10) THEN
124  nz_can = 10
125  else
126  nz_can = 15
127  ENDIF
128  ! fill up heights in canopy
129  dz = zh_rsl/nz_can
130  do i = 1, nz_can
131  zarray(i) = dz*i
132  end do
133  ! fill up heights above canopy
134  dz = (zmeas - zh_rsl)/(nz - nz_can)
135  do i = nz_can + 1, nz
136  zarray(i) = zh_rsl + (i - nz_can)*dz
137  end do
138 
139  ! add key heights (2m and 10m) to zarray
140  ! 2m:
141  DO z = 1, nz
142  dif(z) = abs(zarray(z) - 2)
143  ENDDO
144  idx_2m = minloc(dif, dim=1)
145  zarray(idx_2m) = 2
146  ! 10m:
147  DO z = 1, nz
148  dif(z) = abs(zarray(z) - 10)
149  ENDDO
150  idx_10m = minloc(dif, dim=1)
151  zarray(idx_10m) = 10
152 
153  ! determine index at the canyon top
154  DO z = 1, nz
155  dif(z) = abs(zarray(z) - zh_rsl)
156  ENDDO
157  idx_can = minloc(dif, dim=1)
158  zarray(idx_can) = zh_rsl
159 
160  ! determine index at measurement height
161  DO z = 1, nz
162  ! dif2(z) = ABS(zarray(z) - (zMeas - zd))
163  dif(z) = abs(zarray(z) - zmeas)
164  ENDDO
165  idx_za = minloc(dif, dim=1)
166  zarray(idx_za) = zmeas
167 
168  if (zh_rsl - zd < z0 .or. zh < zh_min) then
169  ! correct parameters if RSL approach doesn't apply for a shallow canyon
170  ! when zh_RSL is too shallow, implying RSL doesn't apply, force RSL correction to zero
171  psihatm_z = 0
172  psihath_z = 0
173  beta = 1.e6
174  !correct RSL-based z0 and zd using Rule of thumb (G&O 1999)
175  ! zd = 0.7*zH_RSL
176  ! z0 = 0.1*zH_RSL
177  !correct RSL-based using SUEWS system-wide values
178  zd = zdm
179  z0 = z0m
180  ! then MOST recovers from RSL correction
181  else
182  !otherwise use RSL approach to calculate correction factors
183  ! Step 3:
184  !
185  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
186  xxm1 = stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
187  xxm1_2 = stab_phi_mom(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
188 
189  xxh1 = stab_phi_heat(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
190  xxh1_2 = stab_phi_heat(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
191 
192  phi_hatmzh = kappa/(2.*beta*xxm1)
193  phi_hathzh = kappa*scc/(2.*beta*xxh1)
194  dphi = xxm1_2 - xxm1
195  dphih = xxh1_2 - xxh1
196  IF (phi_hatmzh > 1.) THEN
197  c2 = 0.5 ! more stable, but less correct
198  c2h = 0.5
199  ELSE
200  c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa) ! if very unstable this might cause some high values of psihat_z
201  c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
202  ENDIF
203  cm = (1.-phi_hatmzh)*exp(c2/2.)
204  ch = (1.-phi_hathzh)*exp(c2h/2.)
205 
206  !
207  ! Step 4: determine psihat at levels above the canopy
208  !
209  psihatm_z = 0.*zarray
210  psihath_z = 0.*zarray
211  ! psihatm_zp = 0.
212  ! psihath_zp = 0
213  DO z = nz - 1, idx_can - 1, -1
214  phimz = stab_phi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
215  phimzp = stab_phi_mom(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
216  phihz = stab_phi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
217  phihzp = stab_phi_heat(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
218 
219  psihatm_z(z) = psihatm_z(z + 1) + dz/2.*phimzp*(cm*exp(-1.*c2*beta*(zarray(z + 1) - zd)/elm)) & !Taylor's approximation for integral
220  /(zarray(z + 1) - zd)
221  psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
222  /(zarray(z) - zd)
223  psihath_z(z) = psihath_z(z + 1) + dz/2.*phihzp*(ch*exp(-1.*c2h*beta*(zarray(z + 1) - zd)/elm)) & !Taylor's approximation for integral
224  /(zarray(z + 1) - zd)
225  psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
226  /(zarray(z) - zd)
227  ENDDO
228  ! psihatm_z=cal_psihatm_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm)
229  ! psihath_z=cal_psihath_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f)
230 
231  end if
232 
233  ! calculate above canopy variables
234  !
235  psimz0 = stab_psi_mom(stabilitymethod, z0/l_mod_rsl)
236  psimza = stab_psi_mom(stabilitymethod, (zmeas - zd)/l_mod_rsl)
237  psihza = stab_psi_heat(stabilitymethod, (zmeas - zd)/l_mod_rsl)
238  ustar_rsl = avu1*kappa/(log((zmeas - zd)/z0) - psimza + psimz0 + psihatm_z(nz))
239  tstar_rsl = -1.*(qh/(avcp*avdens))/ustar_rsl
240  qstar_rsl = -1.*(qe/lv_j_kg*avdens)/ustar_rsl
241  qa_gkg = rh2qa(avrh/100, press_hpa, temp_c)
242 
243  DO z = idx_can, nz
244  psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
245  psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
246  ! dataoutLineURSL(z) = (LOG((zarray(z) - zd)/z0) - psimz + psimz0 - psihatm_z(z) + psihatm_z(idx_can))/kappa ! this is different from Theeuwes et al. (2019 BLM)
247  dataoutlineursl(z) = (log((zarray(z) - zd)/z0) - psimz + psimz0 + psihatm_z(z))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
248  ! dataoutLineURSL(z) = (LOG((zarray(z) - zd)/(zMeas - zd)) - psimz + psimza + psihatm_z(z)-psihatm_z(idx_za))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
249  dataoutlinetrsl(z) = (log((zarray(z) - zd)/(zmeas - zd)) - psihz + psihza + psihath_z(z) - psihath_z(idx_za))/kappa ! eqn. 4 in Theeuwes et al. (2019 BLM)
250  dataoutlineqrsl(z) = (log((zarray(z) - zd)/(zmeas - zd)) - psihz + psihza + psihath_z(z) - psihath_z(idx_za))/kappa
251  ENDDO
252  !
253  ! Step 7
254  ! calculate in canopy variables
255  !
256  t_h = scc*tstar_rsl/(beta*f)
257  q_h = scc*qstar_rsl/(beta*f)
258  DO z = 1, idx_can
259  dataoutlineursl(z) = dataoutlineursl(idx_can)*exp(beta*(zarray(z) - zh_rsl)/elm)
260  dataoutlinetrsl(z) = dataoutlinetrsl(idx_can) + (t_h*exp(beta*f*(zarray(z) - zh_rsl)/elm) - t_h)/tstar_rsl
261  dataoutlineqrsl(z) = dataoutlineqrsl(idx_can) + (q_h*exp(beta*f*(zarray(z) - zh_rsl)/elm) - q_h)/qstar_rsl
262  ENDDO
263 
264  dataoutlineursl = dataoutlineursl*ustar_rsl
265  dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
266  dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
267 
268  dataoutlinersl = (/zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl/)
269 
270  !
271  ! Step 8
272  ! retrieve the diagnostics at key heights
273  !
274  t2_c = dataoutlinetrsl(idx_2m)
275  q2_gkg = dataoutlineqrsl(idx_2m)
276  u10_ms = dataoutlineursl(idx_10m)
277  ! get relative humidity:
278  rh2 = qa2rh(q2_gkg, press_hpa, t2_c)
279 
280 ! print *, 'Wind speed', dataoutLineURSL
281  ! DO z = 1, nz
282  ! print *, dataoutLineTRSL(z)
283  ! ENDDO
284  ! DO z = 1, nz
285  ! print *, zarray(z)
286  ! ENDDO
287 ! print *, 'Temperature', Temp_C, dataoutLineTRSL
288 ! print *, 'qStar', qStar, qe
289 ! print *, 'humidity' , qa_gkg, dataoutLineqRSL*1000.
290  END SUBROUTINE rslprofile
291 
292  recursive function cal_psim_hat(StabilityMethod, z, zh_RSL, L_MOD_RSL, beta, Lc) result(psim_hat_z)
293  ! calculate psi_hat for momentum
294  ! TS, 23 Oct 2019
295  implicit none
296  integer, intent(in) :: StabilityMethod ! stability method
297  real(KIND(1D0)), intent(in) :: z ! height of interest [m]
298  real(KIND(1D0)), intent(in) :: zh_RSL ! canyon depth [m]
299  real(KIND(1D0)), intent(in) :: Lc ! height scale for bluff bodies [m]
300  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
301  real(KIND(1D0)), intent(in) :: L_MOD_RSL ! Obukhov length [m]
302 
303  ! output
304  real(KIND(1D0)) ::psim_hat_z ! psim_hat at height of interest
305 
306  ! internal variables
307  real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
308  real(KIND(1D0)) ::zd ! displacement height used in RSL
309  real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
310  real(KIND(1D0)) ::phim_z ! displacement height used in RSL
311  real(KIND(1D0)) ::phim_zp ! displacement height used in RSL
312  real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
313  real(KIND(1D0)) ::elm ! displacement height used in RSL
314  real(KIND(1D0)) ::xxm1 ! displacement height used in RSL
315  real(KIND(1D0)) ::xxm1_2 ! displacement height used in RSL
316  real(KIND(1D0)) ::dphi ! displacement height used in RSL
317  real(KIND(1D0)) ::phi_hatmZh ! displacement height used in RSL
318  real(KIND(1D0)) ::cm ! displacement height used in RSL
319  real(KIND(1D0)) ::c2 ! displacement height used in RSL
320 
321  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
322  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
323 
324  if (z > 100) then
325  psim_hat_z = 0.
326  return
327  end if
328 
329  zp = 1.01*z ! a height above z
330 
331  zd = zh_rsl - (beta**2.)*lc
332  elm = 2.*beta**3*lc
333 
334  ! phim at Lc
335  phim_lc = stab_phi_mom(stabilitymethod, lc/l_mod_rsl)
336 
337  phim_z = stab_phi_mom(stabilitymethod, (z - zd)/l_mod_rsl)
338  phim_zp = stab_phi_mom(stabilitymethod, (zp - zd)/l_mod_rsl)
339 
340  psim_hat_zp = cal_psim_hat(stabilitymethod, zp, zh_rsl, l_mod_rsl, beta, lc)
341 
342  xxm1 = stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
343  xxm1_2 = stab_phi_mom(stabilitymethod, (zh_rsl - zd + dz)/l_mod_rsl)
344 
345  dphi = (xxm1_2 - xxm1)/dz
346 
347  phi_hatmzh = kappa/(2.*beta*xxm1)
348 
349  IF (phi_hatmzh > 1.) THEN
350  ! more stable, but less correct
351  c2 = 0.5
352  ELSE
353  ! if very unstable this might cause some high values of psihat_z
354  c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa)
355  ENDIF
356 
357  cm = (1.-phi_hatmzh)*exp(c2/2.)
358 
359  !Taylor's approximation for integral
360  psim_hat_z = psim_hat_zp + dz/2.*phim_zp*(cm*exp(-1.*c2*beta*(zp - zd)/elm))/(zp - zd)
361  psim_hat_z = psim_hat_z + dz/2.*phim_z*(cm*exp(-1.*c2*beta*(z - zd)/elm))/(z - zd)
362 
363  end function cal_psim_hat
364 
365  function cal_psihatm_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm) result(psihatm_z)
367  ! calculate psi_hat for momentum
368  ! TS, 23 Oct 2019
369  implicit none
370  integer, intent(in) :: StabilityMethod ! stability method
371  integer, intent(in) :: nz ! number of vertical layers
372  real(KIND(1D0)), DIMENSION(nz), intent(in) :: zarray ! height of interest [m]
373  real(KIND(1D0)), intent(in) :: zh_RSL ! canyon depth [m]
374  real(KIND(1D0)), intent(in) :: Lc ! height scale for bluff bodies [m]
375  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
376  real(KIND(1D0)), intent(in) :: L_MOD_RSL ! Obukhov length [m]
377  real(KIND(1D0)), intent(in) ::zd ! displacement height used in RSL
378  real(KIND(1D0)), intent(in) ::elm ! displacement height used in RSL
379 
380  ! output
381  real(KIND(1D0)), DIMENSION(nz) ::psihatm_z ! psim_hat at height of interest
382 
383  ! internal variables
384  ! real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
385  ! real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
386  real(KIND(1D0)) ::phimz ! displacement height used in RSL
387  real(KIND(1D0)) ::phimzp ! displacement height used in RSL
388  ! real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
389  real(KIND(1D0)) ::xxm1 ! displacement height used in RSL
390  real(KIND(1D0)) ::xxm1_2 ! displacement height used in RSL
391  real(KIND(1D0)) ::dphi ! displacement height used in RSL
392  real(KIND(1D0)) ::phi_hatmZh ! displacement height used in RSL
393  real(KIND(1D0)) ::cm ! displacement height used in RSL
394  real(KIND(1D0)) ::c2 ! displacement height used in RSL
395  REAL(KIND(1d0)), DIMENSION(nz):: dif
396 
397  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
398  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
399 
400  integer::z, idx_can
401 
402  psihatm_z = 0.*zarray
403 
404  ! determine index at the canyon top
405  DO z = 1, nz
406  dif(z) = abs(zarray(z) - zh_rsl)
407  ENDDO
408  idx_can = minloc(dif, dim=1)
409  ! zarray(idx_can) = Zh_RSL
410 
411  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
412  xxm1 = stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
413  xxm1_2 = stab_phi_mom(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
414 
415  phi_hatmzh = kappa/(2.*beta*xxm1)
416  dphi = xxm1_2 - xxm1
417 
418  IF (phi_hatmzh > 1.) THEN
419  c2 = 0.5 ! more stable, but less correct
420  ! c2h = 0.5
421  ELSE
422  c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa) ! if very unstable this might cause some high values of psihat_z
423  ! c2h = (kappa*Scc*(2.+f - (dphih*2.*beta**2.*Lc/xxh1)))/(2.*beta*xxh1 - kappa*Scc)
424  ENDIF
425  cm = (1.-phi_hatmzh)*exp(c2/2.)
426  ! ch = (1.-phi_hathzh)*EXP(c2h/2.)
427 
428  DO z = nz - 1, idx_can - 1, -1
429  phimz = stab_phi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
430  phimzp = stab_phi_mom(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
431 
432  psihatm_z(z) = psihatm_z(z + 1) + dz/2.*phimzp*(cm*exp(-1.*c2*beta*(zarray(z + 1) - zd)/elm)) & !Taylor's approximation for integral
433  /(zarray(z + 1) - zd)
434  psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
435  /(zarray(z) - zd)
436 
437  ENDDO
438 
439  end function cal_psihatm_z
440 
441  function cal_psihath_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f) result(psihath_z)
443  ! calculate psi_hat for momentum
444  ! TS, 23 Oct 2019
445  implicit none
446  integer, intent(in) :: StabilityMethod ! stability method
447  integer, intent(in) :: nz ! number of vertical layers
448 
449  real(KIND(1D0)), DIMENSION(nz), intent(in) :: zarray ! height of interest [m]
450  real(KIND(1D0)), intent(in) :: zh_RSL ! canyon depth [m]
451  real(KIND(1D0)), intent(in) :: Lc ! height scale for bluff bodies [m]
452  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
453  real(KIND(1D0)), intent(in) :: Scc ! parameter in RSL
454  real(KIND(1D0)), intent(in) :: f ! parameter in RSL
455  real(KIND(1D0)), intent(in) :: L_MOD_RSL ! Obukhov length [m]
456  real(KIND(1D0)), intent(in) :: elm ! displacement height used in RSL
457  real(KIND(1D0)), intent(in) ::zd ! displacement height used in RSL
458 
459  ! output
460  real(KIND(1D0)), DIMENSION(nz) ::psihath_z ! psim_hat at height of interest
461 
462  ! internal variables
463  ! real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
464  ! real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
465  real(KIND(1D0)) ::phihz ! displacement height used in RSL
466  real(KIND(1D0)) ::phihzp ! displacement height used in RSL
467  ! real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
468  real(KIND(1D0)) ::xxh1 ! displacement height used in RSL
469  real(KIND(1D0)) ::xxh1_2 ! displacement height used in RSL
470  real(KIND(1D0)) ::dphih ! displacement height used in RSL
471  real(KIND(1D0)) ::phi_hathZh ! displacement height used in RSL
472  real(KIND(1D0)) ::ch ! displacement height used in RSL
473  real(KIND(1D0)) ::c2h ! displacement height used in RSL
474  REAL(KIND(1d0)), DIMENSION(nz):: dif
475 
476  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
477  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
478 
479  integer::z, idx_can
480 
481  psihath_z = 0.*zarray
482 
483  ! determine index at the canyon top
484  DO z = 1, nz
485  dif(z) = abs(zarray(z) - zh_rsl)
486  ENDDO
487  idx_can = minloc(dif, dim=1)
488  ! zarray(idx_can) = Zh_RSL
489 
490  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
491  xxh1 = stab_phi_heat(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
492  xxh1_2 = stab_phi_heat(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
493 
494  phi_hathzh = kappa*scc/(2.*beta*xxh1)
495  dphih = xxh1_2 - xxh1
496 
497  IF (phi_hathzh > 1.) THEN
498  ! c2 = 0.5 ! more stable, but less correct
499  c2h = 0.5
500  ELSE
501  ! c2 = (kappa*(3.-(2.*beta**2.*Lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa) ! if very unstable this might cause some high values of psihat_z
502  c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
503  ENDIF
504  ! cm = (1.-phi_hatmZh)*EXP(c2/2.)
505  ch = (1.-phi_hathzh)*exp(c2h/2.)
506 
507  DO z = nz - 1, idx_can - 1, -1
508  phihz = stab_phi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
509  phihzp = stab_phi_heat(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
510 
511  psihath_z(z) = psihath_z(z + 1) + dz/2.*phihzp*(ch*exp(-1.*c2h*beta*(zarray(z + 1) - zd)/elm)) & !Taylor's approximation for integral
512  /(zarray(z + 1) - zd)
513  psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
514  /(zarray(z) - zd)
515 
516  ENDDO
517 
518  end function cal_psihath_z
519 
520  function rsl_cal_z0(StabilityMethod, zH_RSL, zd, beta, L_MOD_RSL, Lc) result(z0)
521  ! calculate z0 iteratively
522  ! TS, 23 Oct 2019
523  implicit none
524  integer, intent(in) ::StabilityMethod
525  real(KIND(1D0)), intent(in) :: zH_RSL ! canyon depth [m]
526  real(KIND(1D0)), intent(in) :: zd ! displacement height [m]
527  real(KIND(1D0)), intent(in) :: L_MOD_RSL ! Monin Obukhov length[m]
528  real(KIND(1D0)), intent(in) :: Lc ! canyon length scale [m]
529  real(KIND(1D0)), intent(in) :: beta ! height scale for bluff bodies [m]
530 
531  ! output
532  real(KIND(1D0)) ::z0
533 
534  ! internal variables
535  real(KIND(1D0)) ::psimZh, psimz0, z01, psihatm_Zh
536  real(KIND(1D0)) ::err
537  integer ::it
538 
539  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
540  ! REAL(KIND(1d0)), PARAMETER::r = 0.1
541  ! REAL(KIND(1d0)), PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
542 
543  psimzh = stab_psi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
544  psihatm_zh = cal_psim_hat(stabilitymethod, zh_rsl, zh_rsl, l_mod_rsl, beta, lc)
545 
546  !first guess
547  z0 = 0.5
548  ! if (Zh_RSL>Zh_min) then
549  ! if Zh>Zh_min, calculate z0 using RSL method
550  err = 10.
551  ! psimz0 = 0.5
552  it = 1
553  DO WHILE ((err > 0.001) .AND. (it < 10))
554  psimz0 = stab_psi_mom(stabilitymethod, z0/l_mod_rsl)
555  z01 = z0
556  z0 = (zh_rsl - zd)*exp(-1.*kappa/beta)*exp(-1.*psimzh + psimz0)*exp(psihatm_zh)
557  err = abs(z01 - z0)
558  it = it + 1
559  ENDDO
560 
561  end function rsl_cal_z0
562 
563  subroutine rsl_cal_prms( &
564  zh_min, z0m, zdm, &
565  StabilityMethod, zh, L_MOD, sfr, planF, &!input
566  L_MOD_RSL, zH_RSL, Lc, beta, zd, z0, elm, Scc, f)!output
567  ! calculate surface/skin temperature
568  ! TS, 23 Oct 2019
569  implicit none
570  integer, intent(in) :: StabilityMethod ! stability method
571  real(KIND(1D0)), intent(in) :: zh_min ! canyon depth [m]
572  real(KIND(1D0)), intent(in) :: z0m ! roughness length to prescribe if necessary [m]
573  real(KIND(1D0)), intent(in) :: zdm ! displacement height to prescribe if necessary [m]
574  real(KIND(1D0)), intent(in) :: zh ! canyon depth [m]
575  real(KIND(1D0)), intent(in) :: planF ! frontal area index
576  real(KIND(1D0)), intent(in) :: L_MOD ! Obukhov length [m]
577  real(KIND(1D0)), DIMENSION(nsurf), intent(in) :: sfr ! land cover fractions
578 
579  ! output
580  real(KIND(1D0)), intent(out) ::L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
581  real(KIND(1D0)), intent(out) ::zH_RSL ! mean canyon height used in RSL module with thresholds applied
582  real(KIND(1D0)), intent(out) ::Lc ! height scale for bluff bodies [m]
583  real(KIND(1D0)), intent(out) ::beta ! psim_hat at height of interest
584  real(KIND(1D0)), intent(out) ::zd ! displacement height to prescribe if necessary [m]
585  real(KIND(1D0)), intent(out) ::z0 ! roughness length [m]
586  real(KIND(1D0)), intent(out) ::elm ! length scale used in RSL
587  real(KIND(1D0)), intent(out) ::Scc ! parameter in RSL
588  real(KIND(1D0)), intent(out) ::f ! parameter in RSL
589 
590  ! internal variables
591  real(KIND(1D0)) ::sfr_zh
592  real(KIND(1D0)) ::sfr_tr
593  real(KIND(1D0)) ::phim
594  real(KIND(1D0)) ::betaN2
595  real(KIND(1D0)) ::betaHF
596  real(KIND(1D0)) ::betaNL
597 
598  REAL(KIND(1d0)), PARAMETER::planF_low = 1e-6
599  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
600  ! REAL(KIND(1d0)), PARAMETER::z0m= 0.40
601  REAL(KIND(1d0)), PARAMETER::r = 0.1
602  REAL(KIND(1d0)), PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
603  ! REAL(KIND(1d0)), parameter::Zh_min = 0.4! limit for minimum canyon height used in RSL module
604 
605  ! under stable conditions, set a threshold for L_MOD to avoid numerical issues. TS 28 Oct 2019
606  l_mod_rsl = merge(l_mod, max(300., l_mod), l_mod < 0)
607 
608  ! zH_RSL
609  zh_rsl = max(zh, zh_min)
610  ! zH_RSL = zh
611 
612  ! land cover fraction of bluff bodies
613  sfr_zh = sum(sfr([bldgsurf, conifsurf, decidsurf]))
614  ! set a threshold for sfr_zh to avoid numerical difficulties
615  sfr_zh = min(sfr_zh, 0.8)
616 
617  ! land cover fraction of trees
618  sfr_tr = sum(sfr([conifsurf, decidsurf]))
619 
620  ! height scale for buildings !not used? why?
621  ! Lc_build = (1.-sfr(BldgSurf))/planF*Zh_RSL ! Coceal and Belcher 2004 assuming Cd = 2
622 
623  ! height scale for tress
624  ! Lc_tree = 1./(cd_tree*a_tree) ! not used? why?
625 
626  ! height scale for bluff bodies
627  lc = (1.-sfr_zh)/planf*zh_rsl
628 
629  ! phim at Lc
630  phim = stab_phi_mom(stabilitymethod, lc/l_mod_rsl)
631 
632  scc = 0.5 + 0.3*tanh(2.*lc/l_mod_rsl) ! Schmidt number Harman and Finnigan 2008: assuming the same for heat and momemntum
633  f = 0.5*((1.+4.*r*scc)**0.5) - 0.5
634 
635  !
636  ! Step 2:
637  ! Parameterise beta according to Harman 2012 with upper limit of 0.5
638  ! betaN for trees found to be 0.3 and for urban 0.4 linearly interpolate between the two using surface fractions
639  ! betaN2 = 0.30 + (1.-sfr(ConifSurf) - sfr(ConifSurf))*0.1
640  if (sfr_zh > 0) then
641  betan2 = 0.30*sfr_tr/sfr_zh + (sfr_zh - sfr_tr)/sfr_zh*0.4
642  ELSE
643  betan2 = 0.35
644  endif
645 
646  betahf = betan2/phim
647  betanl = (kappa/2.)/phim
648 
649  IF (lc/l_mod_rsl > a2) THEN
650  beta = betahf
651  ELSE
652  beta = betanl + ((betahf - betanl)/(1.+a1*abs(lc/l_mod_rsl - a2)**a3))
653  ENDIF
654 
655  IF (beta > 0.5) THEN
656  beta = 0.5
657  ENDIF
658  zd = zh_rsl - (beta**2.)*lc
659  elm = 2.*beta**3*lc
660 
661  ! calculate z0 iteratively
662  z0 = rsl_cal_z0(stabilitymethod, zh_rsl, zd, beta, l_mod_rsl, lc)
663 
664  ! correct parameters if RSL approach doesn't apply for a shallow canyon
665  ! if (zh_RSL - zd < z0) then
666  ! ! when zh_RSL is too shallow, implying RSL doesn't apply, force RSL correction to zero
667  ! ! psihatm_z=0
668  ! ! psihath_z=0
669  ! beta = 1.e6
670  ! !correct RSL-based z0 and zd using Rule of thumb (G&O 1999)
671  ! zd = 0.7*zH_RSL
672  ! z0 = 0.1*zH_RSL
673  ! zd = zdm
674  ! z0 = z0m
675  ! ! then MOST recovers from RSL correction
676  ! end if
677 
678  end subroutine rsl_cal_prms
679 
680 end module rsl_module
real(kind(1d0)) function stab_phi_mom(StabilityMethod, ZL)
real(kind(1d0)) function, dimension(nz) cal_psihatm_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm)
real(kind(1d0)) function qa2rh(qa_gkg, pres_hPa, Ta_degC)
subroutine rsl_cal_prms(zh_min, z0m, zdm, StabilityMethod, zh, L_MOD, sfr, planF, L_MOD_RSL, zH_RSL, Lc, beta, zd, z0, elm, Scc, f)
real(kind(1d0)) function, dimension(nz) cal_psihath_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f)
integer, parameter nsurf
subroutine cal_stab(StabilityMethod, zzd, z0m, zdm, avU1, Temp_C, QH_init, avdens, avcp, L_MOD, TStar, UStar, zL)
integer, parameter conifsurf
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) function rsl_cal_z0(StabilityMethod, zH_RSL, zd, beta, L_MOD_RSL, Lc)
integer, parameter nz
recursive real(kind(1d0)) function cal_psim_hat(StabilityMethod, z, zh_RSL, L_MOD_RSL, beta, Lc)
real(kind(1d0)) function rh2qa(RH_dec, pres_hPa, Ta_degC)
real(kind(1d0)) function stab_phi_heat(StabilityMethod, ZL)
subroutine rslprofile(Zh, z0m, zdm, L_MOD, sfr, FAI, StabilityMethod, avcp, lv_J_kg, avdens, avU1, Temp_C, avRH, Press_hPa, zMeas, qh, qe, T2_C, q2_gkg, U10_ms, RH2, dataoutLineRSL)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
integer, parameter decidsurf
integer, parameter bldgsurf