SUEWS API Site
Documentation of SUEWS source code
suews_phys_atmmoiststab.f95
Go to the documentation of this file.
2  IMPLICIT NONE
3  REAL(KIND(1d0)), PARAMETER :: neut_limit = 1.e-5 !Limit for neutral stability
4 
5 CONTAINS
6  !.c!! For Lumps Version 2 - no stability calculations
7  ! Latent heat of sublimation when air temperature below zero added. LJ Nov 2012
8  ! explict interface added to all subroutines, TS 08 Aug 2017
9  SUBROUTINE cal_atmmoist( &
10  Temp_C, Press_hPa, avRh, dectime, &! input:
11  lv_J_kg, lvS_J_kg, &! output:
12  es_hPa, Ea_hPa, VPd_hpa, VPD_Pa, dq, dens_dry, avcp, air_dens)
13 
14  USE meteo, ONLY: &
17 
18  IMPLICIT NONE
19  REAL(KIND(1d0))::vap_dens
20 
21  REAL(KIND(1d0)), INTENT(in):: &
22  Temp_C, &
23  Press_hPa, &
24  avRh, dectime
25  REAL(KIND(1d0)), INTENT(out):: &
26  lv_J_kg, &!Latent heat of vaporization in [J kg-1]
27  lvS_J_kg, &!Latent heat of sublimation in J/kg
28  es_hPa, &!Saturation vapour pressure over water in hPa
29  Ea_hPa, &!Vapour pressure of water in hPa
30  VPd_hpa, & !vapour pressure deficit in hPa
31  VPD_Pa, & !vapour pressure deficit in Pa
32  dq, &!Specific humidity deficit in g/kg
33  dens_dry, & !Vap density or absolute humidity (kg/m3)
34  avcp, &!specific heat capacity in J kg-1 K-1
35  air_dens!Air density in kg/m3
36 
37  REAL(KIND(1d0)), PARAMETER:: &
38  ! comp = 0.9995, &
39  ! epsil = 0.62197,& !ratio molecular weight of water vapor/dry air (kg/mol/kg/mol)
40  ! epsil_gkg = 621.97, & !ratio molecular weight of water vapor/dry air in g/kg
41  ! dry_gas = 8.31451,& !Dry gas constant (J/k/mol)
42  ! gas_ct_wat = 461.05,& !Gas constant for water (J/kg/K)
43  ! molar = 0.028965,& !Dry air molar fraction in kg/mol
44  ! molar_wat_vap = 0.0180153,& !Molar fraction of water vapor in kg/mol
45  gas_ct_dry = 8.31451/0.028965, & !j/kg/k=dry_gas/molar
46  gas_ct_wv = 8.31451/0.0180153 !j/kg/k=dry_gas/molar_wat_vap
47  ! waterDens = 999.8395 !Density of water in 0 cel deg
48  INTEGER::from = 1
49 
50  !Saturation vapour pressure over water in hPa
51  es_hpa = sat_vap_press_x(temp_c, press_hpa, from, dectime) ! dectime is more or less unnecessary here
52 
53  !Vapour pressure of water in hPa
54  ea_hpa = avrh/100*es_hpa
55 
56  ! if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'%',Temp_C
57 
58  vpd_hpa = es_hpa - ea_hpa !vapour pressure deficit in hPa
59  vpd_pa = (es_hpa*100) - (ea_hpa*100)!vapour pressure deficit in Pa
60 
61  dq = spec_hum_def(vpd_hpa, press_hpa) !Specific humidity deficit in g/kg
62 
63  !Vap density or absolute humidity (kg/m3)
64  vap_dens = (ea_hpa*100/((temp_c + 273.16)*gas_ct_wv))
65 
66  !density Dry Air Beer(1990) kg/m3
67  dens_dry = ((press_hpa - ea_hpa)*100)/(gas_ct_dry*(273.16 + temp_c))
68 
69  !Air density in kg/m3
70  air_dens = (press_hpa*100)/(gas_ct_dry*(temp_c + 273.16))
71 
72  !Calculate specific heat capacity in J kg-1 K-1
73  avcp = spec_heat_beer(temp_c, avrh, vap_dens, dens_dry)
74 
75  !Latent heat of vaporization in [J kg-1]
76  lv_j_kg = lat_vap(temp_c, ea_hpa, press_hpa, avcp, dectime)
77 
78  !Latent heat of sublimation in J/kg
79  IF (temp_c < 0.000) THEN
80  lvs_j_kg = lat_vapsublim(temp_c, ea_hpa, press_hpa, avcp)
81  ELSE
82  lvs_j_kg = 2834000
83  ENDIF
84  !if(debug)write(*,*)lv_J_kg,Temp_C,'lv2'
85  IF (press_hpa < 900) THEN
86  CALL errorhint(46, 'Function LUMPS_cal_AtmMoist', press_hpa, -55.55d0, -55)
87  ENDIF
88  RETURN
89  END SUBROUTINE cal_atmmoist
90 
91  !.c!! For Lumps Version 2 - no stability calculations
92  !==========================================================
93  ! Last change:
94  ! TS 08 Aug 2017: added explicit interface
95  ! TS 13 Jun 2017: corrections to the integral of stability functions
96  ! MH 12 Apr 2017: Stable limit to exit do-loop
97  ! LJ 25 Nov 2014: Limits for L
98  ! LJ 19 Feb 2010
99  ! SG 27 Mar 2000 4:44 pm
100  ! ust - friction velocity
101  ! L - monin obukhov stability length
102  ! Van Ulden & Holtslag (1985) JCAM: 24: 1196-1207
103 
104  SUBROUTINE cal_stab( &
105  StabilityMethod, &! input
106  zzd, & !Active measurement height (meas. height-displac. height)
107  z0m, & !Aerodynamic roughness length
108  zdm, & !Displacement height
109  avu1, & !Average wind speed
110  temp_c, & !Air temperature
111  qh_init, & !sensible heat flux [W m-2]
112  avdens, & ! air density
113  avcp, & ! heat capacity of air
114  l_mod, & !Obukhov length! output:
115  tstar, & !T*
116  ustar, & !Friction velocity
117  zl)!Stability scale
119  IMPLICIT NONE
120  INTEGER, INTENT(in):: StabilityMethod
121 
122  REAL(KIND(1d0)), INTENT(in)::zzd !Active measurement height (meas. height-displac. height)
123  REAL(KIND(1d0)), INTENT(in)::z0m !Aerodynamic roughness length
124  REAL(KIND(1d0)), INTENT(in)::zdm !Displacement height
125  REAL(KIND(1d0)), INTENT(in)::avU1 !Average wind speed
126  REAL(KIND(1d0)), INTENT(in)::Temp_C !Air temperature
127  REAL(KIND(1d0)), INTENT(in)::QH_init !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
128  REAL(KIND(1d0)), INTENT(in)::avdens ! air density [kg m-3]
129  REAL(KIND(1d0)), INTENT(in)::avcp ! volumetric heat capacity [J m-3 K-1]
130 
131  REAL(KIND(1d0)), INTENT(out)::L_MOD!Obukhov length
132  REAL(KIND(1d0)), INTENT(out)::TStar!T*
133  REAL(KIND(1d0)), INTENT(out)::UStar!Friction velocity
134  REAL(KIND(1d0)), INTENT(out)::zL ! Stability scale
135  ! REAL(KIND(1d0)),INTENT(out)::psim !Stability function of momentum
136 
137  REAL(KIND(1d0))::G_T_K, &
138  KUZ, &
139  LOLD, &
140  psim, &
141  z0l, &
142  psimz0, &
143  H_init, &
144  h
145  REAL(KIND(1d0)), PARAMETER :: &
146  k = 0.4, & !Von Karman's contant
147  grav = 9.80665 !g - gravity - physics today august 1987
148  INTEGER, PARAMETER :: notUsedI = -55
149 
150  INTEGER :: i
151 
152  LOGICAL :: debug = .false.
153 
154  !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
155  h_init = qh_init/(avdens*avcp)
156 
157  IF (debug) WRITE (*, *) stabilitymethod, z0m, avu1, h_init, ustar, l_mod
158  g_t_k = (grav/(temp_c + 273.16))*k !gravity constant/(Temperature*Von Karman Constant)
159  kuz = k*avu1 !Von Karman constant*mean wind speed
160  IF (zzd < 0) CALL errorhint(32, &
161  'Windspeed Ht too low relative to zdm [Stability calc]- values [z-zdm, zdm]', &
162  zzd, zdm, notusedi)
163 
164  ustar = kuz/log(zzd/z0m) !Initial setting of u* and calc. of L_MOD (neutral situation)
165  IF (abs(h_init) < 0.001) THEN ! prevent zero TStar
166  h = 0.001
167  ELSE
168  h = h_init
169  END IF
170  tstar = (-h/ustar)
171  l_mod = (ustar**2)/(g_t_k*tstar)
172 
173  IF (log(zzd/z0m) < 0.001000) THEN
174  ! PRINT*, 1/(z0m-z0m)
175  CALL errorhint(17, 'In stability subroutine, (z-zd) < z0.', zzd, z0m, notusedi)
176  ENDIF
177  i = 1
178  lold = -999.
179  DO WHILE ((abs(lold - l_mod) > 0.01) .AND. (i < 330)) !NT: add error threshold !Iteration starts
180  lold = l_mod
181  zl = zzd/l_mod
182  z0l = z0m/l_mod !z0m roughness length
183 
184  ! IF (z L>2)THEN
185  ! CALL ErrorHint(73,'LUMPS_atmos_functions_stab.f95: stability scale (z/L), UStar',zL,UStar,notUsedI)
186  ! RETURN !MO-theory not necessarily valid above zL>2. Still causes problematic UStar values and correct limit might be 0.3.
187  ! !Needs more investigations.
188  ! END IF
189 
190  psim = stab_psi_mom(stabilitymethod, zl)
191  psimz0 = stab_psi_mom(stabilitymethod, z0l)
192 
193  ustar = kuz/(log(zzd/z0m) - psim + psimz0) !Friction velocity in non-neutral situation
194 
195  tstar = (-h/ustar)
196  l_mod = (ustar**2)/(g_t_k*tstar)
197 
198  ! IF(UStar<0.001000)THEN !If u* too small
199  ! UStar=KUZ/(LOG(Zzd/z0m))
200  ! PRINT*, 'UStar info',UStar,KUZ,(LOG(Zzd/z0m)),Zzd,z0m
201  ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] zl,dectime',zl,dectime,notUsedI)
202  ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] z0l,UStar',z0l,UStar,notUsedI)
203  ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] psim,psimz0',psim,psimz0,notUsedI)
204  ! CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] AVU1,log(zzd/z0m)',AVU1,LOG(zzd/z0m),notUsedI)
205  !
206  ! ! RETURN
207  ! ENDIF
208 
209  i = i + 1
210  ENDDO
211  IF (abs(l_mod) > 1e6) l_mod = l_mod/abs(l_mod)*1e6
212 
213  ! limit zL to be within [-5,5]
214  ! IF (zL < -5 .OR. zL > 5) THEN
215  zl = min(5., max(-5., zl))
216  ! limit other output variables as well as z/L
217  l_mod = zzd/zl
218  z0l = z0m/l_mod
219  psim = stab_psi_mom(stabilitymethod, zl)
220  psimz0 = stab_psi_mom(stabilitymethod, z0l)
221  ! TS 01 Aug 2018: set a low limit at 0.15 m/s (Schumann 1987, BLM)
222  ! to prevent potential issues in other stability-related calcualtions
223  ustar = max(0.15, kuz/(log(zzd/z0m) - psim + psimz0))
224  tstar = (-h/ustar)
225  ! END IF
226 
227  ! TS: limit UStar and TStar to reasonable values
228  ! 02 Aug 2018: set a low limit at 0.15 m/s (Schumann 1987, BLM)
229  ! UStar = MAX(0.15, UStar)
230  tstar = (-h/ustar)
231  IF (ustar < 0.0001) THEN !If u* still too small after iteration, then force quit simulation and write out error info
232  ! UStar=KUZ/(LOG(Zzd/z0m))
233  print *, 'UStar', ustar, kuz, (log(zzd/z0m)), zzd, z0m
234  CALL errorhint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] zl,L_MOD', zl, l_mod, notusedi)
235  CALL errorhint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] z0l,UStar', z0l, ustar, notusedi)
236  CALL errorhint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] psim,psimz0', psim, psimz0, notusedi)
237  CALL errorhint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] AVU1,log(zzd/z0m)', avu1, log(zzd/z0m), notusedi)
238 
239  ! RETURN
240  ENDIF
241 
242  ! if ( L_MOD<-990 ) then
243  ! print*, 'L_MOD too low',L_MOD
244  ! print*, 1/(L_MOD-L_MOD)
245  !
246  ! end if
247 
248  END SUBROUTINE cal_stab
249 
250  !==================================================================
251 
252  FUNCTION stab_psi_mom(StabilityMethod, ZL) RESULT(psim)
253  ! StabilityMethod = 1-4 -
254  ! psim - stability FUNCTION for momentum
255  !Modified by LJ Mar 2010
256  !Input:Used stability method, stability (z-d)/L, zeta (either (z-d)/L or z0/L)
257 
258  ! USE mod_z
259  ! USE mod_k
260 
261  IMPLICIT NONE
262  ! REAL(KIND(1d0)), PARAMETER :: &
263  ! k=0.4,& !Von Karman's contant
264  ! k2=0.16,& !Power of Van Karman's contant
265  ! neut_limit = 0.001000 !Limit for neutral stability
266  ! notUsedI=-55
267 
268  REAL(KIND(1d0)):: piover2, psim, zl, x, x2
269  INTEGER ::StabilityMethod
270 
271  psim = 0
272 
273  piover2 = acos(-1.)/2.
274  !PRINT*,StabilityMethod,zl,"stab_fn_mom:"
275  IF (abs(zl) < neut_limit) THEN
276  psim = 0
277  ELSEIF (zl < -neut_limit) THEN !Unstable
278 
279  IF (stabilitymethod == 1) THEN ! Jensen et al 1984 - Van Ulden & Holtslag (1985) p 1206&
280  psim = ((1.-16.*zl)**0.25) - 1
281  ELSEIF (stabilitymethod == 2) THEN !Dyer (1974)(1-16z/L)**.25' k=0.41 mod. Hogstrom (1988)v15.2
282  x = (1.-(15.2*zl))**0.25
283  x2 = log((1 + (x**2.))/2.)
284  psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
285  ELSEIF (stabilitymethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
286  psim = 0.6*(2)*log((1 + (1 - 16*zl)**0.5)/2)
287  ELSEIF (stabilitymethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
288  x = (1 - 19.3*zl)**(0.25) ! M Nag spotted the wrong exponent, TS corrected this from (1 - 19.3*zl_f)**(-0.25)
289  x2 = log((1 + (x**2.))/2.)
290  psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
291  ELSEIF (stabilitymethod == 7) THEN ! Dyer & Bradley (1982) (1-28z/L)**.25' k=0.4
292  x = (1 - (28.*zl))**0.25 ! NT: changed + to - (bug -> checked reference)
293  x2 = log((1 + x**2.)/2.)
294  psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
295  ELSEIF (stabilitymethod == 5) THEN ! Zilitinkevich & Chalikov (1968) modified Hogstrom (1988)
296  IF (zl >= -0.16) THEN
297  x = 1 + 1.38*zl
298  ELSE
299  x = 0.42*(-1)*zl**0.333
300  ENDIF
301  x2 = log((1 + (x**2.))/2.)
302  psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
303 
304  ELSEIF (stabilitymethod == 6) THEN ! Foken and Skeib (1983)
305  IF (zl >= 0.06) THEN
306  x = 1
307  ELSE
308  x = ((-1)*zl/0.06)**0.25
309  ENDIF
310  x2 = log((1 + (x**2.))/2.)
311  psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
312  ENDIF
313 
314  ELSEIF (zl > neut_limit) THEN !Stable
315 
316  IF (stabilitymethod == 1) THEN !Dyer (1974) k=0.35 x=1+5*zl Mod. Hogstrom (1988)
317  psim = (-4.8)*zl
318  ELSEIF (stabilitymethod == 2) THEN !Van Ulden & Holtslag (1985) p 1206
319  IF (zl > 1000.) THEN
320  zl = 1000.
321  END IF
322  psim = (-17.*(1.-exp(-0.29*zl)))
323  ELSEIF (stabilitymethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
324  ! psim=1+6*zl_f ! this is NOT the integral form but the stability function, TS 13 Jun 2017
325  psim = (-6)*zl ! this is the integral form, TS 13 Jun 2017
326  ELSEIF (stabilitymethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
327  psim = (-6)*log(1 + zl)
328 
329  ENDIF
330  ENDIF
331  RETURN
332  END FUNCTION stab_psi_mom
333 
334  FUNCTION stab_phi_mom(StabilityMethod, ZL) RESULT(phim)
335  ! StabilityMethod = 1-4 -
336  ! phi - stability FUNCTION for momentum
337  !Modified by NT May 2019 !!!!! check if all are correct!
338  !Input:Used stability method, stability (z-d)/L, zeta (either (z-d)/L or z0/L)
339 
340  IMPLICIT NONE
341  ! REAL(KIND(1d0)), PARAMETER :: &
342  ! k=0.4,& !Von Karman's contant
343  ! k2=0.16,& !Power of Van Karman's contant
344  ! neut_limit = 0.001000 !Limit for neutral stability
345  ! notUsedI=-55
346 
347  REAL(KIND(1d0)):: phim, zl
348  INTEGER ::StabilityMethod
349 
350  phim = 0
351 
352  IF (abs(zl) < neut_limit) THEN
353  phim = 0
354  ELSEIF (zl < -neut_limit) THEN !Unstable
355 
356  IF (stabilitymethod == 1) THEN ! Jensen et al 1984 - Van Ulden & Holtslag (1985) p 1206&
357  phim = ((1.-16.*zl)**(-0.25))
358  ELSEIF (stabilitymethod == 2) THEN !Dyer (1974)(1-16z/L)**.25' k=0.41 mod. Hogstrom (1988)v15.2
359  phim = (1.-(15.2*zl))**(-0.25)
360  ELSEIF (stabilitymethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97
361  phim = ((1.-16.*zl)**(-0.25))
362  ELSEIF (stabilitymethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
363  phim = (1.-19.*zl)**(-0.25)
364  ELSEIF (stabilitymethod == 7) THEN ! Dyer & Bradley (1982) (1-28z/L)**.25' k=0.4
365  phim = (1.-(28.*zl))**(-0.25)
366  ELSEIF (stabilitymethod == 5) THEN ! Zilitinkevich & Chalikov (1968) modified Hogstrom (1988)
367  IF (zl >= -0.16) THEN
368  phim = 1 + 1.38*zl
369  ELSE
370  phim = 0.42*(-1)*zl*(-0.333)
371  ENDIF
372  ELSEIF (stabilitymethod == 6) THEN ! Foken and Skeib (1983)
373  IF (zl >= 0.06) THEN
374  phim = 1
375  ELSE
376  phim = ((-1)*zl/0.06)**(-0.25)
377  ENDIF
378  ENDIF
379 
380  ELSEIF (zl > neut_limit) THEN !Stable
381 
382  IF (stabilitymethod == 1) THEN !Dyer (1974) k=0.35 x=1+5*zl Mod. Hogstrom (1988)
383  phim = 1.+(4.8)*zl
384  ELSEIF (stabilitymethod == 2) THEN !Van Ulden & Holtslag (1985) p 1206 ! NT: have no function for phim
385  phim = 1.+(4.8)*zl
386  ELSEIF (stabilitymethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
387  phim = 1 + 6*zl
388  ELSEIF (stabilitymethod == 3) THEN ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97 !!NT: checked
389  phim = 1.+6.*zl/(1.+zl) !!NT: checked reference and updated
390  ENDIF
391  ENDIF
392  RETURN
393  END FUNCTION stab_phi_mom
394  !_______________________________________________________________
395  !
396  ! psih - stability function for heat
397  FUNCTION stab_psi_heat(StabilityMethod, ZL) RESULT(psih)
398  ! USE mod_k
399  IMPLICIT NONE
400  ! REAL(KIND(1d0)), PARAMETER :: &
401  ! k=0.4,& !Von Karman's contant
402  ! k2=0.16,& !Power of Van Karman's contant
403  ! neut_limit = 0.001000 !Limit for neutral stability
404  ! notUsedI=-55
405 
406  REAL(KIND(1d0)):: zl, psih, x
407  INTEGER :: StabilityMethod
408 
409  ! initialisation
410  psih = 0
411  x = 0
412 
413  IF (abs(zl) < neut_limit) THEN !Neutral
414  psih = 0
415  ELSEIF (zl < -neut_limit) THEN ! Unstable
416  IF (stabilitymethod == 3) THEN
417  !campbell & norman eqn 7.26
418  psih = (2)*log((1 + (1 - 16*zl)**0.5)/2)
419  ELSE
420 
421  IF (stabilitymethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
422  x = 0.95*(1.-11.6*zl)**(0.5)
423  ELSEIF (stabilitymethod == 7) THEN
424  x = (1 - (28.*zl))**0.25
425  ELSEIF (stabilitymethod == 2) THEN ! Dyer 1974 X=(1.-(16.*ZL))**(0.5)modified Hosgstrom
426  x = 0.95*(1.-15.2*zl)**0.5
427  ENDIF
428  ! psih = 2*LOG((1 + x**2)/2) ! NT: do not think this is correct
429  psih = 2*log((1 + x)/2)
430  ENDIF
431 
432  ELSE IF (zl > neut_limit) THEN !Stable
433  IF (zl <= 1) THEN ! weak/moderate stable
434  IF (stabilitymethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
435  psih = (-7.8)*zl
436  ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
437  psih = (-4.5)*zl
438  ENDIF
439  ELSE
440  ! adopt the form as Brutasert (1982) eqn 4.58. but following the coeffs. of the above eqns
441  IF (stabilitymethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
442  psih = (-7.8)*(1 + log(zl))
443  ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
444  psih = (-4.5)*(1 + log(zl))
445  ENDIF
446  END IF
447 
448  ENDIF
449 
450  RETURN
451  END FUNCTION stab_psi_heat
452  !_______________________________________________________________
453  !
454  ! Phi - stability function for heat !!!!NT CHECK!!!!
455  FUNCTION stab_phi_heat(StabilityMethod, ZL) RESULT(phih)
456  ! USE mod_k
457  IMPLICIT NONE
458  ! REAL(KIND(1d0)), PARAMETER :: &
459  ! k=0.4,& !Von Karman's contant
460  ! k2=0.16,& !Power of Van Karman's contant
461  ! neut_limit = 0.001000 !Limit for neutral stability
462  ! notUsedI=-55
463 
464  REAL(KIND(1d0)):: zl, phih
465  INTEGER :: StabilityMethod
466 
467  phih = 0
468 
469  IF (abs(zl) < neut_limit) THEN !Neutral
470  phih = 0
471  ELSEIF (zl < -neut_limit) THEN ! Unstable
472  IF (stabilitymethod == 3) THEN
473  !campbell & norman eqn 7.26
474  phih = (1.-16.*zl)**(-0.5)
475  ELSE
476 
477  IF (stabilitymethod == 4) THEN ! Businger et al (1971) modifed Hogstrom (1988)
478  phih = 0.95*(1.-11.6*zl)**(-0.5)
479  ELSEIF (stabilitymethod == 7) THEN
480  phih = (1 - (28.*zl))**(-0.25)
481  ELSEIF (stabilitymethod == 2) THEN ! Dyer 1974 X=(1.-(16.*ZL))**(0.5)modified Hosgstrom
482  phih = 0.95*(1.-15.2*zl)**(-0.5)
483  ENDIF
484  ENDIF
485 
486  ELSE IF (zl > neut_limit) THEN !Stable
487  IF (stabilitymethod == 4) THEN !Businger et al (1971) modifed Hogstrom (1988)
488  ! psih=0.95+(7.8*zl_f) ! this is NOT the integral form but the stability function, TS 13 Jun 2017
489  phih = 1.+7.8*zl ! this is the integral form, TS 13 Jun 2017
490  ELSE !Dyer (1974) psih=(-5)*ZL modifed Hogstrom (1988)
491  phih = 1.+4.5*zl
492  ENDIF
493 
494  ENDIF
495 
496  RETURN
497  END FUNCTION stab_phi_heat
498 
499  !--------------------------------------------------------------------------------
500  ! psys - roughness sublayer correction psi_*
501  !
502  ! Garratt (1980) QJRMS Appendix 1 p 815/816
503 
504  FUNCTION stab_fn_rou(z, zstar) RESULT(psys)
505  IMPLICIT NONE
506  REAL(KIND(1d0))::alpha, zeta, z, psys, zstar, alpha1
507  ! z wind speed height - z_d
508  ! zstar height of the roughness sublayer
509  ! eqn (a4) using alpha=0.5 alpha1=0.7
510  alpha = 0.5
511  alpha1 = 0.7
512  zeta = z/zstar
513  psys = (alpha - 1)*log(zeta) - (alpha*alpha1)*(1 - zeta) - (alpha*alpha1**2) &
514  *(1 - zeta**2)/6.-(alpha*alpha1**3)*(1 - zeta**3)/24.
515  RETURN
516  END FUNCTION stab_fn_rou
517 
518 END MODULE atmmoiststab_module
real(kind(1d0)) function stab_phi_mom(StabilityMethod, ZL)
real(kind(1d0)) function spec_hum_def(vpd_hPa, press_hPa)
real(kind(1d0)) function spec_heat_beer(Temp_C, rh, rho_v, rho_d)
real(kind(1d0)), parameter neut_limit
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)
subroutine cal_stab(StabilityMethod, zzd, z0m, zdm, avU1, Temp_C, QH_init, avdens, avcp, L_MOD, TStar, UStar, zL)
real(kind(1d0)) function stab_fn_rou(z, zstar)
real(kind(1d0)) function lat_vapsublim(Temp_C, Ea_hPa, Press_hPa, cp)
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) function stab_phi_heat(StabilityMethod, ZL)
real(kind(1d0)) function sat_vap_press_x(Temp_c, PRESS_hPa, from, dectime)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
real(kind(1d0)) function lat_vap(Temp_C, Ea_hPa, Press_hPa, cp, dectime)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)