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

Functions/Subroutines

subroutine rslprofile (DiagMethod, Zh, z0m, zdm, z0v, L_MOD, sfr_surf, FAI, PAI, StabilityMethod, RA_h, 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 interp_z (z_x, z, v)
 
real(kind(1d0)) function cal_elm_rsl (beta, Lc)
 
recursive real(kind(1d0)) function cal_psim_hat (StabilityMethod, psihatm_top, psihatm_mid, z_top, z_mid, z_btm, cm, c2, zh_RSL, zd_RSL, L_MOD, beta, elm, Lc)
 
recursive real(kind(1d0)) function cal_psih_hat (StabilityMethod, psihath_top, psihath_mid, z_top, z_mid, z_btm, ch, c2h, zh_RSL, zd_RSL, L_MOD, beta, elm, Lc)
 
real(kind(1d0)) function cal_phim_hat (StabilityMethod, z, zh_RSL, L_MOD, beta, lc)
 
subroutine cal_cm (StabilityMethod, zh_RSL, zd_RSL, Lc, beta, L_MOD, c2, cm)
 
subroutine cal_ch (StabilityMethod, zh_RSL, zd_RSL, Lc, beta, L_MOD, Scc, f, c2h, ch)
 
real(kind(1d0)) function cal_zd_rsl (zh_RSL, beta, Lc)
 
real(kind(1d0)) function cal_z0_rsl (StabilityMethod, zH_RSL, zd_RSL, beta, L_MOD_RSL, psihatm_Zh)
 
subroutine rsl_cal_prms (StabilityMethod, nz_above, z_array, zh, L_MOD, sfr_surf, FAI, PAI, psihatm_array, psihath_array, zH_RSL, L_MOD_RSL, Lc, beta, zd_RSL, z0_RSL, elm, Scc, fx)
 
real(kind(1d0)) function cal_beta_rsl (StabilityMethod, PAI, sfr_tr, lc_over_L)
 
real(kind(1d0)) function cal_beta_lc (stabilityMethod, beta0, lc_over_l)
 

Variables

integer, parameter nz = 30
 

Function/Subroutine Documentation

◆ cal_beta_lc()

real(kind(1d0)) function rsl_module::cal_beta_lc ( integer, intent(in)  stabilityMethod,
real(kind(1d0)), intent(in)  beta0,
real(kind(1d0)), intent(in)  lc_over_l 
)

Definition at line 1141 of file suews_phys_rslprof.f95.

1142 ! TS, 03 Aug 2020:
1143 ! iterative determination of beta depending on Lc/L
1144 ! ref: eqn 10 & 11 in Harman (2012, BLM)
1145 IMPLICIT NONE
1146 INTEGER, INTENT(in) :: StabilityMethod
1147 REAL(KIND(1D0)), INTENT(in) :: beta0
1148 REAL(KIND(1D0)), INTENT(in) :: lc_over_l
1149 REAL(KIND(1D0)) :: beta_x
1150
1151 REAL(KIND(1D0)) :: phim, err, beta_x0
1152
1153 INTEGER :: it
1154
1155 it = 1
1156 phim = 1
1157 err = 1
1158 ! print *, '***********************'
1159 ! print *, 'beta0:', beta0
1160 ! print *, 'Lc/L_MOD:', lc_over_l
1161 DO WHILE ((err > 0.001) .AND. (it < 20))
1162 beta_x0 = beta0/phim
1163 phim = stab_phi_heat(stabilitymethod, (beta_x0**2)*lc_over_l)
1164 ! TODO: how to deal with neutral condition when phim=0? TS 05 Feb 2021
1165 ! here we use beta=0.35 as a temporary workaround but a better solution is still neded.
1166 beta_x = beta0/phim
1167 err = abs(beta_x - beta_x0)
1168 ! print *, it, err, beta_x0, beta_x, phim, lc_over_l
1169 it = it + 1
1170
1171 END DO
1172 ! print *, ''
1173

References atmmoiststab_module::stab_phi_heat().

Referenced by cal_beta_rsl().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_beta_rsl()

real(kind(1d0)) function rsl_module::cal_beta_rsl ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  PAI,
real(kind(1d0)), intent(in)  sfr_tr,
real(kind(1d0)), intent(in)  lc_over_L 
)

Definition at line 1090 of file suews_phys_rslprof.f95.

1091 ! Step 2:
1092 ! Parameterise beta according to Harman 2012 with upper limit of 0.5
1093 IMPLICIT NONE
1094
1095 INTEGER, INTENT(in) :: StabilityMethod ! stability method
1096 REAL(KIND(1D0)), INTENT(in) :: PAI
1097 REAL(KIND(1D0)), INTENT(in) :: sfr_tr
1098 REAL(KIND(1D0)), INTENT(in) :: lc_over_L
1099
1100 ! output
1101 REAL(KIND(1D0)) :: beta
1102
1103 ! internal use
1104 REAL(KIND(1D0)) :: betaHF
1105 REAL(KIND(1D0)) :: betaNL
1106
1107 REAL(KIND(1D0)), PARAMETER :: kappa = 0.4
1108 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
1109 ! real(KIND(1D0)) :: phim_hat
1110 ! real(KIND(1D0)) :: zd_RSL
1111
1112 REAL(KIND(1D0)) :: betaN2
1113 ! INTEGER :: it
1114 ! real(KIND(1D0)) :: err
1115 ! real(KIND(1D0)) :: phim
1116
1117 ! betaN for trees found to be 0.3 and for urban 0.4 linearly interpolate between the two using surface fractions
1118 ! betaN2 = 0.30 + (1.-sfr_surf(ConifSurf) - sfr_surf(ConifSurf))*0.1
1119 IF (pai > 0) THEN
1120 betan2 = 0.30*sfr_tr/pai + (pai - sfr_tr)/pai*0.4
1121 ELSE
1122 betan2 = 0.35
1123 END IF
1124
1125 betahf = cal_beta_lc(stabilitymethod, betan2, lc_over_l)
1126
1127 betanl = cal_beta_lc(stabilitymethod, kappa/2., lc_over_l)
1128
1129 IF (lc_over_l > a2) THEN
1130 beta = betahf
1131 ELSE
1132 beta = betanl + ((betahf - betanl)/(1.+a1*abs(lc_over_l - a2)**a3))
1133 END IF
1134
1135 IF (beta > 0.5) THEN
1136 beta = 0.5
1137 END IF
1138

References allocatearray::a1, allocatearray::a2, allocatearray::a3, and cal_beta_lc().

Referenced by rsl_cal_prms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_ch()

subroutine rsl_module::cal_ch ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  zd_RSL,
real(kind(1d0)), intent(in)  Lc,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), intent(in)  Scc,
real(kind(1d0)), intent(in)  f,
real(kind(1d0)), intent(out)  c2h,
real(kind(1d0)), intent(out)  ch 
)

Definition at line 655 of file suews_phys_rslprof.f95.

656
657 IMPLICIT NONE
658 INTEGER, INTENT(in) :: StabilityMethod ! stability method
659 ! real(KIND(1D0)), intent(in) :: z ! height of interest [m]
660 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
661 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
662 REAL(KIND(1D0)), INTENT(in) :: Scc !
663 REAL(KIND(1D0)), INTENT(in) :: f !
664 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
665 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
666 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
667
668 ! output
669 REAL(KIND(1D0)), INTENT(out) :: ch
670 REAL(KIND(1D0)), INTENT(out) :: c2h ! displacement height used in RSL
671
672 ! internal variables
673 REAL(KIND(1D0)) :: phih_zh ! displacement height used in RSL
674 REAL(KIND(1D0)) :: phih_zhdz ! displacement height used in RSL
675 REAL(KIND(1D0)) :: dphih ! displacement height used in RSL
676 REAL(KIND(1D0)) :: phi_hathZh ! displacement height used in RSL
677
678 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
679 REAL(KIND(1D0)), PARAMETER :: dz = 0.1 !height step
680
681 phih_zh = stab_phi_heat(stabilitymethod, (zh_rsl - zd_rsl)/l_mod)
682 phih_zhdz = stab_phi_heat(stabilitymethod, (zh_rsl - zd_rsl + 1.)/l_mod)
683
684 dphih = phih_zhdz - phih_zh
685 IF (phih_zh /= 0.) THEN
686 phi_hathzh = kappa*scc/(2.*beta*phih_zh)
687 ELSE
688 phi_hathzh = 1.
689 END IF
690
691 IF (phi_hathzh >= 1.) THEN
692 ! more stable, but less correct
693 c2h = 0.5
694 phi_hathzh = 1.
695 ELSE
696 ! if very unstable this might cause some high values of psihat_z
697 c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/phih_zh)))/(2.*beta*phih_zh - kappa*scc)
698 END IF
699 ! force c2h to 0.5 for better stability. TS 14 Jul 2020
700 ! TODO: a more proper threshold needs to be determined
701 c2h = 0.5
702
703 ch = (1.-phi_hathzh)*exp(c2h/2.)
704

References atmmoiststab_module::stab_phi_heat().

Referenced by rsl_cal_prms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_cm()

subroutine rsl_module::cal_cm ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  zd_RSL,
real(kind(1d0)), intent(in)  Lc,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), intent(out)  c2,
real(kind(1d0)), intent(out)  cm 
)

Definition at line 602 of file suews_phys_rslprof.f95.

603
604 IMPLICIT NONE
605 INTEGER, INTENT(in) :: StabilityMethod ! stability method
606 ! real(KIND(1D0)), intent(in) :: z ! height of interest [m]
607 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
608 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
609 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
610 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
611 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
612
613 ! output
614 REAL(KIND(1D0)), INTENT(out) :: c2
615 REAL(KIND(1D0)), INTENT(out) :: cm
616
617 ! internal variables
618 ! real(KIND(1D0)) ::phim_zh
619 REAL(KIND(1D0)) :: phi_hatmZh
620 REAL(KIND(1D0)) :: phim_zh
621 REAL(KIND(1D0)) :: phim_zhdz
622 REAL(KIND(1D0)) :: dphi
623 ! real(KIND(1D0)) ::phi_hatmZh
624
625 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
626 REAL(KIND(1D0)), PARAMETER :: dz = 0.1 !height step
627
628 phim_zh = stab_phi_mom(stabilitymethod, (zh_rsl - zd_rsl)/l_mod)
629 phim_zhdz = stab_phi_mom(stabilitymethod, (zh_rsl - zd_rsl + dz)/l_mod)
630
631 dphi = (phim_zhdz - phim_zh)/dz
632 IF (phim_zh /= 0.) THEN
633 phi_hatmzh = kappa/(2.*beta*phim_zh)
634 ELSE
635 ! neutral condition
636 phi_hatmzh = 1.
637 END IF
638
639 IF (phi_hatmzh >= 1.) THEN
640 ! more stable, but less correct
641 c2 = 0.5
642 phi_hatmzh = 1.
643 ELSE
644 ! if very unstable this might cause some high values of psihat_z
645 c2 = (kappa*(3.-(2.*beta**2.*lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
646 END IF
647 ! force c2 to 0.5 for better stability. TS 14 Jul 2020
648 ! TODO: a more proper threshold needs to be determined
649 c2 = 0.5
650
651 cm = (1.-phi_hatmzh)*exp(c2/2.)
652

References atmmoiststab_module::stab_phi_mom().

Referenced by cal_phim_hat(), and rsl_cal_prms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_elm_rsl()

real(kind(1d0)) function rsl_module::cal_elm_rsl ( real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 414 of file suews_phys_rslprof.f95.

415
416 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
417 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
418
419 ! output
420 REAL(KIND(1D0)) :: elm ! a scaling parameter for RSL
421
422 elm = 2.*beta**3*lc
423

Referenced by cal_phim_hat(), and rsl_cal_prms().

Here is the caller graph for this function:

◆ cal_phim_hat()

real(kind(1d0)) function rsl_module::cal_phim_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,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  lc 
)

Definition at line 577 of file suews_phys_rslprof.f95.

578 IMPLICIT NONE
579 INTEGER, INTENT(in) :: StabilityMethod ! stability method
580 REAL(KIND(1D0)), INTENT(in) :: z
581 REAL(KIND(1D0)), INTENT(in) :: zh_RSL
582 REAL(KIND(1D0)), INTENT(in) :: L_MOD
583 REAL(KIND(1D0)), INTENT(in) :: beta
584 REAL(KIND(1D0)), INTENT(in) :: lc
585 REAL(KIND(1D0)) :: phim_hat
586 REAL(KIND(1D0)) :: zd_RSL
587
588 REAL(KIND(1D0)) :: elm
589 REAL(KIND(1D0)) :: c2
590 REAL(KIND(1D0)) :: cm
591
592 elm = cal_elm_rsl(beta, lc)
593
594 zd_rsl = cal_zd_rsl(zh_rsl, beta, lc)
595
596 CALL cal_cm(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod, c2, cm)
597
598 phim_hat = 1 - cm*exp(-1.*c2*beta*(z - zd_rsl)/elm)
599

References cal_cm(), cal_elm_rsl(), and cal_zd_rsl().

Here is the call graph for this function:

◆ cal_psih_hat()

recursive real(kind(1d0)) function rsl_module::cal_psih_hat ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  psihath_top,
real(kind(1d0)), intent(in)  psihath_mid,
real(kind(1d0)), intent(in)  z_top,
real(kind(1d0)), intent(in)  z_mid,
real(kind(1d0)), intent(in)  z_btm,
real(kind(1d0)), intent(in)  ch,
real(kind(1d0)), intent(in)  c2h,
real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  zd_RSL,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  elm,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 501 of file suews_phys_rslprof.f95.

507 ! TS, 23 Oct 2019: calculate psi_hat for momentum
508 ! TS 30 Jun 2022: revised calculation logic for better calculation performance
509 IMPLICIT NONE
510 INTEGER, INTENT(in) :: StabilityMethod ! stability method
511 REAL(KIND(1D0)), INTENT(in) :: psihath_top ! psihath at z_top [m]
512 REAL(KIND(1D0)), INTENT(in) :: z_top ! height at a top level [m]
513 REAL(KIND(1D0)), INTENT(in) :: psihath_mid ! psihath at z_mid [m]
514 REAL(KIND(1D0)), INTENT(in) :: z_mid ! height at a middle level [m]
515 REAL(KIND(1D0)), INTENT(in) :: z_btm ! height of interest [m]
516 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
517 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! displacement height [m]
518 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
519 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
520 REAL(KIND(1D0)), INTENT(in) :: elm ! parameter in RSL
521 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
522 REAL(KIND(1D0)), INTENT(in) :: ch ! parameter in RSL
523 REAL(KIND(1D0)), INTENT(in) :: c2h ! parameter in RSL
524
525 ! output
526 REAL(KIND(1D0)) :: psihath_btm ! psim_hat at height of interest
527
528 ! internal variables
529 REAL(KIND(1D0)) :: phih_btm ! displacement height used in RSL
530 REAL(KIND(1D0)) :: phih_mid ! displacement height used in RSL
531 ! REAL(KIND(1D0)) :: psihatm_btm ! displacement height used in RSL
532 REAL(KIND(1D0)) :: slope_top ! displacement height used in RSL
533 REAL(KIND(1D0)) :: slope_btm ! displacement height used in RSL
534 REAL(KIND(1D0)) :: z_plus ! displacement height used in RSL
535 REAL(KIND(1D0)) :: psihath_plus ! displacement height used in RSL
536
537 REAL(KIND(1D0)), PARAMETER :: tol = 0.1 ! tolerance for iterative calculations
538 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
539 REAL(KIND(1D0)) :: dz_above
540
541 dz_above = z_mid - z_btm
542 ! single step calculation
543 phih_mid = stab_phi_heat(stabilitymethod, (z_mid - zd_rsl)/l_mod)
544 phih_btm = stab_phi_heat(stabilitymethod, (z_btm - zd_rsl)/l_mod)
545
546 psihath_btm = psihath_mid + dz_above/2.*phih_mid*(ch*exp(-1.*c2h*beta*(z_mid - zd_rsl)/elm)) & !Taylor's approximation for integral
547 /(z_mid - zd_rsl)
548 psihath_btm = psihath_btm + dz_above/2.*phih_btm*(ch*exp(-1.*c2h*beta*(z_btm - zd_rsl)/elm)) &
549 /(z_btm - zd_rsl)
550 IF (dz_above/zd_rsl < 1e-2) THEN
551 RETURN ! psihatm_z will be returned
552 END IF
553
554 IF (abs(psihath_btm) > 1e-3) THEN
555 ! test if recursion is needed by comparing slopes of psihat within the top and mid ranges
556 slope_top = (psihath_top - psihath_mid)/(z_top - z_mid)
557 slope_btm = (psihath_mid - psihath_btm)/(z_mid - z_btm)
558 IF (abs(slope_btm) > (1 + tol)*abs(slope_top)) THEN
559 z_plus = (z_mid + z_btm)/2
560 psihath_plus = cal_psih_hat( &
561 stabilitymethod, &
562 psihath_top, psihath_mid, &
563 z_top, z_mid, z_plus, &
564 ch, c2h, &
565 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
566 psihath_btm = cal_psih_hat( &
567 stabilitymethod, &
568 psihath_mid, psihath_plus, &
569 z_mid, z_plus, z_btm, &
570 ch, c2h, &
571 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
572 END IF
573 END IF
574

References cal_psih_hat(), and atmmoiststab_module::stab_phi_heat().

Referenced by cal_psih_hat(), and rsl_cal_prms().

Here is the call graph for this function:
Here is the caller 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)  psihatm_top,
real(kind(1d0)), intent(in)  psihatm_mid,
real(kind(1d0)), intent(in)  z_top,
real(kind(1d0)), intent(in)  z_mid,
real(kind(1d0)), intent(in)  z_btm,
real(kind(1d0)), intent(in)  cm,
real(kind(1d0)), intent(in)  c2,
real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  zd_RSL,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  elm,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 426 of file suews_phys_rslprof.f95.

432 ! TS, 23 Oct 2019: calculate psi_hat for momentum
433 ! TS 30 Jun 2022: revised calculation logic for better calculation performance
434 IMPLICIT NONE
435 INTEGER, INTENT(in) :: StabilityMethod ! stability method
436 REAL(KIND(1D0)), INTENT(in) :: psihatm_top ! height of interest [m]
437 REAL(KIND(1D0)), INTENT(in) :: z_top ! height of interest [m]
438 REAL(KIND(1D0)), INTENT(in) :: psihatm_mid ! height of interest [m]
439 REAL(KIND(1D0)), INTENT(in) :: z_mid ! height of interest [m]
440 REAL(KIND(1D0)), INTENT(in) :: z_btm ! height of interest [m]
441 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
442 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
443 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
444 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
445 REAL(KIND(1D0)), INTENT(in) :: elm ! parameter in RSL
446 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
447 REAL(KIND(1D0)), INTENT(in) :: cm ! Obukhov length [m]
448 REAL(KIND(1D0)), INTENT(in) :: c2 ! Obukhov length [m]
449
450 ! output
451 REAL(KIND(1D0)) :: psihatm_btm ! psim_hat at height of interest
452
453 ! internal variables
454 REAL(KIND(1D0)) :: phim_btm ! displacement height used in RSL
455 REAL(KIND(1D0)) :: phim_mid ! displacement height used in RSL
456 REAL(KIND(1D0)) :: slope_top ! displacement height used in RSL
457 REAL(KIND(1D0)) :: slope_btm ! displacement height used in RSL
458 REAL(KIND(1D0)) :: z_plus ! displacement height used in RSL
459 REAL(KIND(1D0)) :: psihatm_plus ! displacement height used in RSL
460
461 REAL(KIND(1D0)), PARAMETER :: tol = 0.5 ! tolerance for iterative calculations
462 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
463 REAL(KIND(1D0)) :: dz_above
464
465 dz_above = z_mid - z_btm
466 ! single step calculation
467 phim_mid = stab_phi_mom(stabilitymethod, (z_mid - zd_rsl)/l_mod)
468 phim_btm = stab_phi_mom(stabilitymethod, (z_btm - zd_rsl)/l_mod)
469
470 psihatm_btm = psihatm_mid + dz_above/2.*phim_mid*(cm*exp(-1.*c2*beta*(z_mid - zd_rsl)/elm)) & !Taylor's approximation for integral
471 /(z_mid - zd_rsl)
472 psihatm_btm = psihatm_btm + dz_above/2.*phim_btm*(cm*exp(-1.*c2*beta*(z_btm - zd_rsl)/elm)) &
473 /(z_btm - zd_rsl)
474 IF (dz_above/zd_rsl < 1e-2) THEN
475 RETURN ! psihatm_z will be returned
476 END IF
477
478 IF (abs(psihatm_btm) > 1e-3) THEN
479 ! test if recursion is needed by comparing slopes of psihat within the top and mid ranges
480 slope_top = (psihatm_top - psihatm_mid)/(z_top - z_mid)
481 slope_btm = (psihatm_mid - psihatm_btm)/(z_mid - z_btm)
482 IF (abs(slope_btm) > (1 + tol)*abs(slope_top)) THEN
483 z_plus = (z_mid + z_btm)/2
484 psihatm_plus = cal_psim_hat( &
485 stabilitymethod, &
486 psihatm_top, psihatm_mid, &
487 z_top, z_mid, z_plus, &
488 cm, c2, &
489 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
490 psihatm_btm = cal_psim_hat( &
491 stabilitymethod, &
492 psihatm_mid, psihatm_plus, &
493 z_mid, z_plus, z_btm, &
494 cm, c2, &
495 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
496 END IF
497 END IF
498

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

Referenced by cal_psim_hat(), and rsl_cal_prms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_z0_rsl()

real(kind(1d0)) function rsl_module::cal_z0_rsl ( integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  zH_RSL,
real(kind(1d0)), intent(in)  zd_RSL,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  L_MOD_RSL,
real(kind(1d0)), intent(in)  psihatm_Zh 
)

Definition at line 877 of file suews_phys_rslprof.f95.

878 ! calculate z0 iteratively
879 ! TS, 23 Oct 2019
880 IMPLICIT NONE
881 INTEGER, INTENT(in) :: StabilityMethod
882 REAL(KIND(1D0)), INTENT(in) :: zH_RSL ! canyon depth [m]
883 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! displacement height [m]
884 REAL(KIND(1D0)), INTENT(in) :: L_MOD_RSL ! Monin Obukhov length[m]
885 ! REAL(KIND(1D0)), INTENT(in) :: Lc ! canyon length scale [m]
886 REAL(KIND(1D0)), INTENT(in) :: beta ! height scale for bluff bodies [m]
887 REAL(KIND(1D0)), INTENT(in) :: psihatm_Zh ! psihatm at zH
888
889 ! output
890 REAL(KIND(1D0)) :: z0_RSL
891
892 ! internal variables
893 REAL(KIND(1D0)) :: psimZh, psimz0, z0_RSL_x
894 REAL(KIND(1D0)) :: err
895 INTEGER :: it
896
897 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
898 ! REAL(KIND(1d0)), PARAMETER::r = 0.1
899 ! REAL(KIND(1d0)), PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
900
901 psimzh = stab_psi_mom(stabilitymethod, (zh_rsl - zd_rsl)/l_mod_rsl)
902 ! psihatm_Zh = cal_psim_hat(StabilityMethod, Zh_RSL, zh_RSL, L_MOD_RSL, beta, Lc)
903
904 !first guess
905 z0_rsl = 0.1*zh_rsl
906 err = 10.
907 ! psimz0 = 0.5
908 it = 1
909 DO WHILE ((err > 0.001) .AND. (it < 10))
910 z0_rsl_x = z0_rsl
911 psimz0 = stab_psi_mom(stabilitymethod, z0_rsl_x/l_mod_rsl)
912 z0_rsl = (zh_rsl - zd_rsl)*exp(-1.*kappa/beta)*exp(-1.*psimzh + psimz0)*exp(psihatm_zh)
913 err = abs(z0_rsl_x - z0_rsl)
914 it = it + 1
915 END DO
916
917 ! set limit on z0_RSL for numeric stability
918 z0_rsl = merge(z0_rsl, 1d-2, z0_rsl > 1d-2)
919

References atmmoiststab_module::stab_psi_mom().

Referenced by rsl_cal_prms().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_zd_rsl()

real(kind(1d0)) function rsl_module::cal_zd_rsl ( real(kind(1d0)), intent(in)  zh_RSL,
real(kind(1d0)), intent(in)  beta,
real(kind(1d0)), intent(in)  Lc 
)

Definition at line 862 of file suews_phys_rslprof.f95.

863
864 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
865 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
866 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
867
868 ! output
869 REAL(KIND(1D0)) :: zd_RSL ! zd used in RSL
870
871 zd_rsl = zh_rsl - (beta**2.)*lc
872 !correct negative values using rule of thumb, TS 24 Jun 2020
873 ! if (zd_RSL < 0) zd_RSL = 0.7*Zh_RSL
874

Referenced by cal_phim_hat(), and rsl_cal_prms().

Here is the caller graph for this function:

◆ interp_z()

real(kind(1d0)) function rsl_module::interp_z ( real(kind(1d0)), intent(in)  z_x,
real(kind(1d0)), dimension(nz), intent(in)  z,
real(kind(1d0)), dimension(nz), intent(in)  v 
)

Definition at line 375 of file suews_phys_rslprof.f95.

376
377 REAL(KIND(1D0)), INTENT(in) :: z_x ! height to interpolate at
378 REAL(KIND(1D0)), DIMENSION(nz), INTENT(in) :: z ! heights
379 REAL(KIND(1D0)), DIMENSION(nz), INTENT(in) :: v ! values associated with heights
380
381 ! output
382 REAL(KIND(1D0)) :: v_x ! zd used in RSL
383
384 ! local variables
385 REAL(KIND(1D0)) :: slope ! slope
386 REAL(KIND(1D0)) :: dz ! slope
387 REAL(KIND(1D0)), DIMENSION(nz) :: dif ! slope
388 INTEGER :: idx_low ! vertical index lower than z_x
389 INTEGER :: idx_x ! vertical index lower than z_x
390 INTEGER :: idx_high ! vertical index higher than z_x
391 ! INTEGER :: idx ! vertical index higher than z_x
392 INTEGER, PARAMETER :: nz = 30 ! vertical index higher than z_x
393
394 ! initialise variables
395 idx_x = 0
396
397 dif = z - z_x
398 idx_x = maxloc(dif, 1, abs(dif) < 1.d-6)
399 idx_low = maxloc(dif, 1, dif < 0.)
400 idx_high = minloc(dif, 1, dif > 0.)
401
402 IF (idx_x > 0) THEN
403 ! z_x is one of zarray elements
404 v_x = v(idx_x)
405 ELSE
406 ! linear interpolation is performed
407 dz = z(idx_high) - z(idx_low)
408 slope = (v(idx_high) - v(idx_low))/dz
409 v_x = v(idx_low) + (z_x - z(idx_low))*slope
410 END IF
411

Referenced by rslprofile().

Here is the caller graph for this function:

◆ rsl_cal_prms()

subroutine rsl_module::rsl_cal_prms ( integer, intent(in)  StabilityMethod,
integer, intent(in)  nz_above,
real(kind(1d0)), dimension(nz_above), intent(in)  z_array,
real(kind(1d0)), intent(in)  zh,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(in)  FAI,
real(kind(1d0)), intent(in)  PAI,
real(kind(1d0)), dimension(nz_above), intent(out)  psihatm_array,
real(kind(1d0)), dimension(nz_above), intent(out)  psihath_array,
real(kind(1d0)), intent(out)  zH_RSL,
real(kind(1d0)), intent(out)  L_MOD_RSL,
real(kind(1d0)), intent(out)  Lc,
real(kind(1d0)), intent(out)  beta,
real(kind(1d0)), intent(out)  zd_RSL,
real(kind(1d0)), intent(out)  z0_RSL,
real(kind(1d0)), intent(out)  elm,
real(kind(1d0)), intent(out)  Scc,
real(kind(1d0)), intent(out)  fx 
)

Definition at line 922 of file suews_phys_rslprof.f95.

925
926 IMPLICIT NONE
927 INTEGER, INTENT(in) :: StabilityMethod ! stability method
928 INTEGER, INTENT(IN) :: nz_above ! number of layers above zh
929 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(in) :: z_array ! land cover fractions
930 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(out) :: psihatm_array ! land cover fractions
931 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(out) :: psihath_array ! land cover fractions
932 REAL(KIND(1D0)), INTENT(in) :: zh ! canyon depth [m]
933 REAL(KIND(1D0)), INTENT(in) :: FAI ! frontal area index inlcuding trees
934 ! REAL(KIND(1D0)), INTENT(in) :: FAIBldg ! frontal area index of buildings
935 REAL(KIND(1D0)), INTENT(in) :: PAI ! plan area index inlcuding area of trees
936 ! REAL(KIND(1D0)), INTENT(in) :: porosity_dectr ! porosity of deciduous trees
937 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
938 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! land cover fractions
939
940 ! output
941 ! real(KIND(1D0)), intent(out) ::L_stab ! threshold for Obukhov length under stable conditions
942 ! real(KIND(1D0)), intent(out) ::L_unstab ! threshold for Obukhov length under unstable conditions
943 REAL(KIND(1D0)), INTENT(out) :: L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
944 REAL(KIND(1D0)), INTENT(out) :: zH_RSL ! mean canyon height used in RSL module with thresholds applied
945 ! real(KIND(1D0)), intent(out) ::Lc_stab ! threshold for penetration distance scale under stable conditions
946 ! real(KIND(1D0)), intent(out) ::Lc_unstab ! threshold for penetration distance scale under unstable conditions
947 REAL(KIND(1D0)), INTENT(out) :: Lc ! penetration distance scale for bluff bodies [m]
948 REAL(KIND(1D0)), INTENT(out) :: beta ! psim_hat at height of interest
949 REAL(KIND(1D0)), INTENT(out) :: zd_RSL ! displacement height to prescribe if necessary [m]
950 REAL(KIND(1D0)), INTENT(out) :: z0_RSL ! roughness length [m]
951 REAL(KIND(1D0)), INTENT(out) :: elm ! length scale used in RSL
952 REAL(KIND(1D0)), INTENT(out) :: Scc ! parameter in RSL
953 REAL(KIND(1D0)), INTENT(out) :: fx ! parameter in RSL
954
955 ! internal variables
956 INTEGER :: iz
957 REAL(KIND(1D0)) :: sfr_tr
958 REAL(KIND(1D0)) :: psihatm_Zh
959 ! real(KIND(1D0)) ::L_MOD_RSL_x
960 ! real(KIND(1D0)) ::lc_x
961 REAL(KIND(1D0)) :: lc_over_L
962 REAL(KIND(1D0)) :: z_top, z_mid, z_btm
963 REAL(KIND(1D0)) :: psihatm_top, psihatm_mid, psihatm_btm
964 REAL(KIND(1D0)) :: psihath_top, psihath_mid, psihath_btm
965 REAL(KIND(1D0)) :: cm, ch, c2m, c2h
966 ! real(KIND(1D0)) ::betaHF
967 ! real(KIND(1D0)) ::betaNL
968 REAL(KIND(1D0)) :: Lc_min ! LB Oct2021 - minimum value of Lc
969 REAL(KIND(1D0)) :: dim_bluffbody ! LB Oct2021 - horizontal building dimensions
970
971 REAL(KIND(1D0)), PARAMETER :: cd_tree = 1.2 ! drag coefficient tree canopy !!!!needs adjusting!!!
972 REAL(KIND(1D0)), PARAMETER :: a_tree = 0.05 ! the foliage area per unit volume !!!!needs adjusting!!!
973 ! lv_J_kg = 2.5E6, &! latent heat for water vapor!!! make consistant with rest of code
974 REAL(KIND(1D0)), PARAMETER :: beta_N = 0.40 ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
975 REAL(KIND(1D0)), PARAMETER :: pi = 4.*atan(1.0)
976
977 REAL(KIND(1D0)), PARAMETER :: planF_low = 1e-6
978 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40 ! von karman constant
979 ! REAL(KIND(1d0)), PARAMETER::z0m= 0.40
980 REAL(KIND(1D0)), PARAMETER :: r = 0.1
981 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
982 REAL(KIND(1D0)), PARAMETER :: Zh_min = 0.4 ! limit for minimum canyon height used in RSL module
983 REAL(KIND(1D0)), PARAMETER :: porosity_evetr = 0.32 ! assumed porosity of evergreen trees, ref: Lai et al. (2022), http://dx.doi.org/10.2139/ssrn.4058842
984
985 ! under stable conditions, set a threshold for L_MOD to avoid numerical issues. TS 28 Oct 2019
986 ! L_MOD = merge(L_MOD, 300.d1, L_MOD < 300.)
987
988 ! zH_RSL
989 zh_rsl = max(zh, zh_min)
990
991 ! ! land cover fraction of bluff bodies
992 ! PAI = DOT_PRODUCT(sfr_surf([BldgSurf, ConifSurf, DecidSurf]), [1D0, 1 - porosity_evetr, 1 - porosity_dectr])
993 ! set a threshold for sfr_zh to avoid numerical difficulties
994 ! PAI = min(PAI, 0.8)
995
996 ! land cover fraction of trees
997 sfr_tr = sum(sfr_surf([conifsurf, decidsurf]))
998
999 ! height scale for buildings !not used? why?
1000 ! Lc_build = (1.-sfr_surf(BldgSurf))/FAI*Zh_RSL ! Coceal and Belcher 2004 assuming Cd = 2
1001
1002 ! height scale for tress
1003 ! Lc_tree = 1./(cd_tree*a_tree) ! not used? why?
1004
1005 ! height scale for bluff bodies
1006 ! LB Oct2021 - replaced PAI with sfr(BldgSurf) since the parameter should represent the solid fraction (and trees have negligible solid fraction).
1007 ! TS 18 Jun 2022 - changed sfr_surf(BldgSurf) back to PAI to account for trees with PAI updated by accounting for porosity.
1008 lc = (1.-pai)/fai*zh_rsl
1009
1010 ! set a minimum threshold (of 0.5*Zh_RSL) for Lc to avoid numerical diffulties when FAI is too large (e.g., FAI>10)
1011 ! Lc = MERGE(Lc, 0.5*Zh_RSL, Lc > 0.5*Zh_RSL)
1012 ! LB Oct2021 - set a minimum Lc threshold based on the Lc required to ensure the horizontal length scale associated with changes in canopy geometry (i.e. 3Lc) is greater than a typical street+building unit
1013 ! Note: the horizontal building size and street+building unit size is calculated assuming a regular array of cuboids with the same x and y dimension but with height that can be different
1014 dim_bluffbody = zh_rsl*pai/fai
1015 lc_min = dim_bluffbody*pai**(-0.5)/3.
1016 lc = merge(lc, lc_min, lc > lc_min)
1017
1018 ! a normalised scale with a physcially valid range between [-2,2] (Harman 2012, BLM)
1019 lc_over_l = lc/l_mod
1020 ! lc_over_L = lc/L_MOD_RSL_x
1021 IF (lc_over_l > 0) THEN
1022 lc_over_l = min(2., lc_over_l)
1023 ELSE
1024 lc_over_l = max(-2., lc_over_l)
1025 END IF
1026 ! correct L_MOD_RSL
1027 l_mod_rsl = lc/lc_over_l
1028
1029 ! Step 2:
1030 ! Parameterise beta according to Harman 2012 with upper limit of 0.5
1031 beta = cal_beta_rsl(stabilitymethod, pai, sfr_tr, lc_over_l)
1032
1033 ! Schmidt number Harman and Finnigan 2008: assuming the same for heat and momemntum
1034 scc = 0.5 + 0.3*tanh(2.*lc_over_l)
1035 fx = 0.5*((1.+4.*r*scc)**0.5) - 0.5
1036
1037 ! zero displacement height used in RSL
1038 zd_rsl = cal_zd_rsl(zh_rsl, beta, lc)
1039
1040 ! scaling parameter for RSL
1041 elm = cal_elm_rsl(beta, lc)
1042
1043 CALL cal_ch(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod_rsl, scc, fx, c2h, ch)
1044 CALL cal_cm(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod_rsl, c2m, cm)
1045
1046 ! calculate psihat values at desirable heights
1047 psihatm_top = 0
1048 psihatm_mid = 0
1049 psihatm_array(nz_above) = 0
1050 psihatm_array(nz_above - 1) = 0
1051
1052 psihath_top = 0
1053 psihath_mid = 0
1054 psihath_array(nz_above) = 0
1055 psihath_array(nz_above - 1) = 0
1056
1057 DO iz = nz_above, 3, -1
1058 z_top = z_array(iz)
1059 z_mid = z_array(iz - 1)
1060 z_btm = z_array(iz - 2)
1061
1062 ! momentum
1063 psihatm_btm = cal_psim_hat(stabilitymethod, &
1064 psihatm_top, psihatm_mid, &
1065 z_top, z_mid, z_btm, &
1066 cm, c2m, &
1067 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
1068 psihatm_array(iz - 2) = psihatm_btm
1069 psihatm_top = psihatm_mid
1070 psihatm_mid = psihatm_btm
1071
1072 ! heat
1073 psihath_btm = cal_psih_hat(stabilitymethod, &
1074 psihath_top, psihath_mid, &
1075 z_top, z_mid, z_btm, &
1076 ch, c2h, &
1077 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
1078 psihath_top = psihath_mid
1079 psihath_mid = psihath_btm
1080 psihath_array(iz - 2) = psihath_btm
1081
1082 END DO
1083
1084 ! calculate z0 iteratively
1085 psihatm_zh = psihatm_array(1)
1086 z0_rsl = cal_z0_rsl(stabilitymethod, zh_rsl, zd_rsl, beta, l_mod_rsl, psihatm_zh)
1087

References allocatearray::a3, cal_beta_rsl(), cal_ch(), cal_cm(), cal_elm_rsl(), cal_psih_hat(), cal_psim_hat(), cal_z0_rsl(), cal_zd_rsl(), allocatearray::conifsurf, and allocatearray::decidsurf.

Referenced by rslprofile().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ rslprofile()

subroutine rsl_module::rslprofile ( integer, intent(in)  DiagMethod,
real(kind(1d0)), intent(in)  Zh,
real(kind(1d0)), intent(in)  z0m,
real(kind(1d0)), intent(in)  zdm,
real(kind(1d0)), intent(in)  z0v,
real(kind(1d0)), intent(in)  L_MOD,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), intent(in)  FAI,
real(kind(1d0)), intent(in)  PAI,
integer, intent(in)  StabilityMethod,
real(kind(1d0)), intent(in)  RA_h,
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(ncolumnsdataoutrsl - 5), intent(out)  dataoutLineRSL 
)

Definition at line 13 of file suews_phys_rslprof.f95.

22 !-----------------------------------------------------
23 ! calculates windprofiles using MOST with a RSL-correction
24 ! based on Harman & Finnigan 2007
25 !
26 ! last modified by:
27 ! NT 16 Mar 2019: initial version
28 ! TS 16 Oct 2019: improved consistency in parameters/varaibles within SUEWS
29 ! TODO how to improve the speed of this code
30 !
31 !-----------------------------------------------------
32
33 IMPLICIT NONE
34
35 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! surface fractions [-]
36 REAL(KIND(1D0)), INTENT(in) :: zMeas ! height of atmospheric forcing [m]
37 REAL(KIND(1D0)), INTENT(in) :: avU1 ! Wind speed at forcing height [m s-1]
38 REAL(KIND(1D0)), INTENT(in) :: Temp_C ! Air temperature at forcing height [C]
39 REAL(KIND(1D0)), INTENT(in) :: avRH ! relative humidity at forcing height [-]
40 REAL(KIND(1D0)), INTENT(in) :: Press_hPa ! pressure at forcing height [hPa]
41 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
42 REAL(KIND(1D0)), INTENT(in) :: RA_h ! aerodynamic resistance for heat [s m-1]
43 REAL(KIND(1D0)), INTENT(in) :: avcp ! specific heat capacity [J kg-1 K-1]
44 REAL(KIND(1D0)), INTENT(in) :: lv_J_kg ! Latent heat of vaporization in [J kg-1]
45 REAL(KIND(1D0)), INTENT(in) :: avdens ! air density [kg m-3]
46 REAL(KIND(1D0)), INTENT(in) :: qh ! sensible heat flux [W m-2]
47 REAL(KIND(1D0)), INTENT(in) :: qe ! Latent heat flux [W m-2]
48 REAL(KIND(1D0)), INTENT(in) :: Zh ! Mean building height [m]
49 REAL(KIND(1D0)), INTENT(in) :: z0m ! roughness for momentum [m]
50 REAL(KIND(1D0)), INTENT(in) :: z0v ! roughnesslength for heat [s m-1]
51 REAL(KIND(1D0)), INTENT(in) :: zdm ! zero-plane displacement [m]
52 REAL(KIND(1D0)), INTENT(in) :: FAI ! Frontal area index [-]
53 REAL(KIND(1D0)), INTENT(in) :: PAI ! Plan area index [-]
54 ! REAL(KIND(1D0)), INTENT(in) :: FAIBldg ! Frontal area index of buildings [-]
55 ! REAL(KIND(1D0)), INTENT(in) :: porosity_dectr ! porosity of deciduous trees [-]
56
57 INTEGER, INTENT(in) :: StabilityMethod
58 INTEGER, INTENT(in) :: DiagMethod
59
60 REAL(KIND(1D0)), INTENT(out) :: T2_C ! Air temperature at 2 m [C]
61 REAL(KIND(1D0)), INTENT(out) :: q2_gkg ! Air specific humidity at 2 m [g kg-1]
62 REAL(KIND(1D0)), INTENT(out) :: U10_ms ! wind speed at 10 m [m s-1]
63 REAL(KIND(1D0)), INTENT(out) :: RH2 ! Air relative humidity [-]
64
65 INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
66
67 REAL(KIND(1D0)), PARAMETER :: cd_tree = 1.2 ! drag coefficient tree canopy !!!!needs adjusting!!!
68 REAL(KIND(1D0)), PARAMETER :: a_tree = 0.05 ! the foliage area per unit volume !!!!needs adjusting!!!
69 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40 ! von karman constant
70
71 REAL(KIND(1D0)), PARAMETER :: beta_N = 0.40 ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
72 REAL(KIND(1D0)), PARAMETER :: pi = 4.*atan(1.0), r = 0.1
73 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
74
75 REAL(KIND(1D0)), PARAMETER :: porosity_evetr = 0.32 ! assumed porosity of evergreen trees, ref: Lai et al. (2022), http://dx.doi.org/10.2139/ssrn.4058842
76
77 ! Variables array [z,U,T,q, 12 debug vars]
78 ! z: height array
79 ! U,T,q: wind speed, air temp, specific humidity at z;
80 ! debug vars: see dataoutLineRSL
81 REAL(KIND(1D0)), INTENT(out), DIMENSION(ncolumnsDataOutRSL - 5) :: dataoutLineRSL
82 REAL(KIND(1D0)), DIMENSION(nz) :: psihatm_z
83 REAL(KIND(1D0)), DIMENSION(nz) :: psihath_z
84 ! REAL(KIND(1D0)), DIMENSION(nz) :: dif
85 ! REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z, psihath_z
86 REAL(KIND(1D0)), DIMENSION(nz) :: zarray
87 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineURSL ! wind speed array [m s-1]
88 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineTRSL ! Temperature array [C]
89 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineqRSL ! Specific humidity array [g kg-1]
90
91 REAL(KIND(1D0)) :: z0_RSL ! roughness length from H&F
92 REAL(KIND(1D0)) :: zd_RSL ! zero-plane displacement
93
94 ! REAL(KIND(1d0))::Lc_build, Lc_tree, Lc ! canopy drag length scale
95 REAL(KIND(1D0)) :: Lc ! canopy drag length scale
96 ! REAL(KIND(1d0))::Lc_stab ! threshold of canopy drag length scale under stable conditions
97 ! REAL(KIND(1d0))::Lc_unstab ! threshold of canopy drag length scale under unstable conditions
98 REAL(KIND(1D0)) :: Scc ! Schmidt number for temperature and humidity
99 REAL(KIND(1D0)) :: psimz, psimz0, psimza, psihz, psihz0, psihza ! stability function for momentum
100 ! REAL(KIND(1d0))::betaHF, betaNL, beta, betaN2 ! beta coefficient from Harman 2012
101 REAL(KIND(1D0)) :: beta ! beta coefficient from Harman 2012
102 REAL(KIND(1D0)) :: elm ! mixing length
103 ! REAL(KIND(1d0))::xxm1, xxm1_2, xxh1, xxh1_2, dphi, dphih ! dummy variables for stability functions
104 REAL(KIND(1D0)) :: fx ! H&F'07 and H&F'08 'constants'
105 REAL(KIND(1D0)) :: t_h, q_h ! H&F'08 canopy corrections
106 REAL(KIND(1D0)) :: TStar_RSL ! temperature scale
107 REAL(KIND(1D0)) :: UStar_RSL ! friction velocity used in RSL
108 REAL(KIND(1D0)) :: UStar_heat ! friction velocity derived from RA_h with correction/restriction
109 ! REAL(KIND(1D0)) :: PAI ! plan area index, including areas of roughness elements: buildings and trees
110 ! REAL(KIND(1d0))::sfr_tr ! land cover fraction of trees
111 REAL(KIND(1D0)) :: L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
112 ! real(KIND(1D0))::L_stab ! threshold for Obukhov length under stable conditions
113 ! real(KIND(1D0))::L_unstab ! threshold for Obukhov length under unstable conditions
114
115 REAL(KIND(1D0)) :: zH_RSL ! mean canyon height used in RSL module with thresholds applied
116 REAL(KIND(1D0)) :: dz_above ! height step above canopy
117 REAL(KIND(1D0)) :: dz_can ! height step within canopy
118 ! REAL(KIND(1D0)) :: phi_hatmZh, phim_zh
119 ! REAL(KIND(1d0)), parameter::zH_min = 8! limit for minimum canyon height used in RSL module
120 ! REAL(KIND(1D0)), PARAMETER :: ratio_dz = 1.618 ! ratio between neighbouring height steps
121
122 REAL(KIND(1D0)) :: qa_gkg, qStar_RSL ! specific humidity scale
123 INTEGER :: I, z
124 INTEGER :: nz_can ! number of heights in canyon
125 INTEGER :: nz_above ! number of heights above canyon
126
127 LOGICAL :: flag_RSL ! whether RSL correction is used
128
129 ! CHARACTER(len=1024) :: Errmessage
130 ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
131 ! Step 1: determine if RSL should be used
132 ! Step 2: determine vertical levels used in RSL
133 ! Step 3: calculate the stability dependent H&F constants
134 ! Step 4: determine psihat at levels above the canopy
135 ! Step 5: Calculate z0 iteratively
136 ! Step 6: Calculate mean variables above canopy
137 ! Step 7: Calculate mean variables in canopy
138
139! ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
140! CALL RSL_cal_prms( &
141! StabilityMethod, & !input
142! zh, L_MOD, sfr_surf, FAI, PAI, & !input
143! zH_RSL, L_MOD_RSL, & ! output
144! Lc, beta, zd_RSL, z0_RSL, elm, Scc, fx)
145
146 ! Step 1: determine if RSL should be used
147
148 IF (diagmethod == 0) THEN
149 ! force MOST to be used
150 flag_rsl = .false.
151 ELSEIF (diagmethod == 1) THEN
152 ! force RSL to be used
153 flag_rsl = .true.
154 ELSEIF (diagmethod == 2) THEN
155
156 flag_rsl = &
157 ! zd>0 subject to FAI > beta**2*(1-PAI); also beta<0.5;
158 ! hence the lower limit of FAI below
159 fai > 0.25*(1 - pai) .AND. fai < 0.45 .AND. & ! FAI
160 ! PAI > 0.1 .AND. PAI < 0.61 .AND. & ! PAI
161 pai > 0.1 .AND. pai < 0.68 .AND. & ! PAI
162 zh > 2 ! effective canopy height
163 ! LB Oct2021 - FAI and PAI can be larger than 0.45 and 0.61 respectively -> remove (1.-PAI)/FAI > .021 constraint
164 !(note: it seems wrong anyway - should be 0.87 not 0.021 based on G&O1991 numbers)
165 ELSE
166 ! default is to use MOST
167 flag_rsl = .false.
168 END IF
169
170 !
171 ! ! Step 2
172 ! ! determine vertical levels used in RSL
173 ! Define the height array with consideration of key heights
174 ! set number of heights within canopy
175 IF (zh <= 2) THEN
176 nz_can = 3
177 ELSE IF (zh <= 10) THEN
178 nz_can = 6
179 ELSE
180 nz_can = 15
181 END IF
182 ! fill up heights in canopy
183 dz_can = zh/nz_can
184 DO i = 1, nz_can
185 zarray(i) = dz_can*i
186 END DO
187
188 ! guaranttee 2 m is within the zarray
189 IF (dz_can > 2) zarray(1) = 1.999
190
191 ! fill up heights above canopy
192 nz_above = nz - nz_can
193 dz_above = (zmeas - zh)/nz_above
194 DO i = nz_can + 1, nz
195 zarray(i) = zh + (i - nz_can)*dz_above
196 END DO
197
198 ! ! determine index at the canyon top
199 ! DO z = 1, nz
200 ! dif(z) = ABS(zarray(z) - Zh)
201 ! END DO
202 ! nz_can = MINLOC(dif, DIM=1)
203 ! ! zarray(idx_can+2) = Zh+.1
204 ! ! zarray(idx_can+1) = Zh+.05
205 ! zarray(nz_can) = Zh
206 ! zarray(idx_can-1) = Zh-.1
207
208 ! determine index at measurement height
209 ! nz = nz
210 zarray(nz) = zmeas
211
212 IF (flag_rsl) THEN
213
214 ! use RSL approach to calculate correction factors
215 psihatm_z = 0.*zarray
216 psihath_z = 0.*zarray
217 ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
218 CALL rsl_cal_prms( &
219 stabilitymethod, & !input
220 nz_above, zarray(nz_can + 1:nz), & !input
221 zh, l_mod, sfr_surf, fai, pai, & !input
222 psihatm_z(nz_can + 1:nz), psihath_z(nz_can + 1:nz), & !output
223 zh_rsl, l_mod_rsl, & ! output
224 lc, beta, zd_rsl, z0_rsl, elm, scc, fx)
225
226 ! Step 3: calculate the stability dependent H&F constants
227
228 ! CALL cal_ch(StabilityMethod, zh_RSL, zd_RSL, Lc, beta, L_MOD_RSL, Scc, fx, c2h, ch)
229 ! CALL cal_cm(StabilityMethod, zH_RSL, zd_RSL, Lc, beta, L_MOD_RSL, c2m, cm)
230
231 ! ! Step 4: determine psihat at levels above the canopy
232 ! DO z = nz - 1, idx_can, -1
233 ! phimz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
234 ! phimzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
235 ! phihz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
236 ! phihzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
237
238 ! psihatm_z(z) = psihatm_z(z + 1) + dz_above/2.*phimzp*(cm*EXP(-1.*c2m*beta*(zarray(z + 1) - zd_RSL)/elm)) & !Taylor's approximation for integral
239 ! /(zarray(z + 1) - zd_RSL)
240 ! psihatm_z(z) = psihatm_z(z) + dz_above/2.*phimz*(cm*EXP(-1.*c2m*beta*(zarray(z) - zd_RSL)/elm)) &
241 ! /(zarray(z) - zd_RSL)
242 ! psihath_z(z) = psihath_z(z + 1) + dz_above/2.*phihzp*(ch*EXP(-1.*c2h*beta*(zarray(z + 1) - zd_RSL)/elm)) & !Taylor's approximation for integral
243 ! /(zarray(z + 1) - zd_RSL)
244 ! psihath_z(z) = psihath_z(z) + dz_above/2.*phihz*(ch*EXP(-1.*c2h*beta*(zarray(z) - zd_RSL)/elm)) &
245 ! /(zarray(z) - zd_RSL)
246 ! END DO
247
248 ELSE
249
250 ! correct parameters if RSL approach doesn't apply for scenario of isolated flows
251 ! see Fig 1 of Grimmond and Oke (1999)
252 ! when isolated flow or skimming flow, implying RSL doesn't apply, force RSL correction to zero
253 psihatm_z = 0
254 psihath_z = 0
255
256 ! use L_MOD as in other parts of SUEWS
257 l_mod_rsl = l_mod
258
259 !correct RSL-based using SUEWS system-wide values
260 z0_rsl = z0m
261 zd_rsl = zdm
262 zh_rsl = zh
263
264 lc = -999
265 beta = -999
266 scc = -999
267 fx = -999
268 elm = -999
269
270 ! then MOST recovers from RSL correction
271 END IF
272
273 ! Step 6: Calculate mean variables above canopy
274 !
275 psimz0 = stab_psi_mom(stabilitymethod, z0_rsl/l_mod_rsl)
276 psimza = stab_psi_mom(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
277 psihza = stab_psi_heat(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
278
279 ustar_rsl = avu1*kappa/(log((zmeas - zd_rsl)/z0_rsl) - psimza + psimz0 + psihatm_z(nz))
280
281 ! TS 11 Feb 2021: limit UStar and TStar to reasonable ranges
282 ! under all conditions, min(UStar)==0.001 m s-1 (Jimenez et al 2012, MWR, https://doi.org/10.1175/mwr-d-11-00056.1
283 ustar_rsl = max(0.001, ustar_rsl)
284 ! under convective/unstable condition, min(UStar)==0.15 m s-1: (Schumann 1988, BLM, https://doi.org/10.1007/BF00123019)
285 IF ((zmeas - zd_rsl)/l_mod_rsl < -neut_limit) ustar_rsl = max(0.15, ustar_rsl)
286
287 ! TStar_RSL = -1.*(qh/(avcp*avdens))/UStar_RSL
288 ! qStar_RSL = -1.*(qe/lv_J_kg*avdens)/UStar_RSL
289 IF (flag_rsl) THEN
290 ustar_heat = max(0.15, ustar_rsl)
291 ELSE
292 ! use UStar_heat implied by RA_h using MOST
293 psihz0 = stab_psi_heat(stabilitymethod, z0v/l_mod_rsl)
294 ustar_heat = 1/(kappa*ra_h)*(log((zmeas - zd_rsl)/z0v) - psihza + psihz0)
295 END IF
296 tstar_rsl = -1.*(qh/(avcp*avdens))/ustar_heat
297 IF (qe == 0.0) THEN
298 qstar_rsl = 10.**(-10) ! avoid the situation where qe=0, qstar_RSL=0 and the code breaks LB 21 May 2021
299 ELSE
300 qstar_rsl = -1.*(qe/lv_j_kg*avdens)/ustar_heat
301 END IF
302 qa_gkg = rh2qa(avrh/100, press_hpa, temp_c)
303
304 DO z = nz_can, nz
305 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
306 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
307 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0 + psihatm_z(z))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
308 dataoutlinetrsl(z) = (log((zarray(z) - zd_rsl)/(zmeas - zd_rsl)) - psihz + psihza + psihath_z(z) - psihath_z(nz))/kappa ! eqn. 4 in Theeuwes et al. (2019 BLM)
309 dataoutlineqrsl(z) = dataoutlinetrsl(z)
310 END DO
311 !
312 ! Step 7: calculate in canopy variables
313 !
314 IF (flag_rsl) THEN
315 ! RSL approach: exponential profiles within canopy
316 IF (nz_can > 1) THEN
317 t_h = scc*tstar_rsl/(beta*fx)
318 q_h = scc*qstar_rsl/(beta*fx)
319 DO z = 1, nz_can - 1
320 dataoutlineursl(z) = dataoutlineursl(nz_can)*exp(beta*(zarray(z) - zh_rsl)/elm)
321 dataoutlinetrsl(z) = dataoutlinetrsl(nz_can) + (t_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - t_h)/tstar_rsl
322 dataoutlineqrsl(z) = dataoutlineqrsl(nz_can) + (q_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - q_h)/qstar_rsl
323 END DO
324 END IF
325 ELSE
326 ! MOST approach:
327 DO z = 1, nz_can
328 ! when using MOST, all vertical levels should larger than zd_RSL
329 IF (zarray(z) <= zd_rsl) zarray(z) = 1.01*(zd_rsl + z0_rsl)
330 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
331 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
332 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0)/kappa
333 dataoutlinetrsl(z) = (log((zarray(z) - zd_rsl)/(zmeas - zd_rsl)) - psihz + psihza)/kappa
334 dataoutlineqrsl(z) = dataoutlinetrsl(z)
335 END DO
336 END IF
337
338 ! associate physical quantities to correction profilles
339 dataoutlineursl = dataoutlineursl*ustar_rsl
340 dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
341 dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
342
343 ! construct output line for output file
344 dataoutlinersl = [zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl, &
345 !information for debugging
346 ! L_stab, L_unstab,
347 l_mod_rsl, &
348 zh_rsl, &
349 ! Lc_stab, Lc_unstab,
350 lc, &
351 beta, zd_rsl, z0_rsl, elm, scc, fx, &
352 ustar_rsl, ustar_heat, tstar_rsl, fai, pai, merge(1.d0, 0.d0, flag_rsl) &
353 ]
354
355 !
356 ! Step 8
357 ! retrieve the diagnostics at key heights
358 !
359 IF (flag_rsl) THEN
360 ! RSL approach: diagnostics within canopy, heights are above ground level
361 t2_c = interp_z(2d0, zarray, dataoutlinetrsl)
362 q2_gkg = interp_z(2d0, zarray, dataoutlineqrsl)
363 u10_ms = interp_z(10d0, zarray, dataoutlineursl)
364 ELSE
365 ! MOST approach: diagnostics at heights above zdm+z0m to avoid insane values
366 t2_c = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlinetrsl)
367 q2_gkg = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlineqrsl)
368 u10_ms = interp_z(10d0 + zd_rsl + z0_rsl, zarray, dataoutlineursl)
369 END IF
370 ! get relative humidity:
371 rh2 = qa2rh(q2_gkg, press_hpa, t2_c)
372

References allocatearray::a3, interp_z(), atmmoiststab_module::neut_limit, meteo::qa2rh(), meteo::rh2qa(), rsl_cal_prms(), atmmoiststab_module::stab_psi_heat(), and atmmoiststab_module::stab_psi_mom().

Referenced by suews_driver::suews_cal_main().

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 9 of file suews_phys_rslprof.f95.

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