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, FAI, 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 442 of file suews_phys_rslprof.f95.

References atmmoiststab_module::stab_phi_heat().

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

References atmmoiststab_module::stab_phi_mom().

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

References atmmoiststab_module::stab_phi_mom().

Referenced by rsl_cal_z0().

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

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

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

Referenced by rsl_cal_prms().

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