8 INTEGER,
PARAMETER ::
nz = 30
14 L_MOD, sfr, planF, 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
31 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(in) ::sfr
32 REAL(KIND(1d0)),
INTENT(in):: zMeas
33 REAL(KIND(1d0)),
INTENT(in):: avU1
34 REAL(KIND(1d0)),
INTENT(in):: Temp_C
35 REAL(KIND(1d0)),
INTENT(in):: avRH
36 REAL(KIND(1d0)),
INTENT(in):: Press_hPa
37 REAL(KIND(1d0)),
INTENT(in):: L_MOD
38 REAL(KIND(1d0)),
INTENT(in):: avcp
39 REAL(KIND(1d0)),
INTENT(in):: lv_J_kg
40 REAL(KIND(1d0)),
INTENT(in):: avdens
41 REAL(KIND(1d0)),
INTENT(in):: qh
42 REAL(KIND(1d0)),
INTENT(in):: qe
43 REAL(KIND(1d0)),
INTENT(in):: Zh
44 REAL(KIND(1d0)),
INTENT(in):: z0m
45 REAL(KIND(1d0)),
INTENT(in):: zdm
46 REAL(KIND(1d0)),
INTENT(in):: planF
48 INTEGER,
INTENT(in)::StabilityMethod
50 REAL(KIND(1d0)),
INTENT(out):: T2_C
51 REAL(KIND(1d0)),
INTENT(out):: q2_gkg
52 REAL(KIND(1d0)),
INTENT(out):: U10_ms
53 REAL(KIND(1d0)),
INTENT(out):: RH2
55 INTEGER,
PARAMETER :: nz = 30
57 REAL(KIND(1d0)),
PARAMETER:: cd_tree = 1.2, &
62 pi = 4.*atan(1.0), r = 0.1, &
63 a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
66 REAL(KIND(1d0)),
INTENT(out),
DIMENSION(nz*4):: dataoutLineRSL
67 REAL(KIND(1d0)),
DIMENSION(nz):: psihatm_z
68 REAL(KIND(1d0)),
DIMENSION(nz):: psihath_z
69 REAL(KIND(1d0)),
DIMENSION(nz):: dif
71 REAL(KIND(1d0)),
DIMENSION(nz):: zarray
72 REAL(KIND(1d0)),
DIMENSION(nz):: dataoutLineURSL
73 REAL(KIND(1d0)),
DIMENSION(nz):: dataoutLineTRSL
74 REAL(KIND(1d0)),
DIMENSION(nz):: dataoutLineqRSL
79 REAL(KIND(1d0))::Lc_build, Lc_tree, Lc
81 REAL(KIND(1d0))::phim, psimz, psimZh, psimz0, psimza, phi_hatmZh, phi_hathZh, phimzp, phimz, phihzp, phihz, psihz, psihza
82 REAL(KIND(1d0))::betaHF, betaNL, beta, betaN2
84 REAL(KIND(1d0))::xxm1, xxm1_2, xxh1, xxh1_2, err, z01, dphi, dphih
85 REAL(KIND(1d0))::f, cm, c2, ch, c2h
86 REAL(KIND(1d0))::t_h, q_h
87 REAL(KIND(1d0))::TStar_RSL
88 REAL(KIND(1d0))::UStar_RSL
89 REAL(KIND(1d0))::sfr_zh
90 REAL(KIND(1d0))::sfr_tr
91 REAL(KIND(1d0))::L_MOD_RSL
92 REAL(KIND(1d0))::zH_RSL
94 REAL(KIND(1d0)),
parameter::Zh_min = 0.15
95 REAL(KIND(1d0)),
parameter::ratio_dz = 1.618
97 REAL(KIND(1d0))::qa_gkg, qStar_RSL
98 INTEGER :: I, z, it, idx_can, idx_za, idx_2m, idx_10m
115 stabilitymethod, zh, l_mod, sfr, planf, &
116 l_mod_rsl, zh_rsl, lc, beta, zd, z0, elm, scc, f)
120 IF (zh_rsl <= 2)
THEN 122 ELSE IF (zh_rsl <= 10)
THEN 133 dz = (zmeas - zh_rsl)/(nz - nz_can)
134 do i = nz_can + 1, nz
135 zarray(i) = zh_rsl + (i - nz_can)*dz
141 dif(z) = abs(zarray(z) - 2)
143 idx_2m = minloc(dif, dim=1)
147 dif(z) = abs(zarray(z) - 10)
149 idx_10m = minloc(dif, dim=1)
154 dif(z) = abs(zarray(z) - zh_rsl)
156 idx_can = minloc(dif, dim=1)
157 zarray(idx_can) = zh_rsl
162 dif(z) = abs(zarray(z) - zmeas)
164 idx_za = minloc(dif, dim=1)
165 zarray(idx_za) = zmeas
167 if (zh_rsl - zd < z0 .or. zh < zh_min)
then 185 xxm1 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
186 xxm1_2 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
188 xxh1 =
stab_phi_heat(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
189 xxh1_2 =
stab_phi_heat(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
191 phi_hatmzh = kappa/(2.*beta*xxm1)
192 phi_hathzh = kappa*scc/(2.*beta*xxh1)
194 dphih = xxh1_2 - xxh1
195 IF (phi_hatmzh > 1.)
THEN 199 c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa)
200 c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
202 cm = (1.-phi_hatmzh)*exp(c2/2.)
203 ch = (1.-phi_hathzh)*exp(c2h/2.)
208 psihatm_z = 0.*zarray
209 psihath_z = 0.*zarray
212 DO z = nz - 1, idx_can - 1, -1
213 phimz =
stab_phi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
214 phimzp =
stab_phi_mom(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
215 phihz =
stab_phi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
216 phihzp =
stab_phi_heat(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
218 psihatm_z(z) = psihatm_z(z + 1) + dz/2.*phimzp*(cm*exp(-1.*c2*beta*(zarray(z + 1) - zd)/elm)) &
219 /(zarray(z + 1) - zd)
220 psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
222 psihath_z(z) = psihath_z(z + 1) + dz/2.*phihzp*(ch*exp(-1.*c2h*beta*(zarray(z + 1) - zd)/elm)) &
223 /(zarray(z + 1) - zd)
224 psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
235 psimza =
stab_psi_mom(stabilitymethod, (zmeas - zd)/l_mod_rsl)
236 psihza =
stab_psi_heat(stabilitymethod, (zmeas - zd)/l_mod_rsl)
237 ustar_rsl = avu1*kappa/(log((zmeas - zd)/z0) - psimza + psimz0 + psihatm_z(nz))
238 tstar_rsl = -1.*(qh/(avcp*avdens))/ustar_rsl
239 qstar_rsl = -1.*(qe/lv_j_kg*avdens)/ustar_rsl
240 qa_gkg =
rh2qa(avrh/100, press_hpa, temp_c)
243 psimz =
stab_psi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
244 psihz =
stab_psi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
245 dataoutlineursl(z) = (log((zarray(z) - zd)/z0) - psimz + psimz0 - psihatm_z(z) + psihatm_z(idx_can))/kappa
246 dataoutlineursl(z) = (log((zarray(z) - zd)/z0) - psimz + psimz0 + psihatm_z(z))/kappa
248 dataoutlinetrsl(z) = (log((zarray(z) - zd)/(zmeas - zd)) - psihz + psihza + psihath_z(z) - psihath_z(idx_za))/kappa
249 dataoutlineqrsl(z) = (log((zarray(z) - zd)/(zmeas - zd)) - psihz + psihza + psihath_z(z) - psihath_z(idx_za))/kappa
255 t_h = scc*tstar_rsl/(beta*f)
256 q_h = scc*qstar_rsl/(beta*f)
258 dataoutlineursl(z) = dataoutlineursl(idx_can)*exp(beta*(zarray(z) - zh_rsl)/elm)
259 dataoutlinetrsl(z) = dataoutlinetrsl(idx_can) + (t_h*exp(beta*f*(zarray(z) - zh_rsl)/elm) - t_h)/tstar_rsl
260 dataoutlineqrsl(z) = dataoutlineqrsl(idx_can) + (q_h*exp(beta*f*(zarray(z) - zh_rsl)/elm) - q_h)/qstar_rsl
263 dataoutlineursl = dataoutlineursl*ustar_rsl
264 dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
265 dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
267 dataoutlinersl = (/zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl/)
273 t2_c = dataoutlinetrsl(idx_2m)
274 q2_gkg = dataoutlineqrsl(idx_2m)
275 u10_ms = dataoutlineursl(idx_10m)
277 rh2 =
qa2rh(q2_gkg, press_hpa, t2_c)
291 recursive function cal_psim_hat(StabilityMethod, z, zh_RSL, L_MOD_RSL, beta, Lc)
result(psim_hat_z)
295 integer,
intent(in) :: StabilityMethod
296 real(KIND(1D0)),
intent(in) :: z
297 real(KIND(1D0)),
intent(in) :: zh_RSL
298 real(KIND(1D0)),
intent(in) :: Lc
299 real(KIND(1D0)),
intent(in) :: beta
300 real(KIND(1D0)),
intent(in) :: L_MOD_RSL
303 real(KIND(1D0)) ::psim_hat_z
308 real(KIND(1D0)) ::phim_lc
309 real(KIND(1D0)) ::phim_z
310 real(KIND(1D0)) ::phim_zp
311 real(KIND(1D0)) ::psim_hat_zp
312 real(KIND(1D0)) ::elm
313 real(KIND(1D0)) ::xxm1
314 real(KIND(1D0)) ::xxm1_2
315 real(KIND(1D0)) ::dphi
316 real(KIND(1D0)) ::phi_hatmZh
320 REAL(KIND(1d0)),
PARAMETER::kappa = 0.40
321 REAL(KIND(1d0)),
PARAMETER::dz = 0.1
330 zd = zh_rsl - (beta**2.)*lc
336 phim_z =
stab_phi_mom(stabilitymethod, (z - zd)/l_mod_rsl)
337 phim_zp =
stab_phi_mom(stabilitymethod, (zp - zd)/l_mod_rsl)
339 psim_hat_zp =
cal_psim_hat(stabilitymethod, zp, zh_rsl, l_mod_rsl, beta, lc)
341 xxm1 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
342 xxm1_2 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd + dz)/l_mod_rsl)
344 dphi = (xxm1_2 - xxm1)/dz
346 phi_hatmzh = kappa/(2.*beta*xxm1)
348 IF (phi_hatmzh > 1.)
THEN 353 c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa)
356 cm = (1.-phi_hatmzh)*exp(c2/2.)
359 psim_hat_z = psim_hat_zp + dz/2.*phim_zp*(cm*exp(-1.*c2*beta*(zp - zd)/elm))/(zp - zd)
360 psim_hat_z = psim_hat_z + dz/2.*phim_z*(cm*exp(-1.*c2*beta*(z - zd)/elm))/(z - zd)
364 function cal_psihatm_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm)
result(psihatm_z)
369 integer,
intent(in) :: StabilityMethod
370 integer,
intent(in) :: nz
371 real(KIND(1D0)),
DIMENSION(nz),
intent(in) :: zarray
372 real(KIND(1D0)),
intent(in) :: zh_RSL
373 real(KIND(1D0)),
intent(in) :: Lc
374 real(KIND(1D0)),
intent(in) :: beta
375 real(KIND(1D0)),
intent(in) :: L_MOD_RSL
376 real(KIND(1D0)),
intent(in) ::zd
377 real(KIND(1D0)),
intent(in) ::elm
380 real(KIND(1D0)),
DIMENSION(nz) ::psihatm_z
385 real(KIND(1D0)) ::phimz
386 real(KIND(1D0)) ::phimzp
388 real(KIND(1D0)) ::xxm1
389 real(KIND(1D0)) ::xxm1_2
390 real(KIND(1D0)) ::dphi
391 real(KIND(1D0)) ::phi_hatmZh
394 REAL(KIND(1d0)),
DIMENSION(nz):: dif
396 REAL(KIND(1d0)),
PARAMETER::kappa = 0.40
397 REAL(KIND(1d0)),
PARAMETER::dz = 0.1
401 psihatm_z = 0.*zarray
405 dif(z) = abs(zarray(z) - zh_rsl)
407 idx_can = minloc(dif, dim=1)
411 xxm1 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
412 xxm1_2 =
stab_phi_mom(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
414 phi_hatmzh = kappa/(2.*beta*xxm1)
417 IF (phi_hatmzh > 1.)
THEN 421 c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa)
424 cm = (1.-phi_hatmzh)*exp(c2/2.)
427 DO z = nz - 1, idx_can - 1, -1
428 phimz =
stab_phi_mom(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
429 phimzp =
stab_phi_mom(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
431 psihatm_z(z) = psihatm_z(z + 1) + dz/2.*phimzp*(cm*exp(-1.*c2*beta*(zarray(z + 1) - zd)/elm)) &
432 /(zarray(z + 1) - zd)
433 psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
440 function cal_psihath_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f)
result(psihath_z)
445 integer,
intent(in) :: StabilityMethod
446 integer,
intent(in) :: nz
448 real(KIND(1D0)),
DIMENSION(nz),
intent(in) :: zarray
449 real(KIND(1D0)),
intent(in) :: zh_RSL
450 real(KIND(1D0)),
intent(in) :: Lc
451 real(KIND(1D0)),
intent(in) :: beta
452 real(KIND(1D0)),
intent(in) :: Scc
453 real(KIND(1D0)),
intent(in) :: f
454 real(KIND(1D0)),
intent(in) :: L_MOD_RSL
455 real(KIND(1D0)),
intent(in) :: elm
456 real(KIND(1D0)),
intent(in) ::zd
459 real(KIND(1D0)),
DIMENSION(nz) ::psihath_z
464 real(KIND(1D0)) ::phihz
465 real(KIND(1D0)) ::phihzp
467 real(KIND(1D0)) ::xxh1
468 real(KIND(1D0)) ::xxh1_2
469 real(KIND(1D0)) ::dphih
470 real(KIND(1D0)) ::phi_hathZh
472 real(KIND(1D0)) ::c2h
473 REAL(KIND(1d0)),
DIMENSION(nz):: dif
475 REAL(KIND(1d0)),
PARAMETER::kappa = 0.40
476 REAL(KIND(1d0)),
PARAMETER::dz = 0.1
480 psihath_z = 0.*zarray
484 dif(z) = abs(zarray(z) - zh_rsl)
486 idx_can = minloc(dif, dim=1)
490 xxh1 =
stab_phi_heat(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
491 xxh1_2 =
stab_phi_heat(stabilitymethod, (zh_rsl - zd + 1.)/l_mod_rsl)
493 phi_hathzh = kappa*scc/(2.*beta*xxh1)
494 dphih = xxh1_2 - xxh1
496 IF (phi_hathzh > 1.)
THEN 501 c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
504 ch = (1.-phi_hathzh)*exp(c2h/2.)
506 DO z = nz - 1, idx_can - 1, -1
507 phihz =
stab_phi_heat(stabilitymethod, (zarray(z) - zd)/l_mod_rsl)
508 phihzp =
stab_phi_heat(stabilitymethod, (zarray(z + 1) - zd)/l_mod_rsl)
510 psihath_z(z) = psihath_z(z + 1) + dz/2.*phihzp*(ch*exp(-1.*c2h*beta*(zarray(z + 1) - zd)/elm)) &
511 /(zarray(z + 1) - zd)
512 psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
519 function rsl_cal_z0(StabilityMethod, zH_RSL, zd, beta, L_MOD_RSL, Lc)
result(z0)
523 integer,
intent(in) ::StabilityMethod
524 real(KIND(1D0)),
intent(in) :: zH_RSL
525 real(KIND(1D0)),
intent(in) :: zd
526 real(KIND(1D0)),
intent(in) :: L_MOD_RSL
527 real(KIND(1D0)),
intent(in) :: Lc
528 real(KIND(1D0)),
intent(in) :: beta
534 real(KIND(1D0)) ::psimZh, psimz0, z01, psihatm_Zh
535 real(KIND(1D0)) ::err
538 REAL(KIND(1d0)),
PARAMETER::kappa = 0.40
542 psimzh =
stab_psi_mom(stabilitymethod, (zh_rsl - zd)/l_mod_rsl)
543 psihatm_zh =
cal_psim_hat(stabilitymethod, zh_rsl, zh_rsl, l_mod_rsl, beta, lc)
552 DO WHILE ((err > 0.001) .AND. (it < 10))
555 z0 = (zh_rsl - zd)*exp(-1.*kappa/beta)*exp(-1.*psimzh + psimz0)*exp(psihatm_zh)
564 StabilityMethod, zh, L_MOD, sfr, planF, &!input
565 L_MOD_RSL, zH_RSL, Lc, beta, zd, z0, elm, Scc, f)
569 integer,
intent(in) :: StabilityMethod
570 real(KIND(1D0)),
intent(in) :: zh_min
571 real(KIND(1D0)),
intent(in) :: z0m
572 real(KIND(1D0)),
intent(in) :: zdm
573 real(KIND(1D0)),
intent(in) :: zh
574 real(KIND(1D0)),
intent(in) :: planF
575 real(KIND(1D0)),
intent(in) :: L_MOD
576 real(KIND(1D0)),
DIMENSION(nsurf),
intent(in) :: sfr
579 real(KIND(1D0)),
intent(out) ::L_MOD_RSL
580 real(KIND(1D0)),
intent(out) ::zH_RSL
581 real(KIND(1D0)),
intent(out) ::Lc
582 real(KIND(1D0)),
intent(out) ::beta
583 real(KIND(1D0)),
intent(out) ::zd
584 real(KIND(1D0)),
intent(out) ::z0
585 real(KIND(1D0)),
intent(out) ::elm
586 real(KIND(1D0)),
intent(out) ::Scc
587 real(KIND(1D0)),
intent(out) ::f
590 real(KIND(1D0)) ::sfr_zh
591 real(KIND(1D0)) ::sfr_tr
592 real(KIND(1D0)) ::phim
593 real(KIND(1D0)) ::betaN2
594 real(KIND(1D0)) ::betaHF
595 real(KIND(1D0)) ::betaNL
597 REAL(KIND(1d0)),
PARAMETER::planF_low = 1e-6
598 REAL(KIND(1d0)),
PARAMETER::kappa = 0.40
600 REAL(KIND(1d0)),
PARAMETER::r = 0.1
601 REAL(KIND(1d0)),
PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
605 l_mod_rsl = merge(l_mod, max(300., l_mod), l_mod < 0)
608 zh_rsl = max(zh, zh_min)
614 sfr_zh = min(sfr_zh, 0.8)
626 lc = (1.-sfr_zh)/planf*zh_rsl
631 scc = 0.5 + 0.3*tanh(2.*lc/l_mod_rsl)
632 f = 0.5*((1.+4.*r*scc)**0.5) - 0.5
640 betan2 = 0.30*sfr_tr/sfr_zh + (sfr_zh - sfr_tr)/sfr_zh*0.4
646 betanl = (kappa/2.)/phim
648 IF (lc/l_mod_rsl > a2)
THEN 651 beta = betanl + ((betahf - betanl)/(1.+a1*abs(lc/l_mod_rsl - a2)**a3))
657 zd = zh_rsl - (beta**2.)*lc
661 z0 =
rsl_cal_z0(stabilitymethod, zh_rsl, zd, beta, l_mod_rsl, lc)
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)
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)
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)
subroutine rslprofile(Zh, z0m, zdm, L_MOD, sfr, planF, 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 rsl_cal_z0(StabilityMethod, zH_RSL, zd, beta, L_MOD_RSL, Lc)
real(kind(1d0)) function rh2qa(RH, pres_hPa, Ta_degC)
recursive real(kind(1d0)) function cal_psim_hat(StabilityMethod, z, zh_RSL, L_MOD_RSL, beta, Lc)
real(kind(1d0)) function stab_phi_heat(StabilityMethod, ZL)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
real(kind(1d0)) function qa2rh(qa, pres_hPa, Ta_degC)
integer, parameter decidsurf
integer, parameter bldgsurf