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)
19 REAL(KIND(1d0))::vap_dens
21 REAL(KIND(1d0)),
INTENT(in):: &
25 REAL(KIND(1d0)),
INTENT(out):: &
37 REAL(KIND(1d0)),
PARAMETER:: &
45 gas_ct_dry = 8.31451/0.028965, &
46 gas_ct_wv = 8.31451/0.0180153
54 ea_hpa = avrh/100*es_hpa
58 vpd_hpa = es_hpa - ea_hpa
59 vpd_pa = (es_hpa*100) - (ea_hpa*100)
64 vap_dens = (ea_hpa*100/((temp_c + 273.16)*gas_ct_wv))
67 dens_dry = ((press_hpa - ea_hpa)*100)/(gas_ct_dry*(273.16 + temp_c))
70 air_dens = (press_hpa*100)/(gas_ct_dry*(temp_c + 273.16))
76 lv_j_kg =
lat_vap(temp_c, ea_hpa, press_hpa, avcp, dectime)
79 IF (temp_c < 0.000)
THEN 85 IF (press_hpa < 900)
THEN 86 CALL errorhint(46,
'Function LUMPS_cal_AtmMoist', press_hpa, -55.55d0, -55)
105 StabilityMethod, &! input
106 zzd, & !Active measurement height (meas. height-displac. height)
120 INTEGER,
INTENT(in):: StabilityMethod
122 REAL(KIND(1d0)),
INTENT(in)::zzd
123 REAL(KIND(1d0)),
INTENT(in)::z0m
124 REAL(KIND(1d0)),
INTENT(in)::zdm
125 REAL(KIND(1d0)),
INTENT(in)::avU1
126 REAL(KIND(1d0)),
INTENT(in)::Temp_C
127 REAL(KIND(1d0)),
INTENT(in)::QH_init
128 REAL(KIND(1d0)),
INTENT(in)::avdens
129 REAL(KIND(1d0)),
INTENT(in)::avcp
131 REAL(KIND(1d0)),
INTENT(out)::L_MOD
132 REAL(KIND(1d0)),
INTENT(out)::TStar
133 REAL(KIND(1d0)),
INTENT(out)::UStar
134 REAL(KIND(1d0)),
INTENT(out)::zL
137 REAL(KIND(1d0))::G_T_K, &
145 REAL(KIND(1d0)),
PARAMETER :: &
148 INTEGER,
PARAMETER :: notUsedI = -55
152 LOGICAL :: debug = .false.
155 h_init = qh_init/(avdens*avcp)
157 IF (debug)
WRITE (*, *) stabilitymethod, z0m, avu1, h_init, ustar, l_mod
158 g_t_k = (grav/(temp_c + 273.16))*k
161 'Windspeed Ht too low relative to zdm [Stability calc]- values [z-zdm, zdm]', &
164 ustar = kuz/log(zzd/z0m)
165 IF (abs(h_init) < 0.001)
THEN 171 l_mod = (ustar**2)/(g_t_k*tstar)
173 IF (log(zzd/z0m) < 0.001000)
THEN 175 CALL errorhint(17,
'In stability subroutine, (z-zd) < z0.', zzd, z0m, notusedi)
179 DO WHILE ((abs(lold - l_mod) > 0.01) .AND. (i < 330))
193 ustar = kuz/(log(zzd/z0m) - psim + psimz0)
196 l_mod = (ustar**2)/(g_t_k*tstar)
211 IF (abs(l_mod) > 1e6) l_mod = l_mod/abs(l_mod)*1e6
215 zl = min(5., max(-5., zl))
223 ustar = max(0.15, kuz/(log(zzd/z0m) - psim + psimz0))
231 IF (ustar < 0.0001)
THEN 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)
268 REAL(KIND(1d0)):: piover2, psim, zl, x, x2
269 INTEGER ::StabilityMethod
273 piover2 = acos(-1.)/2.
279 IF (stabilitymethod == 1)
THEN 280 psim = ((1.-16.*zl)**0.25) - 1
281 ELSEIF (stabilitymethod == 2)
THEN 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 286 psim = 0.6*(2)*log((1 + (1 - 16*zl)**0.5)/2)
287 ELSEIF (stabilitymethod == 4)
THEN 288 x = (1 - 19.3*zl)**(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 292 x = (1 - (28.*zl))**0.25
293 x2 = log((1 + x**2.)/2.)
294 psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
295 ELSEIF (stabilitymethod == 5)
THEN 296 IF (zl >= -0.16)
THEN 299 x = 0.42*(-1)*zl**0.333
301 x2 = log((1 + (x**2.))/2.)
302 psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
304 ELSEIF (stabilitymethod == 6)
THEN 308 x = ((-1)*zl/0.06)**0.25
310 x2 = log((1 + (x**2.))/2.)
311 psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
316 IF (stabilitymethod == 1)
THEN 318 ELSEIF (stabilitymethod == 2)
THEN 322 psim = (-17.*(1.-exp(-0.29*zl)))
323 ELSEIF (stabilitymethod == 4)
THEN 326 ELSEIF (stabilitymethod == 3)
THEN 327 psim = (-6)*log(1 + zl)
347 REAL(KIND(1d0)):: phim, zl
348 INTEGER ::StabilityMethod
356 IF (stabilitymethod == 1)
THEN 357 phim = ((1.-16.*zl)**(-0.25))
358 ELSEIF (stabilitymethod == 2)
THEN 359 phim = (1.-(15.2*zl))**(-0.25)
360 ELSEIF (stabilitymethod == 3)
THEN 361 phim = ((1.-16.*zl)**(-0.25))
362 ELSEIF (stabilitymethod == 4)
THEN 363 phim = (1.-19.*zl)**(-0.25)
364 ELSEIF (stabilitymethod == 7)
THEN 365 phim = (1.-(28.*zl))**(-0.25)
366 ELSEIF (stabilitymethod == 5)
THEN 367 IF (zl >= -0.16)
THEN 370 phim = 0.42*(-1)*zl*(-0.333)
372 ELSEIF (stabilitymethod == 6)
THEN 376 phim = ((-1)*zl/0.06)**(-0.25)
382 IF (stabilitymethod == 1)
THEN 384 ELSEIF (stabilitymethod == 2)
THEN 386 ELSEIF (stabilitymethod == 4)
THEN 388 ELSEIF (stabilitymethod == 3)
THEN 389 phim = 1.+6.*zl/(1.+zl)
406 REAL(KIND(1d0)):: zl, psih, x
407 INTEGER :: StabilityMethod
416 IF (stabilitymethod == 3)
THEN 418 psih = (2)*log((1 + (1 - 16*zl)**0.5)/2)
421 IF (stabilitymethod == 4)
THEN 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 426 x = 0.95*(1.-15.2*zl)**0.5
429 psih = 2*log((1 + x)/2)
434 IF (stabilitymethod == 4)
THEN 441 IF (stabilitymethod == 4)
THEN 442 psih = (-7.8)*(1 + log(zl))
444 psih = (-4.5)*(1 + log(zl))
464 REAL(KIND(1d0)):: zl, phih
465 INTEGER :: StabilityMethod
472 IF (stabilitymethod == 3)
THEN 474 phih = (1.-16.*zl)**(-0.5)
477 IF (stabilitymethod == 4)
THEN 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 482 phih = 0.95*(1.-15.2*zl)**(-0.5)
487 IF (stabilitymethod == 4)
THEN 506 REAL(KIND(1d0))::alpha, zeta, z, psys, zstar, alpha1
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.
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)