SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines | Variables
rsl_module Module Reference

Functions/Subroutines

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)
 
recursive real(kind(1d0)) function cal_psim_hat (StabilityMethod, z, zh_RSL, L_MOD_RSL, beta, Lc)
 
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, dimension(nz) cal_psihath_z (StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f)
 
real(kind(1d0)) function rsl_cal_z0 (StabilityMethod, zH_RSL, zd, beta, L_MOD_RSL, Lc)
 
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)
 

Variables

integer, parameter nz = 30
 

Function/Subroutine Documentation

◆ cal_psihath_z()

real(kind(1d0)) function, dimension(nz) rsl_module::cal_psihath_z ( integer, intent(in)  StabilityMethod,
integer, intent(in)  nz,
real(kind(1d0)), dimension(nz), intent(in)  zarray,
real(kind(1d0)), intent(in)  L_MOD_RSL,
real(kind(1d0)), intent(in)  zH_RSL,
real(kind(1d0)), intent(in)  Lc,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  zd,
real(kind(1d0)), intent(in)  elm,
real(kind(1d0)), intent(in)  Scc,
real(kind(1d0)), intent(in)  f 
)

Definition at line 441 of file suews_phys_rslprof.f95.

References atmmoiststab_module::stab_phi_heat().

441 
442  ! calculate psi_hat for momentum
443  ! TS, 23 Oct 2019
444  implicit none
445  integer, intent(in) :: stabilitymethod ! stability method
446  integer, intent(in) :: nz ! number of vertical layers
447 
448  real(KIND(1D0)), DIMENSION(nz), intent(in) :: zarray ! height of interest [m]
449  real(KIND(1D0)), intent(in) :: zh_rsl ! canyon depth [m]
450  real(KIND(1D0)), intent(in) :: lc ! height scale for bluff bodies [m]
451  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
452  real(KIND(1D0)), intent(in) :: scc ! parameter in RSL
453  real(KIND(1D0)), intent(in) :: f ! parameter in RSL
454  real(KIND(1D0)), intent(in) :: l_mod_rsl ! Obukhov length [m]
455  real(KIND(1D0)), intent(in) :: elm ! displacement height used in RSL
456  real(KIND(1D0)), intent(in) ::zd ! displacement height used in RSL
457 
458  ! output
459  real(KIND(1D0)), DIMENSION(nz) ::psihath_z ! psim_hat at height of interest
460 
461  ! internal variables
462  ! real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
463  ! real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
464  real(KIND(1D0)) ::phihz ! displacement height used in RSL
465  real(KIND(1D0)) ::phihzp ! displacement height used in RSL
466  ! real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
467  real(KIND(1D0)) ::xxh1 ! displacement height used in RSL
468  real(KIND(1D0)) ::xxh1_2 ! displacement height used in RSL
469  real(KIND(1D0)) ::dphih ! displacement height used in RSL
470  real(KIND(1D0)) ::phi_hathzh ! displacement height used in RSL
471  real(KIND(1D0)) ::ch ! displacement height used in RSL
472  real(KIND(1D0)) ::c2h ! displacement height used in RSL
473  REAL(KIND(1d0)), DIMENSION(nz):: dif
474 
475  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
476  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
477 
478  integer::z, idx_can
479 
480  psihath_z = 0.*zarray
481 
482  ! determine index at the canyon top
483  DO z = 1, nz
484  dif(z) = abs(zarray(z) - zh_rsl)
485  ENDDO
486  idx_can = minloc(dif, dim=1)
487  ! zarray(idx_can) = Zh_RSL
488 
489  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
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)
492 
493  phi_hathzh = kappa*scc/(2.*beta*xxh1)
494  dphih = xxh1_2 - xxh1
495 
496  IF (phi_hathzh > 1.) THEN
497  ! c2 = 0.5 ! more stable, but less correct
498  c2h = 0.5
499  ELSE
500  ! 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
501  c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
502  ENDIF
503  ! cm = (1.-phi_hatmZh)*EXP(c2/2.)
504  ch = (1.-phi_hathzh)*exp(c2h/2.)
505 
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)
509 
510  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
511  /(zarray(z + 1) - zd)
512  psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
513  /(zarray(z) - zd)
514 
515  ENDDO
516 
real(kind(1d0)) z
Here is the call graph for this function:

◆ cal_psihatm_z()

real(kind(1d0)) function, dimension(nz) rsl_module::cal_psihatm_z ( integer, intent(in)  StabilityMethod,
integer, intent(in)  nz,
real(kind(1d0)), dimension(nz), intent(in)  zarray,
real(kind(1d0)), intent(in)  L_MOD_RSL,
real(kind(1d0)), intent(in)  zH_RSL,
real(kind(1d0)), intent(in)  Lc,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  zd,
real(kind(1d0)), intent(in)  elm 
)

Definition at line 365 of file suews_phys_rslprof.f95.

References atmmoiststab_module::stab_phi_mom().

365 
366  ! calculate psi_hat for momentum
367  ! TS, 23 Oct 2019
368  implicit none
369  integer, intent(in) :: stabilitymethod ! stability method
370  integer, intent(in) :: nz ! number of vertical layers
371  real(KIND(1D0)), DIMENSION(nz), intent(in) :: zarray ! height of interest [m]
372  real(KIND(1D0)), intent(in) :: zh_rsl ! canyon depth [m]
373  real(KIND(1D0)), intent(in) :: lc ! height scale for bluff bodies [m]
374  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
375  real(KIND(1D0)), intent(in) :: l_mod_rsl ! Obukhov length [m]
376  real(KIND(1D0)), intent(in) ::zd ! displacement height used in RSL
377  real(KIND(1D0)), intent(in) ::elm ! displacement height used in RSL
378 
379  ! output
380  real(KIND(1D0)), DIMENSION(nz) ::psihatm_z ! psim_hat at height of interest
381 
382  ! internal variables
383  ! real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
384  ! real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
385  real(KIND(1D0)) ::phimz ! displacement height used in RSL
386  real(KIND(1D0)) ::phimzp ! displacement height used in RSL
387  ! real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
388  real(KIND(1D0)) ::xxm1 ! displacement height used in RSL
389  real(KIND(1D0)) ::xxm1_2 ! displacement height used in RSL
390  real(KIND(1D0)) ::dphi ! displacement height used in RSL
391  real(KIND(1D0)) ::phi_hatmzh ! displacement height used in RSL
392  real(KIND(1D0)) ::cm ! displacement height used in RSL
393  real(KIND(1D0)) ::c2 ! displacement height used in RSL
394  REAL(KIND(1d0)), DIMENSION(nz):: dif
395 
396  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
397  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
398 
399  integer::z, idx_can
400 
401  psihatm_z = 0.*zarray
402 
403  ! determine index at the canyon top
404  DO z = 1, nz
405  dif(z) = abs(zarray(z) - zh_rsl)
406  ENDDO
407  idx_can = minloc(dif, dim=1)
408  ! zarray(idx_can) = Zh_RSL
409 
410  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
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)
413 
414  phi_hatmzh = kappa/(2.*beta*xxm1)
415  dphi = xxm1_2 - xxm1
416 
417  IF (phi_hatmzh > 1.) THEN
418  c2 = 0.5 ! more stable, but less correct
419  ! c2h = 0.5
420  ELSE
421  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
422  ! c2h = (kappa*Scc*(2.+f - (dphih*2.*beta**2.*Lc/xxh1)))/(2.*beta*xxh1 - kappa*Scc)
423  ENDIF
424  cm = (1.-phi_hatmzh)*exp(c2/2.)
425  ! ch = (1.-phi_hathzh)*EXP(c2h/2.)
426 
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)
430 
431  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
432  /(zarray(z + 1) - zd)
433  psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
434  /(zarray(z) - zd)
435 
436  ENDDO
437 
real(kind(1d0)) z
real(kind(1d0)) cm
Here is the call graph for this function:

◆ cal_psim_hat()

recursive real(kind(1d0)) function rsl_module::cal_psim_hat ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  z,
real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  L_MOD_RSL,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 292 of file suews_phys_rslprof.f95.

References atmmoiststab_module::stab_phi_mom().

Referenced by rsl_cal_z0().

292  ! calculate psi_hat for momentum
293  ! TS, 23 Oct 2019
294  implicit none
295  integer, intent(in) :: stabilitymethod ! stability method
296  real(KIND(1D0)), intent(in) :: z ! height of interest [m]
297  real(KIND(1D0)), intent(in) :: zh_rsl ! canyon depth [m]
298  real(KIND(1D0)), intent(in) :: lc ! height scale for bluff bodies [m]
299  real(KIND(1D0)), intent(in) :: beta ! parameter in RSL
300  real(KIND(1D0)), intent(in) :: l_mod_rsl ! Obukhov length [m]
301 
302  ! output
303  real(KIND(1D0)) ::psim_hat_z ! psim_hat at height of interest
304 
305  ! internal variables
306  real(KIND(1D0)) ::zp ! a height above z used for iterative calculations
307  real(KIND(1D0)) ::zd ! displacement height used in RSL
308  real(KIND(1D0)) ::phim_lc ! displacement height used in RSL
309  real(KIND(1D0)) ::phim_z ! displacement height used in RSL
310  real(KIND(1D0)) ::phim_zp ! displacement height used in RSL
311  real(KIND(1D0)) ::psim_hat_zp ! displacement height used in RSL
312  real(KIND(1D0)) ::elm ! displacement height used in RSL
313  real(KIND(1D0)) ::xxm1 ! displacement height used in RSL
314  real(KIND(1D0)) ::xxm1_2 ! displacement height used in RSL
315  real(KIND(1D0)) ::dphi ! displacement height used in RSL
316  real(KIND(1D0)) ::phi_hatmzh ! displacement height used in RSL
317  real(KIND(1D0)) ::cm ! displacement height used in RSL
318  real(KIND(1D0)) ::c2 ! displacement height used in RSL
319 
320  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
321  REAL(KIND(1d0)), PARAMETER::dz = 0.1 !height step
322 
323  if (z > 100) then
324  psim_hat_z = 0.
325  return
326  end if
327 
328  zp = 1.01*z ! a height above z
329 
330  zd = zh_rsl - (beta**2.)*lc
331  elm = 2.*beta**3*lc
332 
333  ! phim at Lc
334  phim_lc = stab_phi_mom(stabilitymethod, lc/l_mod_rsl)
335 
336  phim_z = stab_phi_mom(stabilitymethod, (z - zd)/l_mod_rsl)
337  phim_zp = stab_phi_mom(stabilitymethod, (zp - zd)/l_mod_rsl)
338 
339  psim_hat_zp = cal_psim_hat(stabilitymethod, zp, zh_rsl, l_mod_rsl, beta, lc)
340 
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)
343 
344  dphi = (xxm1_2 - xxm1)/dz
345 
346  phi_hatmzh = kappa/(2.*beta*xxm1)
347 
348  IF (phi_hatmzh > 1.) THEN
349  ! more stable, but less correct
350  c2 = 0.5
351  ELSE
352  ! if very unstable this might cause some high values of psihat_z
353  c2 = (kappa*(3.-(2.*beta**2.*lc/xxm1*dphi)))/(2.*beta*xxm1 - kappa)
354  ENDIF
355 
356  cm = (1.-phi_hatmzh)*exp(c2/2.)
357 
358  !Taylor's approximation for integral
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)
361 
real(kind(1d0)) z
real(kind(1d0)) cm
Here is the call graph for this function:
Here is the caller graph for this function:

◆ rsl_cal_prms()

subroutine rsl_module::rsl_cal_prms ( real(kind(1d0)), intent(in)  zh_min,
real(kind(1d0)), intent(in)  z0m,
real(kind(1d0)), intent(in)  zdm,
integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  zh,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), intent(in)  planF,
real(kind(1d0)), intent(out)  L_MOD_RSL,
real(kind(1d0)), intent(out)  zH_RSL,
real(kind(1d0)), intent(out)  Lc,
real(kind(1d0)), intent(out)  beta,
real(kind(1d0)), intent(out)  zd,
real(kind(1d0)), intent(out)  z0,
real(kind(1d0)), intent(out)  elm,
real(kind(1d0)), intent(out)  Scc,
real(kind(1d0)), intent(out)  f 
)

Definition at line 566 of file suews_phys_rslprof.f95.

References allocatearray::bldgsurf, allocatearray::conifsurf, allocatearray::decidsurf, rsl_cal_z0(), and atmmoiststab_module::stab_phi_mom().

Referenced by rslprofile().

566  ! calculate surface/skin temperature
567  ! TS, 23 Oct 2019
568  implicit none
569  integer, intent(in) :: stabilitymethod ! stability method
570  real(KIND(1D0)), intent(in) :: zh_min ! canyon depth [m]
571  real(KIND(1D0)), intent(in) :: z0m ! roughness length to prescribe if necessary [m]
572  real(KIND(1D0)), intent(in) :: zdm ! displacement height to prescribe if necessary [m]
573  real(KIND(1D0)), intent(in) :: zh ! canyon depth [m]
574  real(KIND(1D0)), intent(in) :: planf ! frontal area index
575  real(KIND(1D0)), intent(in) :: l_mod ! Obukhov length [m]
576  real(KIND(1D0)), DIMENSION(nsurf), intent(in) :: sfr ! land cover fractions
577 
578  ! output
579  real(KIND(1D0)), intent(out) ::l_mod_rsl ! Obukhov length used in RSL module with thresholds applied
580  real(KIND(1D0)), intent(out) ::zh_rsl ! mean canyon height used in RSL module with thresholds applied
581  real(KIND(1D0)), intent(out) ::lc ! height scale for bluff bodies [m]
582  real(KIND(1D0)), intent(out) ::beta ! psim_hat at height of interest
583  real(KIND(1D0)), intent(out) ::zd ! displacement height to prescribe if necessary [m]
584  real(KIND(1D0)), intent(out) ::z0 ! roughness length [m]
585  real(KIND(1D0)), intent(out) ::elm ! length scale used in RSL
586  real(KIND(1D0)), intent(out) ::scc ! parameter in RSL
587  real(KIND(1D0)), intent(out) ::f ! parameter in RSL
588 
589  ! internal variables
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
596 
597  REAL(KIND(1d0)), PARAMETER::planf_low = 1e-6
598  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
599  ! REAL(KIND(1d0)), PARAMETER::z0m= 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.
602  ! REAL(KIND(1d0)), parameter::Zh_min = 0.4! limit for minimum canyon height used in RSL module
603 
604  ! under stable conditions, set a threshold for L_MOD to avoid numerical issues. TS 28 Oct 2019
605  l_mod_rsl = merge(l_mod, max(300., l_mod), l_mod < 0)
606 
607  ! zH_RSL
608  zh_rsl = max(zh, zh_min)
609  ! zH_RSL = zh
610 
611  ! land cover fraction of bluff bodies
612  sfr_zh = sum(sfr([bldgsurf, conifsurf, decidsurf]))
613  ! set a threshold for sfr_zh to avoid numerical difficulties
614  sfr_zh = min(sfr_zh, 0.8)
615 
616  ! land cover fraction of trees
617  sfr_tr = sum(sfr([conifsurf, decidsurf]))
618 
619  ! height scale for buildings !not used? why?
620  ! Lc_build = (1.-sfr(BldgSurf))/planF*Zh_RSL ! Coceal and Belcher 2004 assuming Cd = 2
621 
622  ! height scale for tress
623  ! Lc_tree = 1./(cd_tree*a_tree) ! not used? why?
624 
625  ! height scale for bluff bodies
626  lc = (1.-sfr_zh)/planf*zh_rsl
627 
628  ! phim at Lc
629  phim = stab_phi_mom(stabilitymethod, lc/l_mod_rsl)
630 
631  scc = 0.5 + 0.3*tanh(2.*lc/l_mod_rsl) ! Schmidt number Harman and Finnigan 2008: assuming the same for heat and momemntum
632  f = 0.5*((1.+4.*r*scc)**0.5) - 0.5
633 
634  !
635  ! Step 2:
636  ! Parameterise beta according to Harman 2012 with upper limit of 0.5
637  ! betaN for trees found to be 0.3 and for urban 0.4 linearly interpolate between the two using surface fractions
638  ! betaN2 = 0.30 + (1.-sfr(ConifSurf) - sfr(ConifSurf))*0.1
639  if (sfr_zh > 0) then
640  betan2 = 0.30*sfr_tr/sfr_zh + (sfr_zh - sfr_tr)/sfr_zh*0.4
641  ELSE
642  betan2 = 0.35
643  endif
644 
645  betahf = betan2/phim
646  betanl = (kappa/2.)/phim
647 
648  IF (lc/l_mod_rsl > a2) THEN
649  beta = betahf
650  ELSE
651  beta = betanl + ((betahf - betanl)/(1.+a1*abs(lc/l_mod_rsl - a2)**a3))
652  ENDIF
653 
654  IF (beta > 0.5) THEN
655  beta = 0.5
656  ENDIF
657  zd = zh_rsl - (beta**2.)*lc
658  elm = 2.*beta**3*lc
659 
660  ! calculate z0 iteratively
661  z0 = rsl_cal_z0(stabilitymethod, zh_rsl, zd, beta, l_mod_rsl, lc)
662 
663  ! correct parameters if RSL approach doesn't apply for a shallow canyon
664  ! if (zh_RSL - zd < z0) then
665  ! ! when zh_RSL is too shallow, implying RSL doesn't apply, force RSL correction to zero
666  ! ! psihatm_z=0
667  ! ! psihath_z=0
668  ! beta = 1.e6
669  ! !correct RSL-based z0 and zd using Rule of thumb (G&O 1999)
670  ! zd = 0.7*zH_RSL
671  ! z0 = 0.1*zH_RSL
672  ! zd = zdm
673  ! z0 = z0m
674  ! ! then MOST recovers from RSL correction
675  ! end if
676 
real(kind(1d0)) zdm
real(kind(1d0)) z0m
Here is the call graph for this function:
Here is the caller graph for this function:

◆ rsl_cal_z0()

real(kind(1d0)) function rsl_module::rsl_cal_z0 ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  zH_RSL,
real(kind(1d0)), intent(in)  zd,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  L_MOD_RSL,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 520 of file suews_phys_rslprof.f95.

References cal_psim_hat(), and atmmoiststab_module::stab_psi_mom().

Referenced by rsl_cal_prms().

520  ! calculate z0 iteratively
521  ! TS, 23 Oct 2019
522  implicit none
523  integer, intent(in) ::stabilitymethod
524  real(KIND(1D0)), intent(in) :: zh_rsl ! canyon depth [m]
525  real(KIND(1D0)), intent(in) :: zd ! displacement height [m]
526  real(KIND(1D0)), intent(in) :: l_mod_rsl ! Monin Obukhov length[m]
527  real(KIND(1D0)), intent(in) :: lc ! canyon length scale [m]
528  real(KIND(1D0)), intent(in) :: beta ! height scale for bluff bodies [m]
529 
530  ! output
531  real(KIND(1D0)) ::z0
532 
533  ! internal variables
534  real(KIND(1D0)) ::psimzh, psimz0, z01, psihatm_zh
535  real(KIND(1D0)) ::err
536  integer ::it
537 
538  REAL(KIND(1d0)), PARAMETER::kappa = 0.40
539  ! REAL(KIND(1d0)), PARAMETER::r = 0.1
540  ! REAL(KIND(1d0)), PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
541 
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)
544 
545  !first guess
546  z0 = 0.5
547  ! if (Zh_RSL>Zh_min) then
548  ! if Zh>Zh_min, calculate z0 using RSL method
549  err = 10.
550  ! psimz0 = 0.5
551  it = 1
552  DO WHILE ((err > 0.001) .AND. (it < 10))
553  psimz0 = stab_psi_mom(stabilitymethod, z0/l_mod_rsl)
554  z01 = z0
555  z0 = (zh_rsl - zd)*exp(-1.*kappa/beta)*exp(-1.*psimzh + psimz0)*exp(psihatm_zh)
556  err = abs(z01 - z0)
557  it = it + 1
558  ENDDO
559 
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
integer it
Here is the call graph for this function:
Here is the caller graph for this function:

◆ rslprofile()

subroutine rsl_module::rslprofile ( real(kind(1d0)), intent(in)  Zh,
real(kind(1d0)), intent(in)  z0m,
real(kind(1d0)), intent(in)  zdm,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
real(kind(1d0)), intent(in)  planF,
integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  avcp,
real(kind(1d0)), intent(in)  lv_J_kg,
real(kind(1d0)), intent(in)  avdens,
real(kind(1d0)), intent(in)  avU1,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  avRH,
real(kind(1d0)), intent(in)  Press_hPa,
real(kind(1d0)), intent(in)  zMeas,
real(kind(1d0)), intent(in)  qh,
real(kind(1d0)), intent(in)  qe,
real(kind(1d0)), intent(out)  T2_C,
real(kind(1d0)), intent(out)  q2_gkg,
real(kind(1d0)), intent(out)  U10_ms,
real(kind(1d0)), intent(out)  RH2,
real(kind(1d0)), dimension(nz*4), intent(out)  dataoutLineRSL 
)

Definition at line 19 of file suews_phys_rslprof.f95.

References meteo::qa2rh(), meteo::rh2qa(), rsl_cal_prms(), atmmoiststab_module::stab_phi_heat(), atmmoiststab_module::stab_phi_mom(), atmmoiststab_module::stab_psi_heat(), and atmmoiststab_module::stab_psi_mom().

Referenced by suews_driver::suews_cal_main().

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  !
27  !-----------------------------------------------------
28 
29  IMPLICIT NONE
30 
31  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::sfr! surface fractions [-]
32  REAL(KIND(1d0)), INTENT(in):: zmeas ! height of atmospheric forcing [m]
33  REAL(KIND(1d0)), INTENT(in):: avu1 ! Wind speed at forcing height [m s-1]
34  REAL(KIND(1d0)), INTENT(in):: temp_c ! Air temperature at forcing height [C]
35  REAL(KIND(1d0)), INTENT(in):: avrh ! relative humidity at forcing height [-]
36  REAL(KIND(1d0)), INTENT(in):: press_hpa ! pressure at forcing height [hPa]
37  REAL(KIND(1d0)), INTENT(in):: l_mod ! Obukhov length [m]
38  REAL(KIND(1d0)), INTENT(in):: avcp ! specific heat capacity [J kg-1 K-1]
39  REAL(KIND(1d0)), INTENT(in):: lv_j_kg ! Latent heat of vaporization in [J kg-1]
40  REAL(KIND(1d0)), INTENT(in):: avdens ! air density [kg m-3]
41  REAL(KIND(1d0)), INTENT(in):: qh ! sensible heat flux [W m-2]
42  REAL(KIND(1d0)), INTENT(in):: qe ! Latent heat flux [W m-2]
43  REAL(KIND(1d0)), INTENT(in):: zh ! Mean building height [m]
44  REAL(KIND(1d0)), INTENT(in):: z0m ! Mean building height [m]
45  REAL(KIND(1d0)), INTENT(in):: zdm ! Mean building height [m]
46  REAL(KIND(1d0)), INTENT(in):: planf ! Frontal area index [-]
47 
48  INTEGER, INTENT(in)::stabilitymethod
49 
50  REAL(KIND(1d0)), INTENT(out):: t2_c ! Air temperature at 2 m [C]
51  REAL(KIND(1d0)), INTENT(out):: q2_gkg ! Air specific humidity at 2 m [g kg-1]
52  REAL(KIND(1d0)), INTENT(out):: u10_ms ! wind speed at 10 m [m s-1]
53  REAL(KIND(1d0)), INTENT(out):: rh2 ! Air relative humidity [-]
54 
55  INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
56 
57  REAL(KIND(1d0)), PARAMETER:: cd_tree = 1.2, & ! drag coefficient tree canopy !!!!needs adjusting!!!
58  a_tree = 0.05, & ! the foliage area per unit volume !!!!needs adjusting!!!
59  kappa = 0.40, &! von karman constant
60  ! lv_J_kg = 2.5E6, &! latent heat for water vapor!!! make consistant with rest of code
61  beta_n = 0.40, & ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
62  pi = 4.*atan(1.0), r = 0.1, &
63  a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
64 
65  ! REAL(KIND(1d0)), INTENT(out), DIMENSION(nz*3):: zarrays ! Height array
66  REAL(KIND(1d0)), INTENT(out), DIMENSION(nz*4):: dataoutlinersl ! Variables array (/U,T,q/)
67  REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z
68  REAL(KIND(1d0)), DIMENSION(nz):: psihath_z
69  REAL(KIND(1d0)), DIMENSION(nz):: dif
70  ! REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z, psihath_z
71  REAL(KIND(1d0)), DIMENSION(nz):: zarray
72  REAL(KIND(1d0)), DIMENSION(nz):: dataoutlineursl ! wind speed array [m s-1]
73  REAL(KIND(1d0)), DIMENSION(nz):: dataoutlinetrsl ! Temperature array [C]
74  REAL(KIND(1d0)), DIMENSION(nz):: dataoutlineqrsl ! Specific humidity array [g kg-1]
75 
76  REAL(KIND(1d0))::z0 ! roughness length from H&F
77  REAL(KIND(1d0))::zd ! zero-plane displacement
78 
79  REAL(KIND(1d0))::lc_build, lc_tree, lc ! canopy drag length scale
80  REAL(KIND(1d0))::scc ! Schmidt number for temperature and humidity
81  REAL(KIND(1d0))::phim, psimz, psimzh, psimz0, psimza, phi_hatmzh, phi_hathzh, phimzp, phimz, phihzp, phihz, psihz, psihza ! stability function for momentum
82  REAL(KIND(1d0))::betahf, betanl, beta, betan2 ! beta coefficient from Harman 2012
83  REAL(KIND(1d0))::elm ! mixing length
84  REAL(KIND(1d0))::xxm1, xxm1_2, xxh1, xxh1_2, err, z01, dphi, dphih ! dummy variables for stability functions
85  REAL(KIND(1d0))::f, cm, c2, ch, c2h ! H&F'07 and H&F'08 'constants'
86  REAL(KIND(1d0))::t_h, q_h ! H&F'08 canopy corrections
87  REAL(KIND(1d0))::tstar_rsl ! temperature scale
88  REAL(KIND(1d0))::ustar_rsl ! friction velocity used in RSL
89  REAL(KIND(1d0))::sfr_zh ! land cover fraction of bluff bodies: buildings and trees
90  REAL(KIND(1d0))::sfr_tr ! land cover fraction of trees
91  REAL(KIND(1d0))::l_mod_rsl ! Obukhov length used in RSL module with thresholds applied
92  REAL(KIND(1d0))::zh_rsl ! mean canyon height used in RSL module with thresholds applied
93  REAL(KIND(1d0))::dz! initial height step
94  REAL(KIND(1d0)), parameter::zh_min = 0.15! limit for minimum canyon height used in RSL module
95  REAL(KIND(1d0)), parameter::ratio_dz = 1.618! ratio between neighbouring height steps
96 
97  REAL(KIND(1d0))::qa_gkg, qstar_rsl ! specific humidity scale
98  INTEGER :: i, z, it, idx_can, idx_za, idx_2m, idx_10m
99  INTEGER :: nz_can ! number of heights in canyon
100  !
101  ! Step 1: Calculate grid-cell dependent constants
102  ! Step 2: Calculate Beta (crucial for H&F method)
103  ! Step 3: calculate the stability dependent H&F constants
104  ! Step 4: determine psihat at levels above the canopy
105  ! Step 5: Calculate z0 iteratively
106  ! Step 6: Calculate mean variables above canopy
107  ! Step 7: Calculate mean variables in canopy
108  !
109  ! ! Step 1
110  ! ! Start setting up the parameters
111 
112  call rsl_cal_prms( &
113  zh_min, &
114  z0m, zdm, &
115  stabilitymethod, zh, l_mod, sfr, planf, &!input
116  l_mod_rsl, zh_rsl, lc, beta, zd, z0, elm, scc, f)
117 
118  ! Define the height array with consideration of key heights
119  ! set number of heights within canopy
120  IF (zh_rsl <= 2) THEN
121  nz_can = 5
122  ELSE IF (zh_rsl <= 10) THEN
123  nz_can = 10
124  else
125  nz_can = 15
126  ENDIF
127  ! fill up heights in canopy
128  dz = zh_rsl/nz_can
129  do i = 1, nz_can
130  zarray(i) = dz*i
131  end do
132  ! fill up heights above canopy
133  dz = (zmeas - zh_rsl)/(nz - nz_can)
134  do i = nz_can + 1, nz
135  zarray(i) = zh_rsl + (i - nz_can)*dz
136  end do
137 
138  ! add key heights (2m and 10m) to zarray
139  ! 2m:
140  DO z = 1, nz
141  dif(z) = abs(zarray(z) - 2)
142  ENDDO
143  idx_2m = minloc(dif, dim=1)
144  zarray(idx_2m) = 2
145  ! 10m:
146  DO z = 1, nz
147  dif(z) = abs(zarray(z) - 10)
148  ENDDO
149  idx_10m = minloc(dif, dim=1)
150  zarray(idx_10m) = 10
151 
152  ! determine index at the canyon top
153  DO z = 1, nz
154  dif(z) = abs(zarray(z) - zh_rsl)
155  ENDDO
156  idx_can = minloc(dif, dim=1)
157  zarray(idx_can) = zh_rsl
158 
159  ! determine index at measurement height
160  DO z = 1, nz
161  ! dif2(z) = ABS(zarray(z) - (zMeas - zd))
162  dif(z) = abs(zarray(z) - zmeas)
163  ENDDO
164  idx_za = minloc(dif, dim=1)
165  zarray(idx_za) = zmeas
166 
167  if (zh_rsl - zd < z0 .or. zh < zh_min) then
168  ! correct parameters if RSL approach doesn't apply for a shallow canyon
169  ! when zh_RSL is too shallow, implying RSL doesn't apply, force RSL correction to zero
170  psihatm_z = 0
171  psihath_z = 0
172  beta = 1.e6
173  !correct RSL-based z0 and zd using Rule of thumb (G&O 1999)
174  ! zd = 0.7*zH_RSL
175  ! z0 = 0.1*zH_RSL
176  !correct RSL-based using SUEWS system-wide values
177  zd = zdm
178  z0 = z0m
179  ! then MOST recovers from RSL correction
180  else
181  !otherwise use RSL approach to calculate correction factors
182  ! Step 3:
183  !
184  ! calculate phihatM according to H&F '07 and H&F '08 for heat and humidity
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)
187 
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)
190 
191  phi_hatmzh = kappa/(2.*beta*xxm1)
192  phi_hathzh = kappa*scc/(2.*beta*xxh1)
193  dphi = xxm1_2 - xxm1
194  dphih = xxh1_2 - xxh1
195  IF (phi_hatmzh > 1.) THEN
196  c2 = 0.5 ! more stable, but less correct
197  c2h = 0.5
198  ELSE
199  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
200  c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/xxh1)))/(2.*beta*xxh1 - kappa*scc)
201  ENDIF
202  cm = (1.-phi_hatmzh)*exp(c2/2.)
203  ch = (1.-phi_hathzh)*exp(c2h/2.)
204 
205  !
206  ! Step 4: determine psihat at levels above the canopy
207  !
208  psihatm_z = 0.*zarray
209  psihath_z = 0.*zarray
210  ! psihatm_zp = 0.
211  ! psihath_zp = 0
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)
217 
218  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
219  /(zarray(z + 1) - zd)
220  psihatm_z(z) = psihatm_z(z) + dz/2.*phimz*(cm*exp(-1.*c2*beta*(zarray(z) - zd)/elm)) &
221  /(zarray(z) - zd)
222  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
223  /(zarray(z + 1) - zd)
224  psihath_z(z) = psihath_z(z) + dz/2.*phihz*(ch*exp(-1.*c2h*beta*(zarray(z) - zd)/elm)) &
225  /(zarray(z) - zd)
226  ENDDO
227  ! psihatm_z=cal_psihatm_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm)
228  ! psihath_z=cal_psihath_z(StabilityMethod, nz, zarray, L_MOD_RSL, zH_RSL, Lc, beta, zd, elm, Scc, f)
229 
230  end if
231 
232  ! calculate above canopy variables
233  !
234  psimz0 = stab_psi_mom(stabilitymethod, z0/l_mod_rsl)
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)
241 
242  DO z = idx_can, nz
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 ! this is different from Theeuwes et al. (2019 BLM)
246  dataoutlineursl(z) = (log((zarray(z) - zd)/z0) - psimz + psimz0 + psihatm_z(z))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
247  ! 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)
248  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)
249  dataoutlineqrsl(z) = (log((zarray(z) - zd)/(zmeas - zd)) - psihz + psihza + psihath_z(z) - psihath_z(idx_za))/kappa
250  ENDDO
251  !
252  ! Step 7
253  ! calculate in canopy variables
254  !
255  t_h = scc*tstar_rsl/(beta*f)
256  q_h = scc*qstar_rsl/(beta*f)
257  DO z = 1, idx_can
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
261  ENDDO
262 
263  dataoutlineursl = dataoutlineursl*ustar_rsl
264  dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
265  dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
266 
267  dataoutlinersl = (/zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl/)
268 
269  !
270  ! Step 8
271  ! retrieve the diagnostics at key heights
272  !
273  t2_c = dataoutlinetrsl(idx_2m)
274  q2_gkg = dataoutlineqrsl(idx_2m)
275  u10_ms = dataoutlineursl(idx_10m)
276  ! get relative humidity:
277  rh2 = qa2rh(q2_gkg, press_hpa, t2_c)
278 
279 ! print *, 'Wind speed', dataoutLineURSL
280  ! DO z = 1, nz
281  ! print *, dataoutLineTRSL(z)
282  ! ENDDO
283  ! DO z = 1, nz
284  ! print *, zarray(z)
285  ! ENDDO
286 ! print *, 'Temperature', Temp_C, dataoutLineTRSL
287 ! print *, 'qStar', qStar, qe
288 ! print *, 'humidity' , qa_gkg, dataoutLineqRSL*1000.
real(kind(1d0)) z
real(kind(1d0)) zdm
real(kind(1d0)) lv_j_kg
real(kind(1d0)) cm
real(kind(1d0)) avcp
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
integer it
real(kind(1d0)) z0m
real(kind(1d0)), parameter pi
real(kind(1d0)) avdens
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ nz

integer, parameter rsl_module::nz = 30

Definition at line 8 of file suews_phys_rslprof.f95.

8  INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy