SUEWS API Site
Documentation of SUEWS source code
suews_phys_atmmoiststab.f95
Go to the documentation of this file.
1MODULE atmmoiststab_module
2 IMPLICIT NONE
3 REAL(kind(1d0)), PARAMETER :: neut_limit = 1.e-2 !Limit for neutral stability
4 REAL(kind(1d0)), PARAMETER :: k = 0.4 !Von Karman's contant
5 REAL(kind(1d0)), PARAMETER :: grav = 9.80665 !g - gravity - physics today august 1987
6
7 ! scheme code
8 INTEGER, PARAMETER :: j12 = 2 ! Cheng and Brutsaert (2005)
9 INTEGER, PARAMETER :: k75 = 3 ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
10 INTEGER, PARAMETER :: b71 = 4 ! Businger et al (1971) modifed Hogstrom (1988)
11
12CONTAINS
13 !.c!! For Lumps Version 2 - no stability calculations
14 ! Latent heat of sublimation when air temperature below zero added. LJ Nov 2012
15 ! explict interface added to all subroutines, TS 08 Aug 2017
16 SUBROUTINE cal_atmmoist( &
17 Temp_C, Press_hPa, avRh, dectime, & ! input:
18 lv_J_kg, lvS_J_kg, & ! output:
19 es_hPa, Ea_hPa, VPd_hpa, VPD_Pa, dq, dens_dry, avcp, air_dens)
20
21 USE meteo, ONLY: &
22 sat_vap_press_x, spec_heat_beer, &
23 lat_vap, lat_vapsublim, spec_hum_def
24
25 IMPLICIT NONE
26 REAL(KIND(1D0)) :: vap_dens
27
28 REAL(KIND(1D0)), INTENT(in) :: &
29 Temp_C, &
30 Press_hPa, &
31 avRh, dectime
32 REAL(KIND(1D0)), INTENT(out) :: &
33 lv_J_kg, & !Latent heat of vaporization in [J kg-1]
34 lvS_J_kg, & !Latent heat of sublimation in J/kg
35 es_hPa, & !Saturation vapour pressure over water in hPa
36 Ea_hPa, & !Vapour pressure of water in hPa
37 VPd_hpa, & !vapour pressure deficit in hPa
38 VPD_Pa, & !vapour pressure deficit in Pa
39 dq, & !Specific humidity deficit in g/kg
40 dens_dry, & !Vap density or absolute humidity [kg m-3]
41 avcp, & !specific heat capacity in J kg-1 K-1
42 air_dens !Air density in kg/m3
43
44 REAL(KIND(1D0)), PARAMETER :: &
45 ! comp = 0.9995, &
46 ! epsil = 0.62197,& !ratio molecular weight of water vapor/dry air (kg/mol/kg/mol)
47 ! epsil_gkg = 621.97, & !ratio molecular weight of water vapor/dry air in g/kg
48 ! dry_gas = 8.31451,& !Dry gas constant (J/k/mol)
49 ! gas_ct_wat = 461.05,& !Gas constant for water (J/kg/K)
50 ! molar = 0.028965,& !Dry air molar fraction in kg/mol
51 ! molar_wat_vap = 0.0180153,& !Molar fraction of water vapor in kg/mol
52 gas_ct_dry = 8.31451/0.028965, & !j/kg/k=dry_gas/molar
53 gas_ct_wv = 8.31451/0.0180153 !j/kg/k=dry_gas/molar_wat_vap
54 ! waterDens = 999.8395 !Density of water in 0 cel deg
55 INTEGER :: from = 1
56
57 !Saturation vapour pressure over water in hPa
58 es_hpa = sat_vap_press_x(temp_c, press_hpa, from, dectime) ! dectime is more or less unnecessary here
59
60 !Vapour pressure of water in hPa
61 ea_hpa = avrh/100*es_hpa
62
63 ! if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'%',Temp_C
64
65 vpd_hpa = es_hpa - ea_hpa !vapour pressure deficit in hPa
66 vpd_pa = (es_hpa*100) - (ea_hpa*100) !vapour pressure deficit in Pa
67
68 dq = spec_hum_def(vpd_hpa, press_hpa) !Specific humidity deficit in g/kg
69
70 !Vap density or absolute humidity (kg/m3)
71 vap_dens = (ea_hpa*100/((temp_c + 273.16)*gas_ct_wv))
72
73 !density Dry Air Beer(1990) kg/m3
74 dens_dry = ((press_hpa - ea_hpa)*100)/(gas_ct_dry*(273.16 + temp_c))
75
76 !Air density in kg/m3
77 air_dens = (press_hpa*100)/(gas_ct_dry*(temp_c + 273.16))
78
79 !Calculate specific heat capacity in J kg-1 K-1
80 avcp = spec_heat_beer(temp_c, avrh, vap_dens, dens_dry)
81
82 !Latent heat of vaporization in [J kg-1]
83 lv_j_kg = lat_vap(temp_c, ea_hpa, press_hpa, avcp, dectime)
84
85 !Latent heat of sublimation in J/kg
86 IF (temp_c < 0.000) THEN
87 lvs_j_kg = lat_vapsublim(temp_c, ea_hpa, press_hpa, avcp)
88 ELSE
89 lvs_j_kg = 2834000
90 END IF
91 !if(debug)write(*,*)lv_J_kg,Temp_C,'lv2'
92 IF (press_hpa < 900) THEN
93 CALL errorhint(46, 'Function LUMPS_cal_AtmMoist', press_hpa, -55.55d0, -55)
94 END IF
95 RETURN
96 END SUBROUTINE cal_atmmoist
97
98 !.c!! For Lumps Version 2 - no stability calculations
99 !==========================================================
100 ! Last change:
101 ! TS 08 Aug 2017: added explicit interface
102 ! TS 13 Jun 2017: corrections to the integral of stability functions
103 ! MH 12 Apr 2017: Stable limit to exit do-loop
104 ! LJ 25 Nov 2014: Limits for L
105 ! LJ 19 Feb 2010
106 ! SG 27 Mar 2000 4:44 pm
107 ! ust - friction velocity
108 ! L - monin obukhov stability length
109 ! Van Ulden & Holtslag (1985) JCAM: 24: 1196-1207
110
111 SUBROUTINE cal_stab( &
112 StabilityMethod, & ! input
113 zzd, & !Active measurement height (meas. height-displac. height)
114 z0m, & !Aerodynamic roughness length
115 zdm, & !Displacement height
116 avu1, & !Average wind speed
117 temp_c, & !Air temperature
118 qh_init, & !sensible heat flux [W m-2]
119 avdens, & ! air density
120 avcp, & ! heat capacity of air
121 l_mod, & !Obukhov length! output:
122 tstar, & !T*
123 ustar, & !Friction velocity
124 zl) !Stability scale
125
126 IMPLICIT NONE
127 INTEGER, INTENT(in) :: StabilityMethod
128
129 REAL(KIND(1D0)), INTENT(in) :: zzd !Active measurement height (meas. height-displac. height)
130 REAL(KIND(1D0)), INTENT(in) :: z0m !Aerodynamic roughness length
131 REAL(KIND(1D0)), INTENT(in) :: zdm !Displacement height
132 REAL(KIND(1D0)), INTENT(in) :: avU1 !Average wind speed
133 REAL(KIND(1D0)), INTENT(in) :: Temp_C !Air temperature
134 REAL(KIND(1D0)), INTENT(in) :: QH_init !sensible heat flux [W m-2]
135 REAL(KIND(1D0)), INTENT(in) :: avdens ! air density [kg m-3]
136 REAL(KIND(1D0)), INTENT(in) :: avcp ! volumetric heat capacity [J m-3 K-1]
137
138 REAL(KIND(1D0)), INTENT(out) :: L_MOD !Obukhov length
139 REAL(KIND(1D0)), INTENT(out) :: TStar !T*
140 REAL(KIND(1D0)), INTENT(out) :: UStar !Friction velocity
141 REAL(KIND(1D0)), INTENT(out) :: zL ! Stability scale
142 ! REAL(KIND(1d0)),INTENT(out)::psim !Stability function of momentum
143
144 REAL(KIND(1D0)) :: G_T_K, &
145 KUZ, &
146 LOLD, &
147 psim, &
148 z0l, &
149 psimz0, &
150 H_init, &
151 h
152 INTEGER, PARAMETER :: notUsedI = -55
153
154 INTEGER :: i
155
156 LOGICAL :: debug = .false.
157
158 !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
159 h_init = qh_init/(avdens*avcp)
160
161 IF (debug) WRITE (*, *) stabilitymethod, z0m, avu1, h_init, ustar, l_mod
162 g_t_k = (grav/(temp_c + 273.16))*k !gravity constant/(Temperature*Von Karman Constant)
163 kuz = k*avu1 !Von Karman constant*mean wind speed
164 IF (zzd < 0) CALL errorhint(32, &
165 'Windspeed Ht too low relative to zdm [Stability calc]- values [z-zdm, zdm]', &
166 zzd, zdm, notusedi)
167
168 ustar = kuz/log(zzd/z0m) !Initial setting of u* and calc. of L_MOD (neutral situation)
169 IF (abs(h_init) < 0.001) THEN ! prevent zero TStar
170 h = 0.001
171 ELSE
172 h = h_init
173 END IF
174 tstar = (-h/ustar)
175 l_mod = (ustar**2)/(g_t_k*tstar)
176
177 IF (log(zzd/z0m) < 0.001000) THEN
178 ! PRINT*, 1/(z0m-z0m)
179 CALL errorhint(17, 'In stability subroutine, (z-zd) < z0.', zzd, z0m, notusedi)
180 END IF
181 i = 1
182 lold = -999.
183 z0l = z0m/l_mod !z0m roughness length
184 DO WHILE ((abs(lold - l_mod) > 0.01) .AND. (i < 330)) !NT: add error threshold !Iteration starts
185 lold = l_mod
186 zl = zzd/l_mod
187 z0l = z0m/l_mod !z0m roughness length
188
189 ! IF (zL>2)THEN
190 ! CALL ErrorHint(73,'LUMPS_atmos_functions_stab.f95: stability scale (z/L), UStar',zL,UStar,notUsedI)
191 ! RETURN !MO-theory not necessarily valid above zL>2. Still causes problematic UStar values and correct limit might be 0.3.
192 ! !Needs more investigations.
193 ! END IF
194
195 psim = stab_psi_mom(stabilitymethod, zl)
196 psimz0 = stab_psi_mom(stabilitymethod, z0l)
197
198 ustar = kuz/(log(zzd/z0m) - psim + psimz0) !Friction velocity in non-neutral situation
199
200 tstar = (-h/ustar)
201 l_mod = (ustar**2)/(g_t_k*tstar)
202
203 ! IF(UStar<0.001000)THEN !If u* too small
204 ! UStar=KUZ/(LOG(Zzd/z0m))
205 ! PRINT*, 'UStar info',UStar,KUZ,(LOG(Zzd/z0m)),Zzd,z0m
206 ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] zl,dectime',zl,dectime,notUsedI)
207 ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] z0l,UStar',z0l,UStar,notUsedI)
208 ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] psim,psimz0',psim,psimz0,notUsedI)
209 ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] AVU1,log(zzd/z0m)',AVU1,LOG(zzd/z0m),notUsedI)
210 !
211 ! ! RETURN
212 ! ENDIF
213
214 i = i + 1
215 END DO
216
217 ! limit zL to be within [-5,5]
218 ! IF (zL < -5 .OR. zL > 5) THEN
219 ! zL = MIN(5., MAX(-5., zL))
220 ! limit other output variables as well as z/L
221 ! L_MOD = zzd/zL
222 ! limit L_mod to 3.e4 for consistency with output
223 ! IF (ABS(L_MOD) > 3E4) L_MOD = L_MOD/ABS(L_MOD)*3E4
224 ! z0L = z0m/L_MOD
225 ! psim = stab_psi_mom(StabilityMethod, zL)
226 ! psimz0 = stab_psi_mom(StabilityMethod, z0L)
227
228 ! TS 01 Aug 2018:
229 ! TS: limit UStar and TStar to reasonable values
230 ! to prevent potential issues in other stability-related calcualtions
231 ! 02 Aug 2018: set a low limit at 0.15 m/s (Schumann 1988, BLM, DOI: 10.1007/BF00123019)
232 ! UStar = KUZ/(LOG(Zzd/z0m) - psim + psimz0)
233 ! UStar = MAX(0.15, UStar)
234 ! TStar = (-H/UStar)
235
236 ! IF (UStar < 0.0001) THEN !If u* still too small after iteration, then force quit simulation and write out error info
237 ! ! UStar=KUZ/(LOG(Zzd/z0m))
238 ! PRINT *, 'UStar', UStar, KUZ, (LOG(Zzd/z0m)), Zzd, z0m
239 ! CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] zl,L_MOD', zl, L_MOD, notUsedI)
240 ! CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] z0l,UStar', z0l, UStar, notUsedI)
241 ! CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] psim,psimz0', psim, psimz0, notUsedI)
242 ! CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] AVU1,log(zzd/z0m)', AVU1, LOG(zzd/z0m), notUsedI)
243
244 ! ! RETURN
245 ! ENDIF
246
247 ! TS 11 Feb 2021: limit UStar and TStar to reasonable ranges
248 ! under all conditions, min(UStar)==0.001 m s-1 (Jimenez et al 2012, MWR, https://doi.org/10.1175/mwr-d-11-00056.1
249 ustar = max(0.001, ustar)
250 ! under convective/unstable condition, min(UStar)==0.15 m s-1: (Schumann 1988, BLM, https://doi.org/10.1007/BF00123019)
251 IF (z0l < -neut_limit) ustar = max(0.15, ustar)
252
253 ! update associated variables
254 tstar = (-h/ustar)
255 l_mod = (ustar**2)/(g_t_k*tstar)
256 zl = zzd/l_mod
257
258 ! if ( L_MOD<-990 ) then
259 ! print*, 'L_MOD too low',L_MOD
260 ! print*, 1/(L_MOD-L_MOD)
261 !
262 ! end if
263
264 END SUBROUTINE cal_stab
265
266 !==================================================================
267 ! TS 20210208: restructure stability corrections funtions;
268 ! - the previous implementations are kept down below for cross-validation
269
270 ! integral form for momentum
271 FUNCTION stab_psi_mom(StabilityMethod, ZL) RESULT(psim)
272 IMPLICIT NONE
273
274 REAL(kind(1d0)) :: zl, psim
275 INTEGER :: stabilitymethod
276
277 SELECT CASE (stabilitymethod)
278 CASE (j12)
279 psim = psi_mom_j12(zl)
280
281 CASE (k75)
282 psim = psi_mom_k75(zl)
283
284 CASE (b71)
285 psim = psi_mom_b71(zl)
286
287 END SELECT
288
289 END FUNCTION
290
291 ! integral form for heat
292 FUNCTION stab_psi_heat(StabilityMethod, ZL) RESULT(psih)
293 IMPLICIT NONE
294
295 REAL(kind(1d0)) :: zl, psih
296 INTEGER :: stabilitymethod
297
298 SELECT CASE (stabilitymethod)
299 CASE (j12)
300 psih = psi_heat_j12(zl)
301
302 CASE (k75)
303 psih = psi_heat_k75(zl)
304
305 CASE (b71)
306 psih = psi_heat_b71(zl)
307
308 END SELECT
309
310 END FUNCTION
311
312 ! correction function for momentum
313 FUNCTION stab_phi_mom(StabilityMethod, ZL) RESULT(phim)
314 IMPLICIT NONE
315
316 REAL(kind(1d0)) :: zl, phim
317 INTEGER :: stabilitymethod
318
319 SELECT CASE (stabilitymethod)
320 CASE (j12)
321 phim = phi_mom_j12(zl)
322
323 CASE (k75)
324 phim = phi_mom_k75(zl)
325
326 CASE (b71)
327 phim = phi_mom_b71(zl)
328
329 END SELECT
330
331 END FUNCTION
332
333 ! correction function for heat
334 FUNCTION stab_phi_heat(StabilityMethod, ZL) RESULT(phih)
335 IMPLICIT NONE
336
337 REAL(kind(1d0)) :: zl, phih
338 INTEGER :: stabilitymethod
339
340 SELECT CASE (stabilitymethod)
341 CASE (j12)
342 phih = phi_heat_j12(zl)
343
344 CASE (k75)
345 phih = phi_heat_k75(zl)
346
347 CASE (b71)
348 phih = phi_heat_b71(zl)
349
350 END SELECT
351
352 END FUNCTION
353
354 ! ------------------------------------
355 ! Jimenez et al. (2012) forms for WRF
356 ! https://doi.org/10.1175/mwr-d-11-00056.1
357 FUNCTION psi_mom_j12(ZL) RESULT(psim)
358 ! Jimenez et al. (2012), eqn 17 and 18
359 REAL(kind(1d0)), PARAMETER :: a = 6.1
360 REAL(kind(1d0)), PARAMETER :: b = 2.5
361 REAL(kind(1d0)) :: zl, psim
362
363 IF (abs(zl) < neut_limit) THEN
364 psim = 0
365 ELSEIF (zl < -neut_limit) THEN
366 ! Jimenez et al. (2012), eqn 17
367 psim = psi_mom_g00(zl)
368 ELSEIF (zl > neut_limit) THEN
369 ! Jimenez et al. (2012), eqn 18:
370 ! Cheng and Brutsaert (2005)
371 psim = psi_mom_cb05(zl)
372 END IF
373
374 END FUNCTION
375
376 FUNCTION phi_mom_j12(ZL) RESULT(phim)
377 ! not directly provided by J12
378 ! equations below are derived from related original formulations
379 REAL(kind(1d0)) :: zl, phim
380
381 IF (abs(zl) < neut_limit) THEN
382 phim = 1
383 ELSEIF (zl < -neut_limit) THEN
384 ! derivative of eqn 17 in Jimenez et al. (2012)
385 phim = phi_mom_g00(zl)
386 ELSEIF (zl > neut_limit) THEN
387 ! derivative of eqn 18 in Jimenez et al. (2012)
388 phim = phi_mom_cb05(zl)
389 END IF
390
391 END FUNCTION
392
393 FUNCTION psi_heat_j12(ZL) RESULT(psih)
394 ! Jimenez et al. (2012), eqn 17 and 19
395 REAL(kind(1d0)) :: zl, psih
396
397 IF (abs(zl) < neut_limit) THEN
398 psih = 0
399 ELSEIF (zl < -neut_limit) THEN
400 ! Jimenez et al. (2012), eqn 17
401 psih = psi_heat_g00(zl)
402 ELSEIF (zl > neut_limit) THEN
403 ! Jimenez et al. (2012), eqn 19:
404 ! Cheng and Brutsaert (2005)
405 psih = psi_heat_cb05(zl)
406 END IF
407
408 END FUNCTION
409
410 FUNCTION phi_heat_j12(ZL) RESULT(phih)
411 REAL(kind(1d0)) :: zl, phih
412
413 IF (abs(zl) < neut_limit) THEN
414 phih = 1
415 ELSEIF (zl < -neut_limit) THEN
416 phih = phi_heat_g00(zl)
417 ELSEIF (zl > neut_limit) THEN
418 phih = phi_heat_cb05(zl)
419 END IF
420
421 END FUNCTION
422 ! ------------------------------------
423
424 ! ------------------------------------
425 ! Grachev et al. (2000) forms
426 ! https://doi.org/10.1023/a:1002452529672
427 ! !ONLY DEFINED FOR UNSTABLE CONDITIONS!
428 ! used in J12 under unstable conditions
429 FUNCTION psi_mom_g00(ZL) RESULT(psim)
430 ! Grachev et al. (2000), eqn 18
431 ! determination of `am` see sect. 4.1 of G00
432 REAL(kind(1d0)), PARAMETER :: am = 10
433 REAL(kind(1d0)) :: zl, psim
434 REAL(kind(1d0)) :: psim_k
435 REAL(kind(1d0)) :: psim_c
436
437 IF (abs(zl) < neut_limit) THEN
438 psim = 0
439 ELSEIF (zl < -neut_limit) THEN
440 ! Kansas part:
441 psim_k = psi_mom_b71(zl)
442 ! convective contribution:
443 psim_c = psi_conv(zl, am)
444 ! zL weighted sum:
445 psim = dot_product([psim_k, psim_c], [1d0, zl**2])
446 psim = psim/sum([1d0, zl**2])
447 END IF
448
449 END FUNCTION
450
451 FUNCTION psi_heat_g00(ZL) RESULT(psih)
452 ! Grachev et al. (2000), eqn 18
453 ! determination of `ah` see sect. 4.1 of G00
454 REAL(kind(1d0)), PARAMETER :: ah = 34
455
456 REAL(kind(1d0)) :: zl, psih
457 REAL(kind(1d0)) :: psih_k
458 REAL(kind(1d0)) :: psih_c
459
460 IF (abs(zl) < neut_limit) THEN
461 psih = 0
462 ELSEIF (zl < -neut_limit) THEN
463 ! Kansas part:
464 psih_k = psi_heat_b71(zl)
465 ! convective contribution:
466 psih_c = psi_conv(zl, ah)
467 ! zL weighted sum:
468 psih = dot_product([psih_k, psih_c], [1d0, zl**2])
469 psih = psih/sum([1d0, zl**2])
470
471 END IF
472
473 END FUNCTION
474
475 ! NB: the formulation of phi(zL) below is "correct" by definition:
476 ! phi=1-zL*d(Psi)/d(zL)
477 ! which however shows strange numeric behaviour:
478 ! a jumpy transition from phih_K to phih_C and dives into a negative range at some point
479 ! thus it's NOT adopted here
480 ! instead the following formulations phi_m/phi_h similar to the weighted psi_hm are used:
481
482 FUNCTION phi_mom_g00(ZL) RESULT(phim)
483 REAL(kind(1d0)), PARAMETER :: am = 10
484 REAL(kind(1d0)) :: zl, phim
485 REAL(kind(1d0)) :: phim_k ! Kansas part
486 ! REAL(KIND(1D0)):: psim_K ! Kansas part
487 REAL(kind(1d0)) :: phim_c ! convective part
488 ! REAL(KIND(1D0)):: psim_C ! convective part
489
490 IF (abs(zl) < neut_limit) THEN
491 phim = 1
492 ELSEIF (zl < -neut_limit) THEN
493 ! kansas
494 ! psim_K = psi_mom_B71(ZL)
495 phim_k = phi_mom_b71(zl)
496
497 ! convective
498 ! psim_C = psi_conv(ZL, am)
499 phim_c = phi_conv(zl, am)
500
501 ! by definition correct but not used due to bizzare numeric behaviour
502 ! phim = 1 - zL*dPsi_dzL_G00(zL, psim_K, phim_K, psim_C, phim_C)
503
504 ! weighted sum:
505 phim = dot_product([phim_k, phim_c], [1d0, zl**2])
506 phim = phim/sum([1d0, zl**2])
507
508 END IF
509
510 END FUNCTION
511
512 FUNCTION phi_heat_g00(ZL) RESULT(phih)
513 REAL(kind(1d0)), PARAMETER :: ah = 34
514 REAL(kind(1d0)) :: zl, phih
515 REAL(kind(1d0)) :: phih_k ! Kansas part
516 ! REAL(KIND(1D0)):: psih_K ! Kansas part
517 REAL(kind(1d0)) :: phih_c ! convective part
518 ! REAL(KIND(1D0)):: psih_C ! convective part
519
520 IF (abs(zl) < neut_limit) THEN
521 phih = 1
522 ELSEIF (zl < -neut_limit) THEN
523 ! kansas
524 ! psih_K = psi_heat_B71(ZL)
525 phih_k = phi_heat_b71(zl)
526
527 ! convective
528 ! psih_C = psi_conv(ZL, ah)
529 phih_c = phi_conv(zl, ah)
530
531 phih = dot_product([phih_k, phih_c], [1d0, zl**2])
532 phih = phih/sum([1d0, zl**2])
533
534 END IF
535
536 END FUNCTION
537
538 FUNCTION psi_conv(ZL, ax) RESULT(psiC)
539 ! Grachev et al. (2000), eqn 12
540 REAL(kind(1d0)) :: zl, psic, y, ax
541 REAL(kind(1d0)), PARAMETER :: pi = acos(-1.)
542
543 y = (1 - ax*zl)**(1/3.)
544 psic = 3./2*log(y**2 + y + 1/3.) - sqrt(3.)*atan(2*y + 1/sqrt(3.)) + pi/sqrt(3.)
545
546 END FUNCTION
547
548 FUNCTION phi_conv(ZL, ax) RESULT(phiC)
549 ! Grachev et al. (2000), eqn 15
550 REAL(kind(1d0)) :: zl, phic, ax, dpsi_dzl
551 REAL(kind(1d0)), PARAMETER :: pi = acos(-1.)
552
553 ! d(psi_conv)/d(zL)
554 dpsi_dzl = (9 + 4*(-3 + sqrt(3.))*ax*zl - 9*(1 - ax*zl)**(1/3.))/ &
555 (6.*zl*(2 - 2*ax*zl + (1 - ax*zl)**(2/3.)))
556
557 ! phi=1-zL*d(Psi)/d(zL)
558 phic = 1 - zl*dpsi_dzl
559
560 END FUNCTION
561
562 FUNCTION dpsi_dzl_g00(ZL, psiK, phiK, psiC, phiC) RESULT(dPsi)
563 ! Grachev et al. (2000), eqn 15
564 ! note: Psi is a zL-weighted sum of Kansas and convective components
565 ! derived from d(Psi)/d(zL)=d((psiK+zL**2*psiC)/(1+zL**2))/d(zL)
566
567 REAL(kind(1d0)) :: zl, dpsi, psic, phic, psik, phik
568 REAL(kind(1d0)) :: x_zl, x_psic, x_phic, x_psik, x_phik
569
570 ! Kansas part
571 x_psik = -(2*psik*zl)/(1 + zl**2)**2
572 x_phik = -phik/(zl*(1 + zl**2))
573
574 ! convective contribution
575 x_psic = (2*psic*zl)/(1 + zl**2)**2
576 x_phic = -(phic*zl)/(1 + zl**2)
577
578 ! zL related
579 x_zl = 1/zl
580
581 ! totoal
582 dpsi = x_psic + x_psik + x_phik + x_phic + x_zl
583
584 END FUNCTION
585
586 ! ------------------------------------
587
588 ! ------------------------------------
589 ! Cheng and Brutsaert (2005) forms
590 ! https://doi.org/10.1007/s10546-004-1425-4
591 ! !ONLY DEFINED FOR STABLE CONDITIONS!
592 FUNCTION psi_cb05(ZL, k1, k2) RESULT(psi)
593 ! Cheng and Brutsaert (2005), eqn 21/23
594 REAL(kind(1d0)) :: zl, psi, k1, k2
595 psi = -k1*log(zl + (1 + zl**k2)**(1/k2))
596 END FUNCTION
597
598 FUNCTION psi_mom_cb05(ZL) RESULT(psim)
599 ! Cheng and Brutsaert (2005), eqn 21
600 REAL(kind(1d0)), PARAMETER :: a = 6.1
601 REAL(kind(1d0)), PARAMETER :: b = 2.5
602 REAL(kind(1d0)) :: zl, psim
603
604 IF (abs(zl) < neut_limit) THEN
605 psim = 0
606 ELSEIF (zl > neut_limit) THEN
607 psim = psi_cb05(zl, a, b)
608 END IF
609
610 END FUNCTION
611
612 FUNCTION psi_heat_cb05(ZL) RESULT(psih)
613 ! Cheng and Brutsaert (2005), eqn 23
614 REAL(kind(1d0)), PARAMETER :: c = 5.3
615 REAL(kind(1d0)), PARAMETER :: d = 1.1
616 REAL(kind(1d0)) :: zl, psih
617
618 IF (abs(zl) < neut_limit) THEN
619 psih = 0
620 ELSEIF (zl > neut_limit) THEN
621 psih = psi_cb05(zl, c, d)
622 END IF
623
624 END FUNCTION
625
626 FUNCTION phi_cb05(ZL, k1, k2) RESULT(phi)
627 ! Cheng and Brutsaert (2005), eqn 22/24
628 REAL(kind(1d0)) :: zl, phi, k1, k2, zlk2
629 zlk2 = zl**k2
630 phi = 1 + k1*( &
631 (zl + zlk2*(1 + zlk2)**((1 - k2)/k2)) &
632 /(zl + (1 + zlk2)**(1/k2)) &
633 )
634 END FUNCTION
635
636 FUNCTION phi_mom_cb05(ZL) RESULT(phim)
637 ! Cheng and Brutsaert (2005), eqn 22
638 REAL(kind(1d0)), PARAMETER :: a = 6.1
639 REAL(kind(1d0)), PARAMETER :: b = 2.5
640 REAL(kind(1d0)) :: zl, phim
641
642 IF (abs(zl) < neut_limit) THEN
643 phim = 1
644 ELSEIF (zl > neut_limit) THEN
645 phim = phi_cb05(zl, a, b)
646 END IF
647
648 END FUNCTION
649
650 FUNCTION phi_heat_cb05(ZL) RESULT(phih)
651 ! Cheng and Brutsaert (2005), eqn 24
652 REAL(kind(1d0)), PARAMETER :: c = 5.3
653 REAL(kind(1d0)), PARAMETER :: d = 1.1
654 REAL(kind(1d0)) :: zl, phih, zld
655
656 IF (abs(zl) < neut_limit) THEN
657 phih = 1
658 ELSEIF (zl > neut_limit) THEN
659 zld = zl**d
660 phih = phi_cb05(zl, c, d)
661 END IF
662
663 END FUNCTION
664 ! ------------------------------------
665
666 ! ------------------------------------
667 ! Kondo (1975) forms
668 ! https://doi.org/10.1007/bf00232256
669 FUNCTION phi_mom_k75(ZL) RESULT(phim)
670 ! Kondo (1975) adopted by Campbell & Norman eqn 7.22 and 7.23 p 97
671 REAL(kind(1d0)) :: zl, phim
672
673 IF (abs(zl) < neut_limit) THEN
674 phim = 1
675 ELSEIF (zl < -neut_limit) THEN
676 phim = 1/(1 - 16*zl)**.25
677 ELSEIF (zl > neut_limit) THEN
678 phim = 1 + 6*zl/(1 + zl)
679 END IF
680
681 END FUNCTION
682
683 FUNCTION phi_heat_k75(ZL) RESULT(phih)
684 ! Kondo (1975) adopted by Campbell & Norman eqn 7.22 and 7.23 p 97
685 REAL(kind(1d0)) :: zl, phih
686
687 IF (abs(zl) < neut_limit) THEN
688 phih = 1
689 ELSEIF (zl < -neut_limit) THEN
690 phih = phi_mom_k75(zl)**2
691 ELSEIF (zl > neut_limit) THEN
692 phih = phi_mom_k75(zl)
693 END IF
694
695 END FUNCTION
696
697 FUNCTION psi_mom_k75(ZL) RESULT(psim)
698 ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 and 7.27 p 97
699 REAL(kind(1d0)) :: zl, psim
700
701 IF (abs(zl) < neut_limit) THEN
702 psim = 0
703 ELSEIF (zl < -neut_limit) THEN
704 psim = 0.6*psi_heat_k75(zl)
705 ELSEIF (zl > neut_limit) THEN
706 psim = psi_heat_k75(zl)
707 END IF
708
709 END FUNCTION
710
711 FUNCTION psi_heat_k75(ZL) RESULT(psih)
712 ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 and 7.27 p 97
713 REAL(kind(1d0)) :: zl, psih
714
715 IF (abs(zl) < neut_limit) THEN
716 psih = 0
717 ELSEIF (zl < -neut_limit) THEN
718 psih = (2)*log((1 + (1 - 16*zl)**0.5)/2)
719 ELSEIF (zl > neut_limit) THEN
720 psih = (-6)*log(1 + zl)
721 END IF
722
723 END FUNCTION
724 ! ------------------------------------
725
726 ! ------------------------------------
727 ! Businger et al. (1971) modified by Hogstrom (1988)
728 ! DOI: https://doi.org/10.1007/bf00119875
729 FUNCTION phi_mom_b71(ZL) RESULT(phim)
730 ! Businger et al. (1971)
731 ! modified by Hogstrom (1988): see Table VI
732 REAL(kind(1d0)) :: zl, phim
733
734 IF (abs(zl) < neut_limit) THEN
735 phim = 1
736 ELSEIF (zl < -neut_limit) THEN
737 phim = (1.-19.3*zl)**(-0.25)
738 ELSEIF (zl > neut_limit) THEN
739 phim = 1 + 6*zl
740 END IF
741
742 END FUNCTION
743
744 FUNCTION phi_heat_b71(ZL) RESULT(phih)
745 ! Businger et al. (1971)
746 ! modified by Hogstrom (1988): see Table VII
747 REAL(kind(1d0)) :: zl, phih
748
749 IF (abs(zl) < neut_limit) THEN
750 phih = 1
751 ELSEIF (zl < -neut_limit) THEN
752 phih = 0.95*(1.-11.6*zl)**(-0.5)
753 ELSEIF (zl > neut_limit) THEN
754 phih = 0.95 + 7.8*zl
755 END IF
756
757 END FUNCTION
758
759 FUNCTION psi_mom_b71(ZL) RESULT(psim)
760 ! integral form of phi_mom_B71
761 REAL(kind(1d0)), PARAMETER :: piover2 = acos(-1.)/2.
762 REAL(kind(1d0)) :: zl, psim, x, x2
763
764 IF (abs(zl) < neut_limit) THEN
765 psim = 0
766 ELSEIF (zl < -neut_limit) THEN
767 x = (1 - 19.3*zl)**(0.25)
768 x2 = log((1 + (x**2.))/2.)
769 psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
770
771 ELSEIF (zl > neut_limit) THEN
772 psim = (-6)*zl
773 END IF
774
775 END FUNCTION
776
777 FUNCTION psi_heat_b71(ZL) RESULT(psih)
778 ! integral form of phi_heat_B71
779 REAL(kind(1d0)) :: zl, psih, x
780
781 IF (abs(zl) < neut_limit) THEN
782 psih = 0
783 ELSEIF (zl < -neut_limit) THEN
784 x = 0.95*(1.-11.6*zl)**(0.5)
785 psih = 2*log((1 + x)/2)
786 ELSEIF (zl > neut_limit) THEN
787 IF (zl < 1) THEN
788 psih = (-7.8)*zl
789 ELSE
790 psih = (-7.8)*(1 + log(zl))
791 END IF
792 END IF
793
794 END FUNCTION
795 ! ------------------------------------
796
797 !==================================================================
798 ! previous implementations
799 ! TS 10 Feb 2021: commented after checing the new ones above match the ones below
800 ! FUNCTION stab_psi_mom(StabilityMethod, ZL) RESULT(psim)
801 ! ! StabilityMethod = 1-4 -
802 ! ! psim - stability FUNCTION for momentum
803 ! !Modified by LJ Mar 2010
804 ! !Input:Used stability method, stability (z-d)/L, zeta (either (z-d)/L or z0/L)
805
806 ! ! USE mod_z
807 ! ! USE mod_k
808
809 ! IMPLICIT NONE
810
811 ! REAL(KIND(1D0)):: piover2, psim, zl, x, x2
812 ! INTEGER ::StabilityMethod
813
814 ! psim = 0
815
816 ! PIOVER2 = ACOS(-1.)/2.
817 ! !PRINT*,StabilityMethod,zl,"stab_fn_mom:"
818 ! IF (ABS(zL) < neut_limit) THEN
819 ! psim = 0
820 ! ELSEIF (zL < -neut_limit) THEN !Unstable
821
822 ! IF (StabilityMethod == 1) THEN ! Jensen et al 1984 - Van Ulden & Holtslag (1985) p 1206
823 ! psim = ((1.-16.*zl)**0.25) - 1
824
825 ! ELSEIF (StabilityMethod == 2) THEN !Dyer (1974)(1-16z/L)**.25' k=0.41 mod. Hogstrom (1988)v15.2
826 ! X = (1.-(15.2*zl))**0.25
827 ! X2 = LOG((1 + (X**2.))/2.)
828 ! psim = (2.*LOG((1 + X)/2.)) + X2 - (2.*ATAN(X)) + PIOVER2
829
830 ! ELSEIF (StabilityMethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
831 ! psim = 0.6*(2)*LOG((1 + (1 - 16*zl)**0.5)/2)
832
833 ! ELSEIF (StabilityMethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
834 ! x = (1 - 19.3*zl)**(0.25) ! M Nag spotted the wrong exponent, TS corrected this from (1 - 19.3*zl_f)**(-0.25)
835 ! X2 = LOG((1 + (X**2.))/2.)
836 ! psim = (2.*LOG((1 + X)/2.)) + X2 - (2.*ATAN(X)) + PIOVER2
837
838 ! ! TS 20210208: commented out as not mentioned in manual
839 ! ! ELSEIF (StabilityMethod == 5) THEN ! Zilitinkevich & Chalikov (1968) modified Hogstrom (1988)
840 ! ! IF (zl >= -0.16) THEN
841 ! ! x = 1 + 1.38*zl
842 ! ! ELSE
843 ! ! x = 0.42*(-1)*zl**0.333
844 ! ! END IF
845 ! ! X2 = LOG((1 + (X**2.))/2.)
846 ! ! psim = (2.*LOG((1 + X)/2.)) + X2 - (2.*ATAN(X)) + PIOVER2
847
848 ! ! ELSEIF (StabilityMethod == 6) THEN ! Foken and Skeib (1983)
849 ! ! IF (zl >= 0.06) THEN
850 ! ! x = 1
851 ! ! ELSE
852 ! ! x = ((-1)*zl/0.06)**0.25
853 ! ! END IF
854 ! ! X2 = LOG((1 + (X**2.))/2.)
855 ! ! psim = (2.*LOG((1 + X)/2.)) + X2 - (2.*ATAN(X)) + PIOVER2
856
857 ! ! ELSEIF (StabilityMethod == 7) THEN ! Dyer & Bradley (1982) (1-28z/L)**.25' k=0.4
858 ! ! X = (1 - (28.*zl))**0.25 ! NT: changed + to - (bug -> checked reference)
859 ! ! X2 = LOG((1 + X**2.)/2.)
860 ! ! psim = (2.*LOG((1 + X)/2.)) + X2 - (2.*ATAN(X)) + PIOVER2
861 ! END IF
862
863 ! ELSEIF (zL > neut_limit) THEN !Stable
864
865 ! IF (StabilityMethod == 1) THEN !Dyer (1974) k=0.35 x=1+5*zl Mod. Hogstrom (1988)
866 ! psim = (-4.8)*zl
867 ! ELSEIF (StabilityMethod == 2) THEN !Van Ulden & Holtslag (1985) p 1206
868 ! IF (zl > 1000.) THEN
869 ! zl = 1000.
870 ! END IF
871 ! psim = (-17.*(1.-EXP(-0.29*zl)))
872 ! ELSEIF (StabilityMethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
873 ! psim = (-6)*LOG(1 + zl)
874 ! ELSEIF (StabilityMethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
875 ! ! psim=1+6*zl_f ! this is NOT the integral form but the stability function, TS 13 Jun 2017
876 ! psim = (-6)*zl ! this is the integral form, TS 13 Jun 2017
877
878 ! END IF
879 ! END IF
880 ! RETURN
881 ! END FUNCTION stab_psi_mom
882
883 ! FUNCTION stab_phi_mom(StabilityMethod, ZL) RESULT(phim)
884 ! ! StabilityMethod = 1-4 -
885 ! ! phi - stability FUNCTION for momentum
886 ! !Modified by NT May 2019 !!!!! check if all are correct!
887 ! !Input:Used stability method, stability (z-d)/L, zeta (either (z-d)/L or z0/L)
888
889 ! IMPLICIT NONE
890
891 ! REAL(KIND(1D0)):: phim, zl
892 ! INTEGER ::StabilityMethod
893
894 ! phim = 1
895
896 ! IF (ABS(zL) < neut_limit) THEN
897 ! phim = 1
898 ! ELSEIF (zL < -neut_limit) THEN !Unstable
899
900 ! IF (StabilityMethod == 1) THEN ! Jensen et al 1984 - Van Ulden & Holtslag (1985) p 1206&
901 ! phim = ((1.-16.*zl)**(-0.25))
902 ! ELSEIF (StabilityMethod == 2) THEN !Dyer (1974)(1-16z/L)**.25' k=0.41 mod. Hogstrom (1988)v15.2
903 ! phim = (1.-(15.2*zl))**(-0.25)
904 ! ELSEIF (StabilityMethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
905 ! phim = ((1.-16.*zl)**(-0.25))
906 ! ELSEIF (StabilityMethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
907 ! phim = (1.-19.3*zl)**(-0.25)
908
909 ! ! TS 20210208: methods 5–7 are commented as they are not mentioned in the manual
910 ! ! ELSEIF (StabilityMethod == 5) THEN ! Zilitinkevich & Chalikov (1968) modified Hogstrom (1988)
911 ! ! IF (zl >= -0.16) THEN
912 ! ! phim = 1 + 1.38*zl
913 ! ! ELSE
914 ! ! phim = 0.42*(-1)*zl*(-0.333)
915 ! ! END IF
916 ! ! ELSEIF (StabilityMethod == 6) THEN ! Foken and Skeib (1983)
917 ! ! IF (zl >= 0.06) THEN
918 ! ! phim = 1
919 ! ! ELSE
920 ! ! phim = ((-1)*zl/0.06)**(-0.25)
921 ! ! END IF
922 ! ! ELSEIF (StabilityMethod == 7) THEN ! Dyer & Bradley (1982) (1-28z/L)**.25' k=0.4
923 ! ! phim = (1.-(28.*zl))**(-0.25)
924 ! END IF
925
926 ! ELSEIF (zL > neut_limit) THEN !Stable
927
928 ! IF (StabilityMethod == 1) THEN !Dyer (1974) k=0.35 x=1+5*zl Mod. Hogstrom (1988)
929 ! phim = 1.+(4.8)*zl
930 ! ELSEIF (StabilityMethod == 2) THEN !Van Ulden & Holtslag (1985) p 1206 ! NT: have no function for phim
931 ! phim = 1.+(4.8)*zl
932 ! ELSEIF (StabilityMethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97 !!NT: checked
933 ! phim = 1.+6.*zl/(1.+zl) !!NT: checked reference and updated
934 ! ELSEIF (StabilityMethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
935 ! phim = 1 + 6*zl
936 ! END IF
937 ! END IF
938 ! RETURN
939 ! END FUNCTION stab_phi_mom
940 ! !_______________________________________________________________
941 ! !
942 ! ! psih - stability function for heat
943 ! FUNCTION stab_psi_heat(StabilityMethod, ZL) RESULT(psih)
944 ! ! USE mod_k
945 ! IMPLICIT NONE
946
947 ! REAL(KIND(1D0)):: zl, psih, x
948 ! INTEGER :: StabilityMethod
949
950 ! ! initialisation
951 ! psih = 0
952 ! x = 0
953
954 ! IF (ABS(zl) < neut_limit) THEN !Neutral
955 ! psih = 0
956 ! ELSEIF (zL < -neut_limit) THEN ! Unstable
957 ! IF (StabilityMethod == 3) THEN
958 ! !campbell & norman eqn 7.26
959 ! psih = (2)*LOG((1 + (1 - 16*zl)**0.5)/2)
960 ! ELSE
961 ! IF (StabilityMethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
962 ! x = 0.95*(1.-11.6*zl)**(0.5)
963 ! ELSEIF (StabilityMethod == 2) THEN ! Dyer 1974 X=(1.-(16.*ZL))**(0.5)modified Hosgstrom
964 ! x = 0.95*(1.-15.2*zl)**0.5
965
966 ! ! TS 20210208: methed=7 removed as not mentioned in manual
967 ! ! ELSEIF (StabilityMethod == 7) THEN
968 ! ! x = (1 - (28.*ZL))**0.25
969 ! END IF
970 ! ! psih = 2*LOG((1 + x**2)/2) ! NT: do not think this is correct
971 ! psih = 2*LOG((1 + x)/2)
972 ! END IF
973
974 ! ELSE IF (zL > neut_limit) THEN !Stable
975 ! IF (StabilityMethod == 3) THEN
976 ! !campbell & norman eqn 7.26
977 ! psih = (-6)*LOG(1 + zl)
978 ! ELSE
979 ! IF (zL <= 1) THEN ! weak/moderate stable
980 ! IF (StabilityMethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
981 ! psih = (-7.8)*zl
982 ! ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
983 ! psih = (-4.5)*zl
984 ! END IF
985 ! ELSE
986 ! ! adopt the form as Brutsaert (1982) eqn 4.58. but following the coeffs. of the above eqns
987 ! IF (StabilityMethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
988 ! psih = (-7.8)*(1 + LOG(zl))
989 ! ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
990 ! psih = (-4.5)*(1 + LOG(zl))
991 ! END IF
992 ! END IF
993 ! END IF
994
995 ! END IF
996
997 ! RETURN
998 ! END FUNCTION stab_psi_heat
999 ! !_______________________________________________________________
1000 ! !
1001 ! ! Phi - stability function for heat !!!!NT CHECK!!!!
1002 ! FUNCTION stab_phi_heat(StabilityMethod, ZL) RESULT(phih)
1003
1004 ! IMPLICIT NONE
1005
1006 ! REAL(KIND(1D0)):: zl, phih
1007 ! INTEGER :: StabilityMethod
1008
1009 ! phih = 1
1010
1011 ! IF (ABS(zl) < neut_limit) THEN !Neutral
1012 ! phih = 1
1013 ! ELSEIF (zL < -neut_limit) THEN ! Unstable
1014 ! IF (StabilityMethod == 3) THEN
1015 ! !campbell & norman eqn 7.26
1016 ! phih = (1.-16.*zl)**(-0.5)
1017 ! ELSEIF (StabilityMethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
1018 ! phih = 0.95*(1.-11.6*zl)**(-0.5)
1019 ! ELSEIF (StabilityMethod == 7) THEN
1020 ! phih = (1 - (28.*ZL))**(-0.25)
1021 ! ELSEIF (StabilityMethod == 2) THEN ! Dyer 1974 X=(1.-(16.*ZL))**(0.5)modified Hosgstrom
1022 ! phih = 0.95*(1.-15.2*zl)**(-0.5)
1023 ! END IF
1024
1025 ! ELSE IF (zL > neut_limit) THEN !Stable
1026 ! IF (StabilityMethod == 3) THEN
1027 ! !campbell & norman eqn 7.26
1028 ! phih = 1 + 6*zl/(1 + zL)
1029 ! ELSE
1030 ! IF (StabilityMethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
1031 ! phih = 0.95 + 7.8*zl ! this is NOT the integral form but the stability function, TS 13 Jun 2017
1032 ! ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
1033 ! phih = 1.+4.5*zl
1034 ! END IF
1035
1036 ! END IF
1037 ! END IF
1038
1039 ! RETURN
1040 ! END FUNCTION stab_phi_heat
1041
1042 !--------------------------------------------------------------------------------
1043 ! psys - roughness sublayer correction psi_*
1044 !
1045 ! Garratt (1980) QJRMS Appendix 1 p 815/816
1046
1047 ! FUNCTION stab_fn_rou(z, zstar) RESULT(psys)
1048 ! IMPLICIT NONE
1049 ! REAL(KIND(1d0))::alpha, zeta, z, psys, zstar, alpha1
1050 ! ! z wind speed height - z_d
1051 ! ! zstar height of the roughness sublayer
1052 ! ! eqn (a4) using alpha=0.5 alpha1=0.7
1053 ! alpha = 0.5
1054 ! alpha1 = 0.7
1055 ! zeta = z/zstar
1056 ! psys = (alpha - 1)*LOG(zeta) - (alpha*alpha1)*(1 - zeta) - (alpha*alpha1**2) &
1057 ! *(1 - zeta**2)/6.-(alpha*alpha1**3)*(1 - zeta**3)/24.
1058 ! RETURN
1059 ! END FUNCTION stab_fn_rou
1060
1061END MODULE atmmoiststab_module
atmmoiststab_module
Definition suews_phys_atmmoiststab.f95:1
atmmoiststab_module::phi_heat_k75
real(kind(1d0)) function phi_heat_k75(zl)
Definition suews_phys_atmmoiststab.f95:684
atmmoiststab_module::phi_cb05
real(kind(1d0)) function phi_cb05(zl, k1, k2)
Definition suews_phys_atmmoiststab.f95:627
atmmoiststab_module::psi_conv
real(kind(1d0)) function psi_conv(zl, ax)
Definition suews_phys_atmmoiststab.f95:539
atmmoiststab_module::dpsi_dzl_g00
real(kind(1d0)) function dpsi_dzl_g00(zl, psik, phik, psic, phic)
Definition suews_phys_atmmoiststab.f95:563
atmmoiststab_module::psi_heat_k75
real(kind(1d0)) function psi_heat_k75(zl)
Definition suews_phys_atmmoiststab.f95:712
atmmoiststab_module::phi_mom_k75
real(kind(1d0)) function phi_mom_k75(zl)
Definition suews_phys_atmmoiststab.f95:670
atmmoiststab_module::cal_atmmoist
subroutine cal_atmmoist(temp_c, press_hpa, avrh, dectime, lv_j_kg, lvs_j_kg, es_hpa, ea_hpa, vpd_hpa, vpd_pa, dq, dens_dry, avcp, air_dens)
Definition suews_phys_atmmoiststab.f95:20
atmmoiststab_module::phi_mom_j12
real(kind(1d0)) function phi_mom_j12(zl)
Definition suews_phys_atmmoiststab.f95:377
atmmoiststab_module::k
real(kind(1d0)), parameter k
Definition suews_phys_atmmoiststab.f95:4
atmmoiststab_module::j12
integer, parameter j12
Definition suews_phys_atmmoiststab.f95:8
atmmoiststab_module::psi_heat_j12
real(kind(1d0)) function psi_heat_j12(zl)
Definition suews_phys_atmmoiststab.f95:394
atmmoiststab_module::stab_psi_mom
real(kind(1d0)) function stab_psi_mom(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:272
atmmoiststab_module::psi_mom_b71
real(kind(1d0)) function psi_mom_b71(zl)
Definition suews_phys_atmmoiststab.f95:760
atmmoiststab_module::b71
integer, parameter b71
Definition suews_phys_atmmoiststab.f95:10
atmmoiststab_module::psi_mom_k75
real(kind(1d0)) function psi_mom_k75(zl)
Definition suews_phys_atmmoiststab.f95:698
atmmoiststab_module::phi_mom_g00
real(kind(1d0)) function phi_mom_g00(zl)
Definition suews_phys_atmmoiststab.f95:483
atmmoiststab_module::phi_mom_b71
real(kind(1d0)) function phi_mom_b71(zl)
Definition suews_phys_atmmoiststab.f95:730
atmmoiststab_module::psi_mom_cb05
real(kind(1d0)) function psi_mom_cb05(zl)
Definition suews_phys_atmmoiststab.f95:599
atmmoiststab_module::cal_stab
subroutine cal_stab(stabilitymethod, zzd, z0m, zdm, avu1, temp_c, qh_init, avdens, avcp, l_mod, tstar, ustar, zl)
Definition suews_phys_atmmoiststab.f95:125
atmmoiststab_module::psi_mom_g00
real(kind(1d0)) function psi_mom_g00(zl)
Definition suews_phys_atmmoiststab.f95:430
atmmoiststab_module::phi_heat_b71
real(kind(1d0)) function phi_heat_b71(zl)
Definition suews_phys_atmmoiststab.f95:745
atmmoiststab_module::stab_phi_mom
real(kind(1d0)) function stab_phi_mom(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:314
atmmoiststab_module::grav
real(kind(1d0)), parameter grav
Definition suews_phys_atmmoiststab.f95:5
atmmoiststab_module::psi_heat_g00
real(kind(1d0)) function psi_heat_g00(zl)
Definition suews_phys_atmmoiststab.f95:452
atmmoiststab_module::psi_mom_j12
real(kind(1d0)) function psi_mom_j12(zl)
Definition suews_phys_atmmoiststab.f95:358
atmmoiststab_module::phi_heat_g00
real(kind(1d0)) function phi_heat_g00(zl)
Definition suews_phys_atmmoiststab.f95:513
atmmoiststab_module::phi_heat_j12
real(kind(1d0)) function phi_heat_j12(zl)
Definition suews_phys_atmmoiststab.f95:411
atmmoiststab_module::stab_phi_heat
real(kind(1d0)) function stab_phi_heat(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:335
atmmoiststab_module::psi_cb05
real(kind(1d0)) function psi_cb05(zl, k1, k2)
Definition suews_phys_atmmoiststab.f95:593
atmmoiststab_module::psi_heat_cb05
real(kind(1d0)) function psi_heat_cb05(zl)
Definition suews_phys_atmmoiststab.f95:613
atmmoiststab_module::phi_mom_cb05
real(kind(1d0)) function phi_mom_cb05(zl)
Definition suews_phys_atmmoiststab.f95:637
atmmoiststab_module::neut_limit
real(kind(1d0)), parameter neut_limit
Definition suews_phys_atmmoiststab.f95:3
atmmoiststab_module::k75
integer, parameter k75
Definition suews_phys_atmmoiststab.f95:9
atmmoiststab_module::phi_conv
real(kind(1d0)) function phi_conv(zl, ax)
Definition suews_phys_atmmoiststab.f95:549
atmmoiststab_module::phi_heat_cb05
real(kind(1d0)) function phi_heat_cb05(zl)
Definition suews_phys_atmmoiststab.f95:651
atmmoiststab_module::psi_heat_b71
real(kind(1d0)) function psi_heat_b71(zl)
Definition suews_phys_atmmoiststab.f95:778
atmmoiststab_module::stab_psi_heat
real(kind(1d0)) function stab_psi_heat(stabilitymethod, zl)
Definition suews_phys_atmmoiststab.f95:293
meteo::spec_hum_def
real(kind(1d0)) function spec_hum_def(vpd_hpa, press_hpa)
Definition suews_util_meteo.f95:236
meteo::spec_heat_beer
real(kind(1d0)) function spec_heat_beer(temp_c, rh, rho_v, rho_d)
Definition suews_util_meteo.f95:245
meteo::lat_vap
real(kind(1d0)) function lat_vap(temp_c, ea_hpa, press_hpa, cp, dectime)
Definition suews_util_meteo.f95:278
meteo::sat_vap_press_x
real(kind(1d0)) function sat_vap_press_x(temp_c, press_hpa, from, dectime)
Definition suews_util_meteo.f95:149
meteo::lat_vapsublim
real(kind(1d0)) function lat_vapsublim(temp_c, ea_hpa, press_hpa, cp)
Definition suews_util_meteo.f95:349
errorhint
subroutine errorhint(errh, problemfile, value, value2, valuei)
Definition suews_ctrl_error.f95:2
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0