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)
 
subroutine rslprofile_dts (diagmethod, zh, z0m, zdm, z0v, l_mod, sfr_paved, sfr_bldg, sfr_evetr, sfr_dectr, sfr_grass, sfr_bsoil, sfr_water, 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 1515 of file suews_phys_rslprof.f95.

1516 ! TS, 03 Aug 2020:
1517 ! iterative determination of beta depending on Lc/L
1518 ! ref: eqn 10 & 11 in Harman (2012, BLM)
1519 IMPLICIT NONE
1520 INTEGER, INTENT(in) :: StabilityMethod
1521 REAL(KIND(1D0)), INTENT(in) :: beta0
1522 REAL(KIND(1D0)), INTENT(in) :: lc_over_l
1523 REAL(KIND(1D0)) :: beta_x
1524
1525 REAL(KIND(1D0)) :: phim, err, beta_x0
1526
1527 INTEGER :: it
1528
1529 it = 1
1530 phim = 1
1531 err = 1
1532 ! print *, '***********************'
1533 ! print *, 'beta0:', beta0
1534 ! print *, 'Lc/L_MOD:', lc_over_l
1535 DO WHILE ((err > 0.001) .AND. (it < 20))
1536 beta_x0 = beta0/phim
1537 phim = stab_phi_heat(stabilitymethod, (beta_x0**2)*lc_over_l)
1538 ! TODO: how to deal with neutral condition when phim=0? TS 05 Feb 2021
1539 ! here we use beta=0.35 as a temporary workaround but a better solution is still neded.
1540 beta_x = beta0/phim
1541 err = abs(beta_x - beta_x0)
1542 ! print *, it, err, beta_x0, beta_x, phim, lc_over_l
1543 it = it + 1
1544
1545 END DO
1546 ! print *, ''
1547

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

1465 ! Step 2:
1466 ! Parameterise beta according to Harman 2012 with upper limit of 0.5
1467 IMPLICIT NONE
1468
1469 INTEGER, INTENT(in) :: StabilityMethod ! stability method
1470 REAL(KIND(1D0)), INTENT(in) :: PAI
1471 REAL(KIND(1D0)), INTENT(in) :: sfr_tr
1472 REAL(KIND(1D0)), INTENT(in) :: lc_over_L
1473
1474 ! output
1475 REAL(KIND(1D0)) :: beta
1476
1477 ! internal use
1478 REAL(KIND(1D0)) :: betaHF
1479 REAL(KIND(1D0)) :: betaNL
1480
1481 REAL(KIND(1D0)), PARAMETER :: kappa = 0.4
1482 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
1483 ! real(KIND(1D0)) :: phim_hat
1484 ! real(KIND(1D0)) :: zd_RSL
1485
1486 REAL(KIND(1D0)) :: betaN2
1487 ! INTEGER :: it
1488 ! real(KIND(1D0)) :: err
1489 ! real(KIND(1D0)) :: phim
1490
1491 ! betaN for trees found to be 0.3 and for urban 0.4 linearly interpolate between the two using surface fractions
1492 ! betaN2 = 0.30 + (1.-sfr_surf(ConifSurf) - sfr_surf(ConifSurf))*0.1
1493 IF (pai > 0) THEN
1494 betan2 = 0.30*sfr_tr/pai + (pai - sfr_tr)/pai*0.4
1495 ELSE
1496 betan2 = 0.35
1497 END IF
1498
1499 betahf = cal_beta_lc(stabilitymethod, betan2, lc_over_l)
1500
1501 betanl = cal_beta_lc(stabilitymethod, kappa/2., lc_over_l)
1502
1503 IF (lc_over_l > a2) THEN
1504 beta = betahf
1505 ELSE
1506 beta = betanl + ((betahf - betanl)/(1.+a1*abs(lc_over_l - a2)**a3))
1507 END IF
1508
1509 IF (beta > 0.5) THEN
1510 beta = 0.5
1511 END IF
1512

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

1030
1031 IMPLICIT NONE
1032 INTEGER, INTENT(in) :: StabilityMethod ! stability method
1033 ! real(KIND(1D0)), intent(in) :: z ! height of interest [m]
1034 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
1035 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
1036 REAL(KIND(1D0)), INTENT(in) :: Scc !
1037 REAL(KIND(1D0)), INTENT(in) :: f !
1038 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
1039 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
1040 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
1041
1042 ! output
1043 REAL(KIND(1D0)), INTENT(out) :: ch
1044 REAL(KIND(1D0)), INTENT(out) :: c2h ! displacement height used in RSL
1045
1046 ! internal variables
1047 REAL(KIND(1D0)) :: phih_zh ! displacement height used in RSL
1048 REAL(KIND(1D0)) :: phih_zhdz ! displacement height used in RSL
1049 REAL(KIND(1D0)) :: dphih ! displacement height used in RSL
1050 REAL(KIND(1D0)) :: phi_hathZh ! displacement height used in RSL
1051
1052 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
1053 REAL(KIND(1D0)), PARAMETER :: dz = 0.1 !height step
1054
1055 phih_zh = stab_phi_heat(stabilitymethod, (zh_rsl - zd_rsl)/l_mod)
1056 phih_zhdz = stab_phi_heat(stabilitymethod, (zh_rsl - zd_rsl + 1.)/l_mod)
1057
1058 dphih = phih_zhdz - phih_zh
1059 IF (phih_zh /= 0.) THEN
1060 phi_hathzh = kappa*scc/(2.*beta*phih_zh)
1061 ELSE
1062 phi_hathzh = 1.
1063 END IF
1064
1065 IF (phi_hathzh >= 1.) THEN
1066 ! more stable, but less correct
1067 c2h = 0.5
1068 phi_hathzh = 1.
1069 ELSE
1070 ! if very unstable this might cause some high values of psihat_z
1071 c2h = (kappa*scc*(2.+f - (dphih*2.*beta**2.*lc/phih_zh)))/(2.*beta*phih_zh - kappa*scc)
1072 END IF
1073 ! force c2h to 0.5 for better stability. TS 14 Jul 2020
1074 ! TODO: a more proper threshold needs to be determined
1075 c2h = 0.5
1076
1077 ch = (1.-phi_hathzh)*exp(c2h/2.)
1078

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

977
978 IMPLICIT NONE
979 INTEGER, INTENT(in) :: StabilityMethod ! stability method
980 ! real(KIND(1D0)), intent(in) :: z ! height of interest [m]
981 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
982 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
983 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
984 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
985 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
986
987 ! output
988 REAL(KIND(1D0)), INTENT(out) :: c2
989 REAL(KIND(1D0)), INTENT(out) :: cm
990
991 ! internal variables
992 ! real(KIND(1D0)) ::phim_zh
993 REAL(KIND(1D0)) :: phi_hatmZh
994 REAL(KIND(1D0)) :: phim_zh
995 REAL(KIND(1D0)) :: phim_zhdz
996 REAL(KIND(1D0)) :: dphi
997 ! real(KIND(1D0)) ::phi_hatmZh
998
999 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
1000 REAL(KIND(1D0)), PARAMETER :: dz = 0.1 !height step
1001
1002 phim_zh = stab_phi_mom(stabilitymethod, (zh_rsl - zd_rsl)/l_mod)
1003 phim_zhdz = stab_phi_mom(stabilitymethod, (zh_rsl - zd_rsl + dz)/l_mod)
1004
1005 dphi = (phim_zhdz - phim_zh)/dz
1006 IF (phim_zh /= 0.) THEN
1007 phi_hatmzh = kappa/(2.*beta*phim_zh)
1008 ELSE
1009 ! neutral condition
1010 phi_hatmzh = 1.
1011 END IF
1012
1013 IF (phi_hatmzh >= 1.) THEN
1014 ! more stable, but less correct
1015 c2 = 0.5
1016 phi_hatmzh = 1.
1017 ELSE
1018 ! if very unstable this might cause some high values of psihat_z
1019 c2 = (kappa*(3.-(2.*beta**2.*lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
1020 END IF
1021 ! force c2 to 0.5 for better stability. TS 14 Jul 2020
1022 ! TODO: a more proper threshold needs to be determined
1023 c2 = 0.5
1024
1025 cm = (1.-phi_hatmzh)*exp(c2/2.)
1026

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

789
790 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
791 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
792
793 ! output
794 REAL(KIND(1D0)) :: elm ! a scaling parameter for RSL
795
796 elm = 2.*beta**3*lc
797

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

952 IMPLICIT NONE
953 INTEGER, INTENT(in) :: StabilityMethod ! stability method
954 REAL(KIND(1D0)), INTENT(in) :: z
955 REAL(KIND(1D0)), INTENT(in) :: zh_RSL
956 REAL(KIND(1D0)), INTENT(in) :: L_MOD
957 REAL(KIND(1D0)), INTENT(in) :: beta
958 REAL(KIND(1D0)), INTENT(in) :: lc
959 REAL(KIND(1D0)) :: phim_hat
960 REAL(KIND(1D0)) :: zd_RSL
961
962 REAL(KIND(1D0)) :: elm
963 REAL(KIND(1D0)) :: c2
964 REAL(KIND(1D0)) :: cm
965
966 elm = cal_elm_rsl(beta, lc)
967
968 zd_rsl = cal_zd_rsl(zh_rsl, beta, lc)
969
970 CALL cal_cm(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod, c2, cm)
971
972 phim_hat = 1 - cm*exp(-1.*c2*beta*(z - zd_rsl)/elm)
973

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

881 ! TS, 23 Oct 2019: calculate psi_hat for momentum
882 ! TS 30 Jun 2022: revised calculation logic for better calculation performance
883 IMPLICIT NONE
884 INTEGER, INTENT(in) :: StabilityMethod ! stability method
885 REAL(KIND(1D0)), INTENT(in) :: psihath_top ! psihath at z_top [m]
886 REAL(KIND(1D0)), INTENT(in) :: z_top ! height at a top level [m]
887 REAL(KIND(1D0)), INTENT(in) :: psihath_mid ! psihath at z_mid [m]
888 REAL(KIND(1D0)), INTENT(in) :: z_mid ! height at a middle level [m]
889 REAL(KIND(1D0)), INTENT(in) :: z_btm ! height of interest [m]
890 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
891 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! displacement height [m]
892 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
893 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
894 REAL(KIND(1D0)), INTENT(in) :: elm ! parameter in RSL
895 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
896 REAL(KIND(1D0)), INTENT(in) :: ch ! parameter in RSL
897 REAL(KIND(1D0)), INTENT(in) :: c2h ! parameter in RSL
898
899 ! output
900 REAL(KIND(1D0)) :: psihath_btm ! psim_hat at height of interest
901
902 ! internal variables
903 REAL(KIND(1D0)) :: phih_btm ! displacement height used in RSL
904 REAL(KIND(1D0)) :: phih_mid ! displacement height used in RSL
905 ! REAL(KIND(1D0)) :: psihatm_btm ! displacement height used in RSL
906 REAL(KIND(1D0)) :: slope_top ! displacement height used in RSL
907 REAL(KIND(1D0)) :: slope_btm ! displacement height used in RSL
908 REAL(KIND(1D0)) :: z_plus ! displacement height used in RSL
909 REAL(KIND(1D0)) :: psihath_plus ! displacement height used in RSL
910
911 REAL(KIND(1D0)), PARAMETER :: tol = 0.1 ! tolerance for iterative calculations
912 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
913 REAL(KIND(1D0)) :: dz_above
914
915 dz_above = z_mid - z_btm
916 ! single step calculation
917 phih_mid = stab_phi_heat(stabilitymethod, (z_mid - zd_rsl)/l_mod)
918 phih_btm = stab_phi_heat(stabilitymethod, (z_btm - zd_rsl)/l_mod)
919
920 psihath_btm = psihath_mid + dz_above/2.*phih_mid*(ch*exp(-1.*c2h*beta*(z_mid - zd_rsl)/elm)) & !Taylor's approximation for integral
921 /(z_mid - zd_rsl)
922 psihath_btm = psihath_btm + dz_above/2.*phih_btm*(ch*exp(-1.*c2h*beta*(z_btm - zd_rsl)/elm)) &
923 /(z_btm - zd_rsl)
924 IF (dz_above/zd_rsl < 1e-2) THEN
925 RETURN ! psihatm_z will be returned
926 END IF
927
928 IF (abs(psihath_btm) > 1e-3) THEN
929 ! test if recursion is needed by comparing slopes of psihat within the top and mid ranges
930 slope_top = (psihath_top - psihath_mid)/(z_top - z_mid)
931 slope_btm = (psihath_mid - psihath_btm)/(z_mid - z_btm)
932 IF (abs(slope_btm) > (1 + tol)*abs(slope_top)) THEN
933 z_plus = (z_mid + z_btm)/2
934 psihath_plus = cal_psih_hat( &
935 stabilitymethod, &
936 psihath_top, psihath_mid, &
937 z_top, z_mid, z_plus, &
938 ch, c2h, &
939 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
940 psihath_btm = cal_psih_hat( &
941 stabilitymethod, &
942 psihath_mid, psihath_plus, &
943 z_mid, z_plus, z_btm, &
944 ch, c2h, &
945 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
946 END IF
947 END IF
948

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

806 ! TS, 23 Oct 2019: calculate psi_hat for momentum
807 ! TS 30 Jun 2022: revised calculation logic for better calculation performance
808 IMPLICIT NONE
809 INTEGER, INTENT(in) :: StabilityMethod ! stability method
810 REAL(KIND(1D0)), INTENT(in) :: psihatm_top ! height of interest [m]
811 REAL(KIND(1D0)), INTENT(in) :: z_top ! height of interest [m]
812 REAL(KIND(1D0)), INTENT(in) :: psihatm_mid ! height of interest [m]
813 REAL(KIND(1D0)), INTENT(in) :: z_mid ! height of interest [m]
814 REAL(KIND(1D0)), INTENT(in) :: z_btm ! height of interest [m]
815 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
816 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! canyon depth [m]
817 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
818 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
819 REAL(KIND(1D0)), INTENT(in) :: elm ! parameter in RSL
820 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
821 REAL(KIND(1D0)), INTENT(in) :: cm ! Obukhov length [m]
822 REAL(KIND(1D0)), INTENT(in) :: c2 ! Obukhov length [m]
823
824 ! output
825 REAL(KIND(1D0)) :: psihatm_btm ! psim_hat at height of interest
826
827 ! internal variables
828 REAL(KIND(1D0)) :: phim_btm ! displacement height used in RSL
829 REAL(KIND(1D0)) :: phim_mid ! displacement height used in RSL
830 REAL(KIND(1D0)) :: slope_top ! displacement height used in RSL
831 REAL(KIND(1D0)) :: slope_btm ! displacement height used in RSL
832 REAL(KIND(1D0)) :: z_plus ! displacement height used in RSL
833 REAL(KIND(1D0)) :: psihatm_plus ! displacement height used in RSL
834
835 REAL(KIND(1D0)), PARAMETER :: tol = 0.5 ! tolerance for iterative calculations
836 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
837 REAL(KIND(1D0)) :: dz_above
838
839 dz_above = z_mid - z_btm
840 ! single step calculation
841 phim_mid = stab_phi_mom(stabilitymethod, (z_mid - zd_rsl)/l_mod)
842 phim_btm = stab_phi_mom(stabilitymethod, (z_btm - zd_rsl)/l_mod)
843
844 psihatm_btm = psihatm_mid + dz_above/2.*phim_mid*(cm*exp(-1.*c2*beta*(z_mid - zd_rsl)/elm)) & !Taylor's approximation for integral
845 /(z_mid - zd_rsl)
846 psihatm_btm = psihatm_btm + dz_above/2.*phim_btm*(cm*exp(-1.*c2*beta*(z_btm - zd_rsl)/elm)) &
847 /(z_btm - zd_rsl)
848 IF (dz_above/zd_rsl < 1e-2) THEN
849 RETURN ! psihatm_z will be returned
850 END IF
851
852 IF (abs(psihatm_btm) > 1e-3) THEN
853 ! test if recursion is needed by comparing slopes of psihat within the top and mid ranges
854 slope_top = (psihatm_top - psihatm_mid)/(z_top - z_mid)
855 slope_btm = (psihatm_mid - psihatm_btm)/(z_mid - z_btm)
856 IF (abs(slope_btm) > (1 + tol)*abs(slope_top)) THEN
857 z_plus = (z_mid + z_btm)/2
858 psihatm_plus = cal_psim_hat( &
859 stabilitymethod, &
860 psihatm_top, psihatm_mid, &
861 z_top, z_mid, z_plus, &
862 cm, c2, &
863 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
864 psihatm_btm = cal_psim_hat( &
865 stabilitymethod, &
866 psihatm_mid, psihatm_plus, &
867 z_mid, z_plus, z_btm, &
868 cm, c2, &
869 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
870 END IF
871 END IF
872

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

1252 ! calculate z0 iteratively
1253 ! TS, 23 Oct 2019
1254 IMPLICIT NONE
1255 INTEGER, INTENT(in) :: StabilityMethod
1256 REAL(KIND(1D0)), INTENT(in) :: zH_RSL ! canyon depth [m]
1257 REAL(KIND(1D0)), INTENT(in) :: zd_RSL ! displacement height [m]
1258 REAL(KIND(1D0)), INTENT(in) :: L_MOD_RSL ! Monin Obukhov length[m]
1259 ! REAL(KIND(1D0)), INTENT(in) :: Lc ! canyon length scale [m]
1260 REAL(KIND(1D0)), INTENT(in) :: beta ! height scale for bluff bodies [m]
1261 REAL(KIND(1D0)), INTENT(in) :: psihatm_Zh ! psihatm at zH
1262
1263 ! output
1264 REAL(KIND(1D0)) :: z0_RSL
1265
1266 ! internal variables
1267 REAL(KIND(1D0)) :: psimZh, psimz0, z0_RSL_x
1268 REAL(KIND(1D0)) :: err
1269 INTEGER :: it
1270
1271 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40
1272 ! REAL(KIND(1d0)), PARAMETER::r = 0.1
1273 ! REAL(KIND(1d0)), PARAMETER::a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1.
1274
1275 psimzh = stab_psi_mom(stabilitymethod, (zh_rsl - zd_rsl)/l_mod_rsl)
1276 ! psihatm_Zh = cal_psim_hat(StabilityMethod, Zh_RSL, zh_RSL, L_MOD_RSL, beta, Lc)
1277
1278 !first guess
1279 z0_rsl = 0.1*zh_rsl
1280 err = 10.
1281 ! psimz0 = 0.5
1282 it = 1
1283 DO WHILE ((err > 0.001) .AND. (it < 10))
1284 z0_rsl_x = z0_rsl
1285 psimz0 = stab_psi_mom(stabilitymethod, z0_rsl_x/l_mod_rsl)
1286 z0_rsl = (zh_rsl - zd_rsl)*exp(-1.*kappa/beta)*exp(-1.*psimzh + psimz0)*exp(psihatm_zh)
1287 err = abs(z0_rsl_x - z0_rsl)
1288 it = it + 1
1289 END DO
1290
1291 ! set limit on z0_RSL for numeric stability
1292 z0_rsl = merge(z0_rsl, 1d-2, z0_rsl > 1d-2)
1293

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

1237
1238 REAL(KIND(1D0)), INTENT(in) :: zh_RSL ! canyon depth [m]
1239 REAL(KIND(1D0)), INTENT(in) :: Lc ! height scale for bluff bodies [m]
1240 REAL(KIND(1D0)), INTENT(in) :: beta ! parameter in RSL
1241
1242 ! output
1243 REAL(KIND(1D0)) :: zd_RSL ! zd used in RSL
1244
1245 zd_rsl = zh_rsl - (beta**2.)*lc
1246 !correct negative values using rule of thumb, TS 24 Jun 2020
1247 ! if (zd_RSL < 0) zd_RSL = 0.7*Zh_RSL
1248

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

750
751 REAL(KIND(1D0)), INTENT(in) :: z_x ! height to interpolate at
752 REAL(KIND(1D0)), DIMENSION(nz), INTENT(in) :: z ! heights
753 REAL(KIND(1D0)), DIMENSION(nz), INTENT(in) :: v ! values associated with heights
754
755 ! output
756 REAL(KIND(1D0)) :: v_x ! zd used in RSL
757
758 ! local variables
759 REAL(KIND(1D0)) :: slope ! slope
760 REAL(KIND(1D0)) :: dz ! slope
761 REAL(KIND(1D0)), DIMENSION(nz) :: dif ! slope
762 INTEGER :: idx_low ! vertical index lower than z_x
763 INTEGER :: idx_x ! vertical index lower than z_x
764 INTEGER :: idx_high ! vertical index higher than z_x
765 ! INTEGER :: idx ! vertical index higher than z_x
766 INTEGER, PARAMETER :: nz = 30 ! vertical index higher than z_x
767
768 ! initialise variables
769 idx_x = 0
770
771 dif = z - z_x
772 idx_x = maxloc(dif, 1, abs(dif) < 1.d-6)
773 idx_low = maxloc(dif, 1, dif < 0.)
774 idx_high = minloc(dif, 1, dif > 0.)
775
776 IF (idx_x > 0) THEN
777 ! z_x is one of zarray elements
778 v_x = v(idx_x)
779 ELSE
780 ! linear interpolation is performed
781 dz = z(idx_high) - z(idx_low)
782 slope = (v(idx_high) - v(idx_low))/dz
783 v_x = v(idx_low) + (z_x - z(idx_low))*slope
784 END IF
785

Referenced by rslprofile(), and rslprofile_dts().

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

1299
1300 IMPLICIT NONE
1301 INTEGER, INTENT(in) :: StabilityMethod ! stability method
1302 INTEGER, INTENT(IN) :: nz_above ! number of layers above zh
1303 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(in) :: z_array ! land cover fractions
1304 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(out) :: psihatm_array ! land cover fractions
1305 REAL(KIND(1D0)), DIMENSION(nz_above), INTENT(out) :: psihath_array ! land cover fractions
1306 REAL(KIND(1D0)), INTENT(in) :: zh ! canyon depth [m]
1307 REAL(KIND(1D0)), INTENT(in) :: FAI ! frontal area index inlcuding trees
1308 ! REAL(KIND(1D0)), INTENT(in) :: FAIBldg ! frontal area index of buildings
1309 REAL(KIND(1D0)), INTENT(in) :: PAI ! plan area index inlcuding area of trees
1310 ! REAL(KIND(1D0)), INTENT(in) :: porosity_dectr ! porosity of deciduous trees
1311 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
1312 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf ! land cover fractions
1313
1314 ! output
1315 ! real(KIND(1D0)), intent(out) ::L_stab ! threshold for Obukhov length under stable conditions
1316 ! real(KIND(1D0)), intent(out) ::L_unstab ! threshold for Obukhov length under unstable conditions
1317 REAL(KIND(1D0)), INTENT(out) :: L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
1318 REAL(KIND(1D0)), INTENT(out) :: zH_RSL ! mean canyon height used in RSL module with thresholds applied
1319 ! real(KIND(1D0)), intent(out) ::Lc_stab ! threshold for penetration distance scale under stable conditions
1320 ! real(KIND(1D0)), intent(out) ::Lc_unstab ! threshold for penetration distance scale under unstable conditions
1321 REAL(KIND(1D0)), INTENT(out) :: Lc ! penetration distance scale for bluff bodies [m]
1322 REAL(KIND(1D0)), INTENT(out) :: beta ! psim_hat at height of interest
1323 REAL(KIND(1D0)), INTENT(out) :: zd_RSL ! displacement height to prescribe if necessary [m]
1324 REAL(KIND(1D0)), INTENT(out) :: z0_RSL ! roughness length [m]
1325 REAL(KIND(1D0)), INTENT(out) :: elm ! length scale used in RSL
1326 REAL(KIND(1D0)), INTENT(out) :: Scc ! parameter in RSL
1327 REAL(KIND(1D0)), INTENT(out) :: fx ! parameter in RSL
1328
1329 ! internal variables
1330 INTEGER :: iz
1331 REAL(KIND(1D0)) :: sfr_tr
1332 REAL(KIND(1D0)) :: psihatm_Zh
1333 ! real(KIND(1D0)) ::L_MOD_RSL_x
1334 ! real(KIND(1D0)) ::lc_x
1335 REAL(KIND(1D0)) :: lc_over_L
1336 REAL(KIND(1D0)) :: z_top, z_mid, z_btm
1337 REAL(KIND(1D0)) :: psihatm_top, psihatm_mid, psihatm_btm
1338 REAL(KIND(1D0)) :: psihath_top, psihath_mid, psihath_btm
1339 REAL(KIND(1D0)) :: cm, ch, c2m, c2h
1340 ! real(KIND(1D0)) ::betaHF
1341 ! real(KIND(1D0)) ::betaNL
1342 REAL(KIND(1D0)) :: Lc_min ! LB Oct2021 - minimum value of Lc
1343 REAL(KIND(1D0)) :: dim_bluffbody ! LB Oct2021 - horizontal building dimensions
1344
1345 REAL(KIND(1D0)), PARAMETER :: cd_tree = 1.2 ! drag coefficient tree canopy !!!!needs adjusting!!!
1346 REAL(KIND(1D0)), PARAMETER :: a_tree = 0.05 ! the foliage area per unit volume !!!!needs adjusting!!!
1347 ! lv_J_kg = 2.5E6, &! latent heat for water vapor!!! make consistant with rest of code
1348 REAL(KIND(1D0)), PARAMETER :: beta_N = 0.40 ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
1349 REAL(KIND(1D0)), PARAMETER :: pi = 4.*atan(1.0)
1350
1351 REAL(KIND(1D0)), PARAMETER :: planF_low = 1e-6
1352 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40 ! von karman constant
1353 ! REAL(KIND(1d0)), PARAMETER::z0m= 0.40
1354 REAL(KIND(1D0)), PARAMETER :: r = 0.1
1355 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
1356 REAL(KIND(1D0)), PARAMETER :: Zh_min = 0.4 ! limit for minimum canyon height used in RSL module
1357 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
1358
1359 ! under stable conditions, set a threshold for L_MOD to avoid numerical issues. TS 28 Oct 2019
1360 ! L_MOD = merge(L_MOD, 300.d1, L_MOD < 300.)
1361
1362 ! zH_RSL
1363 zh_rsl = max(zh, zh_min)
1364
1365 ! ! land cover fraction of bluff bodies
1366 ! PAI = DOT_PRODUCT(sfr_surf([BldgSurf, ConifSurf, DecidSurf]), [1D0, 1 - porosity_evetr, 1 - porosity_dectr])
1367 ! set a threshold for sfr_zh to avoid numerical difficulties
1368 ! PAI = min(PAI, 0.8)
1369
1370 ! land cover fraction of trees
1371 sfr_tr = sum(sfr_surf([conifsurf, decidsurf]))
1372
1373 ! height scale for buildings !not used? why?
1374 ! Lc_build = (1.-sfr_surf(BldgSurf))/FAI*Zh_RSL ! Coceal and Belcher 2004 assuming Cd = 2
1375
1376 ! height scale for tress
1377 ! Lc_tree = 1./(cd_tree*a_tree) ! not used? why?
1378
1379 ! height scale for bluff bodies
1380 ! LB Oct2021 - replaced PAI with sfr(BldgSurf) since the parameter should represent the solid fraction (and trees have negligible solid fraction).
1381 ! TS 18 Jun 2022 - changed sfr_surf(BldgSurf) back to PAI to account for trees with PAI updated by accounting for porosity.
1382 lc = (1.-pai)/fai*zh_rsl
1383
1384 ! set a minimum threshold (of 0.5*Zh_RSL) for Lc to avoid numerical diffulties when FAI is too large (e.g., FAI>10)
1385 ! Lc = MERGE(Lc, 0.5*Zh_RSL, Lc > 0.5*Zh_RSL)
1386 ! 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
1387 ! 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
1388 dim_bluffbody = zh_rsl*pai/fai
1389 lc_min = dim_bluffbody*pai**(-0.5)/3.
1390 lc = merge(lc, lc_min, lc > lc_min)
1391
1392 ! a normalised scale with a physcially valid range between [-2,2] (Harman 2012, BLM)
1393 lc_over_l = lc/l_mod
1394 ! lc_over_L = lc/L_MOD_RSL_x
1395 IF (lc_over_l > 0) THEN
1396 lc_over_l = min(2., lc_over_l)
1397 ELSE
1398 lc_over_l = max(-2., lc_over_l)
1399 END IF
1400 ! correct L_MOD_RSL
1401 l_mod_rsl = lc/lc_over_l
1402
1403 ! Step 2:
1404 ! Parameterise beta according to Harman 2012 with upper limit of 0.5
1405 beta = cal_beta_rsl(stabilitymethod, pai, sfr_tr, lc_over_l)
1406
1407 ! Schmidt number Harman and Finnigan 2008: assuming the same for heat and momemntum
1408 scc = 0.5 + 0.3*tanh(2.*lc_over_l)
1409 fx = 0.5*((1.+4.*r*scc)**0.5) - 0.5
1410
1411 ! zero displacement height used in RSL
1412 zd_rsl = cal_zd_rsl(zh_rsl, beta, lc)
1413
1414 ! scaling parameter for RSL
1415 elm = cal_elm_rsl(beta, lc)
1416
1417 CALL cal_ch(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod_rsl, scc, fx, c2h, ch)
1418 CALL cal_cm(stabilitymethod, zh_rsl, zd_rsl, lc, beta, l_mod_rsl, c2m, cm)
1419
1420 ! calculate psihat values at desirable heights
1421 psihatm_top = 0
1422 psihatm_mid = 0
1423 psihatm_array(nz_above) = 0
1424 psihatm_array(nz_above - 1) = 0
1425
1426 psihath_top = 0
1427 psihath_mid = 0
1428 psihath_array(nz_above) = 0
1429 psihath_array(nz_above - 1) = 0
1430
1431 DO iz = nz_above, 3, -1
1432 z_top = z_array(iz)
1433 z_mid = z_array(iz - 1)
1434 z_btm = z_array(iz - 2)
1435
1436 ! momentum
1437 psihatm_btm = cal_psim_hat(stabilitymethod, &
1438 psihatm_top, psihatm_mid, &
1439 z_top, z_mid, z_btm, &
1440 cm, c2m, &
1441 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
1442 psihatm_array(iz - 2) = psihatm_btm
1443 psihatm_top = psihatm_mid
1444 psihatm_mid = psihatm_btm
1445
1446 ! heat
1447 psihath_btm = cal_psih_hat(stabilitymethod, &
1448 psihath_top, psihath_mid, &
1449 z_top, z_mid, z_btm, &
1450 ch, c2h, &
1451 zh_rsl, zd_rsl, l_mod, beta, elm, lc)
1452 psihath_top = psihath_mid
1453 psihath_mid = psihath_btm
1454 psihath_array(iz - 2) = psihath_btm
1455
1456 END DO
1457
1458 ! calculate z0 iteratively
1459 psihatm_zh = psihatm_array(1)
1460 z0_rsl = cal_z0_rsl(stabilitymethod, zh_rsl, zd_rsl, beta, l_mod_rsl, psihatm_zh)
1461

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(), and rslprofile_dts().

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 ! FAI < 0.45 .AND. & ! remove the lower limit of FAI
161 ! PAI > 0.1 .AND. PAI < 0.61 .AND. & ! PAI
162 pai > 0.1 .AND. pai < 0.68 .AND. & ! PAI
163 zh > 2 ! effective canopy height
164 ! LB Oct2021 - FAI and PAI can be larger than 0.45 and 0.61 respectively -> remove (1.-PAI)/FAI > .021 constraint
165 !(note: it seems wrong anyway - should be 0.87 not 0.021 based on G&O1991 numbers)
166 ELSE
167 ! default is to use MOST
168 flag_rsl = .false.
169 END IF
170
171 !
172 ! ! Step 2
173 ! ! determine vertical levels used in RSL
174 ! Define the height array with consideration of key heights
175 ! set number of heights within canopy
176 IF (zh <= 2) THEN
177 nz_can = 3
178 ELSE IF (zh <= 10) THEN
179 nz_can = 6
180 ELSE
181 nz_can = 15
182 END IF
183 ! fill up heights in canopy
184 dz_can = zh/nz_can
185 DO i = 1, nz_can
186 zarray(i) = dz_can*i
187 END DO
188
189 ! guaranttee 2 m is within the zarray
190 IF (dz_can > 2) zarray(1) = 1.999
191
192 ! fill up heights above canopy
193 nz_above = nz - nz_can
194 dz_above = (zmeas - zh)/nz_above
195 DO i = nz_can + 1, nz
196 zarray(i) = zh + (i - nz_can)*dz_above
197 END DO
198
199 ! ! determine index at the canyon top
200 ! DO z = 1, nz
201 ! dif(z) = ABS(zarray(z) - Zh)
202 ! END DO
203 ! nz_can = MINLOC(dif, DIM=1)
204 ! ! zarray(idx_can+2) = Zh+.1
205 ! ! zarray(idx_can+1) = Zh+.05
206 ! zarray(nz_can) = Zh
207 ! zarray(idx_can-1) = Zh-.1
208
209 ! determine index at measurement height
210 ! nz = nz
211 zarray(nz) = zmeas
212
213 IF (flag_rsl) THEN
214
215 ! use RSL approach to calculate correction factors
216 psihatm_z = 0.*zarray
217 psihath_z = 0.*zarray
218 ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
219 CALL rsl_cal_prms( &
220 stabilitymethod, & !input
221 nz_above, zarray(nz_can + 1:nz), & !input
222 zh, l_mod, sfr_surf, fai, pai, & !input
223 psihatm_z(nz_can + 1:nz), psihath_z(nz_can + 1:nz), & !output
224 zh_rsl, l_mod_rsl, & ! output
225 lc, beta, zd_rsl, z0_rsl, elm, scc, fx)
226
227 ! Step 3: calculate the stability dependent H&F constants
228
229 ! CALL cal_ch(StabilityMethod, zh_RSL, zd_RSL, Lc, beta, L_MOD_RSL, Scc, fx, c2h, ch)
230 ! CALL cal_cm(StabilityMethod, zH_RSL, zd_RSL, Lc, beta, L_MOD_RSL, c2m, cm)
231
232 ! ! Step 4: determine psihat at levels above the canopy
233 ! DO z = nz - 1, idx_can, -1
234 ! phimz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
235 ! phimzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
236 ! phihz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
237 ! phihzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
238
239 ! 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
240 ! /(zarray(z + 1) - zd_RSL)
241 ! psihatm_z(z) = psihatm_z(z) + dz_above/2.*phimz*(cm*EXP(-1.*c2m*beta*(zarray(z) - zd_RSL)/elm)) &
242 ! /(zarray(z) - zd_RSL)
243 ! 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
244 ! /(zarray(z + 1) - zd_RSL)
245 ! psihath_z(z) = psihath_z(z) + dz_above/2.*phihz*(ch*EXP(-1.*c2h*beta*(zarray(z) - zd_RSL)/elm)) &
246 ! /(zarray(z) - zd_RSL)
247 ! END DO
248
249 ELSE
250
251 ! correct parameters if RSL approach doesn't apply for scenario of isolated flows
252 ! see Fig 1 of Grimmond and Oke (1999)
253 ! when isolated flow or skimming flow, implying RSL doesn't apply, force RSL correction to zero
254 psihatm_z = 0
255 psihath_z = 0
256
257 ! use L_MOD as in other parts of SUEWS
258 l_mod_rsl = l_mod
259
260 !correct RSL-based using SUEWS system-wide values
261 z0_rsl = z0m
262 zd_rsl = zdm
263 zh_rsl = zh
264
265 lc = -999
266 beta = -999
267 scc = -999
268 fx = -999
269 elm = -999
270
271 ! then MOST recovers from RSL correction
272 END IF
273
274 ! Step 6: Calculate mean variables above canopy
275 !
276 psimz0 = stab_psi_mom(stabilitymethod, z0_rsl/l_mod_rsl)
277 psimza = stab_psi_mom(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
278 psihza = stab_psi_heat(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
279
280 ustar_rsl = avu1*kappa/(log((zmeas - zd_rsl)/z0_rsl) - psimza + psimz0 + psihatm_z(nz))
281
282 ! TS 11 Feb 2021: limit UStar and TStar to reasonable ranges
283 ! 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
284 ustar_rsl = max(0.001, ustar_rsl)
285 ! under convective/unstable condition, min(UStar)==0.15 m s-1: (Schumann 1988, BLM, https://doi.org/10.1007/BF00123019)
286 IF ((zmeas - zd_rsl)/l_mod_rsl < -neut_limit) ustar_rsl = max(0.15, ustar_rsl)
287
288 ! TStar_RSL = -1.*(qh/(avcp*avdens))/UStar_RSL
289 ! qStar_RSL = -1.*(qe/lv_J_kg*avdens)/UStar_RSL
290 IF (flag_rsl) THEN
291 ustar_heat = max(0.15, ustar_rsl)
292 ELSE
293 ! use UStar_heat implied by RA_h using MOST
294 psihz0 = stab_psi_heat(stabilitymethod, z0v/l_mod_rsl)
295 ustar_heat = 1/(kappa*ra_h)*(log((zmeas - zd_rsl)/z0v) - psihza + psihz0)
296 END IF
297 tstar_rsl = -1.*(qh/(avcp*avdens))/ustar_heat
298 IF (qe == 0.0) THEN
299 qstar_rsl = 10.**(-10) ! avoid the situation where qe=0, qstar_RSL=0 and the code breaks LB 21 May 2021
300 ELSE
301 qstar_rsl = -1.*(qe/lv_j_kg*avdens)/ustar_heat
302 END IF
303 qa_gkg = rh2qa(avrh/100, press_hpa, temp_c)
304
305 DO z = nz_can, nz
306 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
307 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
308 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0 + psihatm_z(z))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
309 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)
310 dataoutlineqrsl(z) = dataoutlinetrsl(z)
311 END DO
312 !
313 ! Step 7: calculate in canopy variables
314 !
315 IF (flag_rsl) THEN
316 ! RSL approach: exponential profiles within canopy
317 IF (nz_can > 1) THEN
318 t_h = scc*tstar_rsl/(beta*fx)
319 q_h = scc*qstar_rsl/(beta*fx)
320 DO z = 1, nz_can - 1
321 dataoutlineursl(z) = dataoutlineursl(nz_can)*exp(beta*(zarray(z) - zh_rsl)/elm)
322 dataoutlinetrsl(z) = dataoutlinetrsl(nz_can) + (t_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - t_h)/tstar_rsl
323 dataoutlineqrsl(z) = dataoutlineqrsl(nz_can) + (q_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - q_h)/qstar_rsl
324 END DO
325 END IF
326 ELSE
327 ! MOST approach:
328 DO z = 1, nz_can
329 ! when using MOST, all vertical levels should larger than zd_RSL
330 IF (zarray(z) <= zd_rsl) zarray(z) = 1.01*(zd_rsl + z0_rsl)
331 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
332 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
333 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0)/kappa
334 dataoutlinetrsl(z) = (log((zarray(z) - zd_rsl)/(zmeas - zd_rsl)) - psihz + psihza)/kappa
335 dataoutlineqrsl(z) = dataoutlinetrsl(z)
336 END DO
337 END IF
338
339 ! associate physical quantities to correction profilles
340 dataoutlineursl = dataoutlineursl*ustar_rsl
341 dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
342 dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
343
344 ! construct output line for output file
345 dataoutlinersl = [zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl, &
346 !information for debugging
347 ! L_stab, L_unstab,
348 l_mod_rsl, &
349 zh_rsl, &
350 ! Lc_stab, Lc_unstab,
351 lc, &
352 beta, zd_rsl, z0_rsl, elm, scc, fx, &
353 ustar_rsl, ustar_heat, tstar_rsl, fai, pai, merge(1.d0, 0.d0, flag_rsl) &
354 ]
355
356 !
357 ! Step 8
358 ! retrieve the diagnostics at key heights
359 !
360 IF (flag_rsl) THEN
361 ! RSL approach: diagnostics within canopy, heights are above ground level
362 t2_c = interp_z(2d0, zarray, dataoutlinetrsl)
363 q2_gkg = interp_z(2d0, zarray, dataoutlineqrsl)
364 u10_ms = interp_z(10d0, zarray, dataoutlineursl)
365 ELSE
366 ! MOST approach: diagnostics at heights above zdm+z0m to avoid insane values
367 t2_c = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlinetrsl)
368 q2_gkg = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlineqrsl)
369 u10_ms = interp_z(10d0 + zd_rsl + z0_rsl, zarray, dataoutlineursl)
370 END IF
371 ! get relative humidity:
372 rh2 = qa2rh(q2_gkg, press_hpa, t2_c)
373

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:

◆ rslprofile_dts()

subroutine rsl_module::rslprofile_dts ( 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)), intent(in) sfr_paved,
real(kind(1d0)), intent(in) sfr_bldg,
real(kind(1d0)), intent(in) sfr_evetr,
real(kind(1d0)), intent(in) sfr_dectr,
real(kind(1d0)), intent(in) sfr_grass,
real(kind(1d0)), intent(in) sfr_bsoil,
real(kind(1d0)), intent(in) sfr_water,
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 376 of file suews_phys_rslprof.f95.

387 !-----------------------------------------------------
388 ! calculates windprofiles using MOST with a RSL-correction
389 ! based on Harman & Finnigan 2007
390 !
391 ! last modified by:
392 ! NT 16 Mar 2019: initial version
393 ! TS 16 Oct 2019: improved consistency in parameters/varaibles within SUEWS
394 ! TODO how to improve the speed of this code
395 !
396 !-----------------------------------------------------
397
398 IMPLICIT NONE
399
400 REAL(KIND(1D0)), INTENT(IN) :: sfr_paved
401 REAL(KIND(1D0)), INTENT(IN) :: sfr_bldg
402 REAL(KIND(1D0)), INTENT(IN) :: sfr_evetr
403 REAL(KIND(1D0)), INTENT(IN) :: sfr_dectr
404 REAL(KIND(1D0)), INTENT(IN) :: sfr_grass
405 REAL(KIND(1D0)), INTENT(IN) :: sfr_bsoil
406 REAL(KIND(1D0)), INTENT(IN) :: sfr_water
407 REAL(KIND(1D0)), DIMENSION(nsurf) :: sfr_surf !surface fraction ratio [-]
408
409 REAL(KIND(1D0)), INTENT(in) :: zMeas ! height of atmospheric forcing [m]
410 REAL(KIND(1D0)), INTENT(in) :: avU1 ! Wind speed at forcing height [m s-1]
411 REAL(KIND(1D0)), INTENT(in) :: Temp_C ! Air temperature at forcing height [C]
412 REAL(KIND(1D0)), INTENT(in) :: avRH ! relative humidity at forcing height [-]
413 REAL(KIND(1D0)), INTENT(in) :: Press_hPa ! pressure at forcing height [hPa]
414 REAL(KIND(1D0)), INTENT(in) :: L_MOD ! Obukhov length [m]
415 REAL(KIND(1D0)), INTENT(in) :: RA_h ! aerodynamic resistance for heat [s m-1]
416 REAL(KIND(1D0)), INTENT(in) :: avcp ! specific heat capacity [J kg-1 K-1]
417 REAL(KIND(1D0)), INTENT(in) :: lv_J_kg ! Latent heat of vaporization in [J kg-1]
418 REAL(KIND(1D0)), INTENT(in) :: avdens ! air density [kg m-3]
419 REAL(KIND(1D0)), INTENT(in) :: qh ! sensible heat flux [W m-2]
420 REAL(KIND(1D0)), INTENT(in) :: qe ! Latent heat flux [W m-2]
421 REAL(KIND(1D0)), INTENT(in) :: Zh ! Mean building height [m]
422 REAL(KIND(1D0)), INTENT(in) :: z0m ! roughness for momentum [m]
423 REAL(KIND(1D0)), INTENT(in) :: z0v ! roughnesslength for heat [s m-1]
424 REAL(KIND(1D0)), INTENT(in) :: zdm ! zero-plane displacement [m]
425 REAL(KIND(1D0)), INTENT(in) :: FAI ! Frontal area index [-]
426 REAL(KIND(1D0)), INTENT(in) :: PAI ! Plan area index [-]
427 ! REAL(KIND(1D0)), INTENT(in) :: FAIBldg ! Frontal area index of buildings [-]
428 ! REAL(KIND(1D0)), INTENT(in) :: porosity_dectr ! porosity of deciduous trees [-]
429
430 INTEGER, INTENT(in) :: StabilityMethod
431 INTEGER, INTENT(in) :: DiagMethod
432
433 REAL(KIND(1D0)), INTENT(out) :: T2_C ! Air temperature at 2 m [C]
434 REAL(KIND(1D0)), INTENT(out) :: q2_gkg ! Air specific humidity at 2 m [g kg-1]
435 REAL(KIND(1D0)), INTENT(out) :: U10_ms ! wind speed at 10 m [m s-1]
436 REAL(KIND(1D0)), INTENT(out) :: RH2 ! Air relative humidity [-]
437
438 INTEGER, PARAMETER :: nz = 30 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
439
440 REAL(KIND(1D0)), PARAMETER :: cd_tree = 1.2 ! drag coefficient tree canopy !!!!needs adjusting!!!
441 REAL(KIND(1D0)), PARAMETER :: a_tree = 0.05 ! the foliage area per unit volume !!!!needs adjusting!!!
442 REAL(KIND(1D0)), PARAMETER :: kappa = 0.40 ! von karman constant
443
444 REAL(KIND(1D0)), PARAMETER :: beta_N = 0.40 ! H&F beta coefficient in neutral conditions from Theeuwes et al., 2019 BLM
445 REAL(KIND(1D0)), PARAMETER :: pi = 4.*atan(1.0), r = 0.1
446 REAL(KIND(1D0)), PARAMETER :: a1 = 4., a2 = -0.1, a3 = 1.5, a4 = -1. ! constraints to determine beta
447
448 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
449
450 ! Variables array [z,U,T,q, 12 debug vars]
451 ! z: height array
452 ! U,T,q: wind speed, air temp, specific humidity at z;
453 ! debug vars: see dataoutLineRSL
454 REAL(KIND(1D0)), INTENT(out), DIMENSION(ncolumnsDataOutRSL - 5) :: dataoutLineRSL
455 REAL(KIND(1D0)), DIMENSION(nz) :: psihatm_z
456 REAL(KIND(1D0)), DIMENSION(nz) :: psihath_z
457 ! REAL(KIND(1D0)), DIMENSION(nz) :: dif
458 ! REAL(KIND(1d0)), DIMENSION(nz):: psihatm_z, psihath_z
459 REAL(KIND(1D0)), DIMENSION(nz) :: zarray
460 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineURSL ! wind speed array [m s-1]
461 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineTRSL ! Temperature array [C]
462 REAL(KIND(1D0)), DIMENSION(nz) :: dataoutLineqRSL ! Specific humidity array [g kg-1]
463
464 REAL(KIND(1D0)) :: z0_RSL ! roughness length from H&F
465 REAL(KIND(1D0)) :: zd_RSL ! zero-plane displacement
466
467 ! REAL(KIND(1d0))::Lc_build, Lc_tree, Lc ! canopy drag length scale
468 REAL(KIND(1D0)) :: Lc ! canopy drag length scale
469 ! REAL(KIND(1d0))::Lc_stab ! threshold of canopy drag length scale under stable conditions
470 ! REAL(KIND(1d0))::Lc_unstab ! threshold of canopy drag length scale under unstable conditions
471 REAL(KIND(1D0)) :: Scc ! Schmidt number for temperature and humidity
472 REAL(KIND(1D0)) :: psimz, psimz0, psimza, psihz, psihz0, psihza ! stability function for momentum
473 ! REAL(KIND(1d0))::betaHF, betaNL, beta, betaN2 ! beta coefficient from Harman 2012
474 REAL(KIND(1D0)) :: beta ! beta coefficient from Harman 2012
475 REAL(KIND(1D0)) :: elm ! mixing length
476 ! REAL(KIND(1d0))::xxm1, xxm1_2, xxh1, xxh1_2, dphi, dphih ! dummy variables for stability functions
477 REAL(KIND(1D0)) :: fx ! H&F'07 and H&F'08 'constants'
478 REAL(KIND(1D0)) :: t_h, q_h ! H&F'08 canopy corrections
479 REAL(KIND(1D0)) :: TStar_RSL ! temperature scale
480 REAL(KIND(1D0)) :: UStar_RSL ! friction velocity used in RSL
481 REAL(KIND(1D0)) :: UStar_heat ! friction velocity derived from RA_h with correction/restriction
482 ! REAL(KIND(1D0)) :: PAI ! plan area index, including areas of roughness elements: buildings and trees
483 ! REAL(KIND(1d0))::sfr_tr ! land cover fraction of trees
484 REAL(KIND(1D0)) :: L_MOD_RSL ! Obukhov length used in RSL module with thresholds applied
485 ! real(KIND(1D0))::L_stab ! threshold for Obukhov length under stable conditions
486 ! real(KIND(1D0))::L_unstab ! threshold for Obukhov length under unstable conditions
487
488 REAL(KIND(1D0)) :: zH_RSL ! mean canyon height used in RSL module with thresholds applied
489 REAL(KIND(1D0)) :: dz_above ! height step above canopy
490 REAL(KIND(1D0)) :: dz_can ! height step within canopy
491 ! REAL(KIND(1D0)) :: phi_hatmZh, phim_zh
492 ! REAL(KIND(1d0)), parameter::zH_min = 8! limit for minimum canyon height used in RSL module
493 ! REAL(KIND(1D0)), PARAMETER :: ratio_dz = 1.618 ! ratio between neighbouring height steps
494
495 REAL(KIND(1D0)) :: qa_gkg, qStar_RSL ! specific humidity scale
496 INTEGER :: I, z
497 INTEGER :: nz_can ! number of heights in canyon
498 INTEGER :: nz_above ! number of heights above canyon
499
500 LOGICAL :: flag_RSL ! whether RSL correction is used
501
502 ! CHARACTER(len=1024) :: Errmessage
503 ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
504 ! Step 1: determine if RSL should be used
505 ! Step 2: determine vertical levels used in RSL
506 ! Step 3: calculate the stability dependent H&F constants
507 ! Step 4: determine psihat at levels above the canopy
508 ! Step 5: Calculate z0 iteratively
509 ! Step 6: Calculate mean variables above canopy
510 ! Step 7: Calculate mean variables in canopy
511
512! ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
513! CALL RSL_cal_prms( &
514! StabilityMethod, & !input
515! zh, L_MOD, sfr_surf, FAI, PAI, & !input
516! zH_RSL, L_MOD_RSL, & ! output
517! Lc, beta, zd_RSL, z0_RSL, elm, Scc, fx)
518
519 ! Step 1: determine if RSL should be used
520 sfr_surf = [sfr_paved, sfr_bldg, sfr_evetr, sfr_dectr, sfr_grass, sfr_bsoil, sfr_water]
521
522 IF (diagmethod == 0) THEN
523 ! force MOST to be used
524 flag_rsl = .false.
525 ELSEIF (diagmethod == 1) THEN
526 ! force RSL to be used
527 flag_rsl = .true.
528 ELSEIF (diagmethod == 2) THEN
529
530 flag_rsl = &
531 ! zd>0 subject to FAI > beta**2*(1-PAI); also beta<0.5;
532 ! hence the lower limit of FAI below
533 ! FAI > 0.25*(1 - PAI) .AND. FAI < 0.45 .AND. & ! FAI
534 ! PAI > 0.1 .AND. PAI < 0.61 .AND. & ! PAI
535 pai > 0.1 .AND. pai < 0.68 .AND. & ! PAI
536 zh > 2 ! effective canopy height
537 ! LB Oct2021 - FAI and PAI can be larger than 0.45 and 0.61 respectively -> remove (1.-PAI)/FAI > .021 constraint
538 !(note: it seems wrong anyway - should be 0.87 not 0.021 based on G&O1991 numbers)
539 ELSE
540 ! default is to use MOST
541 flag_rsl = .false.
542 END IF
543
544 !
545 ! ! Step 2
546 ! ! determine vertical levels used in RSL
547 ! Define the height array with consideration of key heights
548 ! set number of heights within canopy
549 IF (zh <= 2) THEN
550 nz_can = 3
551 ELSE IF (zh <= 10) THEN
552 nz_can = 6
553 ELSE
554 nz_can = 15
555 END IF
556 ! fill up heights in canopy
557 dz_can = zh/nz_can
558 DO i = 1, nz_can
559 zarray(i) = dz_can*i
560 END DO
561
562 ! guaranttee 2 m is within the zarray
563 IF (dz_can > 2) zarray(1) = 1.999
564
565 ! fill up heights above canopy
566 nz_above = nz - nz_can
567 dz_above = (zmeas - zh)/nz_above
568 DO i = nz_can + 1, nz
569 zarray(i) = zh + (i - nz_can)*dz_above
570 END DO
571
572 ! ! determine index at the canyon top
573 ! DO z = 1, nz
574 ! dif(z) = ABS(zarray(z) - Zh)
575 ! END DO
576 ! nz_can = MINLOC(dif, DIM=1)
577 ! ! zarray(idx_can+2) = Zh+.1
578 ! ! zarray(idx_can+1) = Zh+.05
579 ! zarray(nz_can) = Zh
580 ! zarray(idx_can-1) = Zh-.1
581
582 ! determine index at measurement height
583 ! nz = nz
584 zarray(nz) = zmeas
585
586 IF (flag_rsl) THEN
587
588 ! use RSL approach to calculate correction factors
589 psihatm_z = 0.*zarray
590 psihath_z = 0.*zarray
591 ! Step 0: Calculate grid-cell dependent constants and Beta (crucial for H&F method)
592 CALL rsl_cal_prms( &
593 stabilitymethod, & !input
594 nz_above, zarray(nz_can + 1:nz), & !input
595 zh, l_mod, sfr_surf, fai, pai, & !input
596 psihatm_z(nz_can + 1:nz), psihath_z(nz_can + 1:nz), & !output
597 zh_rsl, l_mod_rsl, & ! output
598 lc, beta, zd_rsl, z0_rsl, elm, scc, fx)
599
600 ! Step 3: calculate the stability dependent H&F constants
601
602 ! CALL cal_ch(StabilityMethod, zh_RSL, zd_RSL, Lc, beta, L_MOD_RSL, Scc, fx, c2h, ch)
603 ! CALL cal_cm(StabilityMethod, zH_RSL, zd_RSL, Lc, beta, L_MOD_RSL, c2m, cm)
604
605 ! ! Step 4: determine psihat at levels above the canopy
606 ! DO z = nz - 1, idx_can, -1
607 ! phimz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
608 ! phimzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
609 ! phihz = stab_phi_heat(StabilityMethod, (zarray(z) - zd_RSL)/L_MOD_RSL)
610 ! phihzp = stab_phi_heat(StabilityMethod, (zarray(z + 1) - zd_RSL)/L_MOD_RSL)
611
612 ! 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
613 ! /(zarray(z + 1) - zd_RSL)
614 ! psihatm_z(z) = psihatm_z(z) + dz_above/2.*phimz*(cm*EXP(-1.*c2m*beta*(zarray(z) - zd_RSL)/elm)) &
615 ! /(zarray(z) - zd_RSL)
616 ! 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
617 ! /(zarray(z + 1) - zd_RSL)
618 ! psihath_z(z) = psihath_z(z) + dz_above/2.*phihz*(ch*EXP(-1.*c2h*beta*(zarray(z) - zd_RSL)/elm)) &
619 ! /(zarray(z) - zd_RSL)
620 ! END DO
621
622 ELSE
623
624 ! correct parameters if RSL approach doesn't apply for scenario of isolated flows
625 ! see Fig 1 of Grimmond and Oke (1999)
626 ! when isolated flow or skimming flow, implying RSL doesn't apply, force RSL correction to zero
627 psihatm_z = 0
628 psihath_z = 0
629
630 ! use L_MOD as in other parts of SUEWS
631 l_mod_rsl = l_mod
632
633 !correct RSL-based using SUEWS system-wide values
634 z0_rsl = z0m
635 zd_rsl = zdm
636 zh_rsl = zh
637
638 lc = -999
639 beta = -999
640 scc = -999
641 fx = -999
642 elm = -999
643
644 ! then MOST recovers from RSL correction
645 END IF
646
647 ! Step 6: Calculate mean variables above canopy
648 !
649 psimz0 = stab_psi_mom(stabilitymethod, z0_rsl/l_mod_rsl)
650 psimza = stab_psi_mom(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
651 psihza = stab_psi_heat(stabilitymethod, (zmeas - zd_rsl)/l_mod_rsl)
652
653 ustar_rsl = avu1*kappa/(log((zmeas - zd_rsl)/z0_rsl) - psimza + psimz0 + psihatm_z(nz))
654
655 ! TS 11 Feb 2021: limit UStar and TStar to reasonable ranges
656 ! 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
657 ustar_rsl = max(0.001, ustar_rsl)
658 ! under convective/unstable condition, min(UStar)==0.15 m s-1: (Schumann 1988, BLM, https://doi.org/10.1007/BF00123019)
659 IF ((zmeas - zd_rsl)/l_mod_rsl < -neut_limit) ustar_rsl = max(0.15, ustar_rsl)
660
661 ! TStar_RSL = -1.*(qh/(avcp*avdens))/UStar_RSL
662 ! qStar_RSL = -1.*(qe/lv_J_kg*avdens)/UStar_RSL
663 IF (flag_rsl) THEN
664 ustar_heat = max(0.15, ustar_rsl)
665 ELSE
666 ! use UStar_heat implied by RA_h using MOST
667 psihz0 = stab_psi_heat(stabilitymethod, z0v/l_mod_rsl)
668 ustar_heat = 1/(kappa*ra_h)*(log((zmeas - zd_rsl)/z0v) - psihza + psihz0)
669 END IF
670 tstar_rsl = -1.*(qh/(avcp*avdens))/ustar_heat
671 IF (qe == 0.0) THEN
672 qstar_rsl = 10.**(-10) ! avoid the situation where qe=0, qstar_RSL=0 and the code breaks LB 21 May 2021
673 ELSE
674 qstar_rsl = -1.*(qe/lv_j_kg*avdens)/ustar_heat
675 END IF
676 qa_gkg = rh2qa(avrh/100, press_hpa, temp_c)
677
678 DO z = nz_can, nz
679 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
680 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
681 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0 + psihatm_z(z))/kappa ! eqn. 3 in Theeuwes et al. (2019 BLM)
682 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)
683 dataoutlineqrsl(z) = dataoutlinetrsl(z)
684 END DO
685 !
686 ! Step 7: calculate in canopy variables
687 !
688 IF (flag_rsl) THEN
689 ! RSL approach: exponential profiles within canopy
690 IF (nz_can > 1) THEN
691 t_h = scc*tstar_rsl/(beta*fx)
692 q_h = scc*qstar_rsl/(beta*fx)
693 DO z = 1, nz_can - 1
694 dataoutlineursl(z) = dataoutlineursl(nz_can)*exp(beta*(zarray(z) - zh_rsl)/elm)
695 dataoutlinetrsl(z) = dataoutlinetrsl(nz_can) + (t_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - t_h)/tstar_rsl
696 dataoutlineqrsl(z) = dataoutlineqrsl(nz_can) + (q_h*exp(beta*fx*(zarray(z) - zh_rsl)/elm) - q_h)/qstar_rsl
697 END DO
698 END IF
699 ELSE
700 ! MOST approach:
701 DO z = 1, nz_can
702 ! when using MOST, all vertical levels should larger than zd_RSL
703 IF (zarray(z) <= zd_rsl) zarray(z) = 1.01*(zd_rsl + z0_rsl)
704 psimz = stab_psi_mom(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
705 psihz = stab_psi_heat(stabilitymethod, (zarray(z) - zd_rsl)/l_mod_rsl)
706 dataoutlineursl(z) = (log((zarray(z) - zd_rsl)/z0_rsl) - psimz + psimz0)/kappa
707 dataoutlinetrsl(z) = (log((zarray(z) - zd_rsl)/(zmeas - zd_rsl)) - psihz + psihza)/kappa
708 dataoutlineqrsl(z) = dataoutlinetrsl(z)
709 END DO
710 END IF
711
712 ! associate physical quantities to correction profilles
713 dataoutlineursl = dataoutlineursl*ustar_rsl
714 dataoutlinetrsl = dataoutlinetrsl*tstar_rsl + temp_c
715 dataoutlineqrsl = (dataoutlineqrsl*qstar_rsl + qa_gkg/1000.)*1000.
716
717 ! construct output line for output file
718 dataoutlinersl = [zarray, dataoutlineursl, dataoutlinetrsl, dataoutlineqrsl, &
719 !information for debugging
720 ! L_stab, L_unstab,
721 l_mod_rsl, &
722 zh_rsl, &
723 ! Lc_stab, Lc_unstab,
724 lc, &
725 beta, zd_rsl, z0_rsl, elm, scc, fx, &
726 ustar_rsl, ustar_heat, tstar_rsl, fai, pai, merge(1.d0, 0.d0, flag_rsl) &
727 ]
728
729 !
730 ! Step 8
731 ! retrieve the diagnostics at key heights
732 !
733 IF (flag_rsl) THEN
734 ! RSL approach: diagnostics within canopy, heights are above ground level
735 t2_c = interp_z(2d0, zarray, dataoutlinetrsl)
736 q2_gkg = interp_z(2d0, zarray, dataoutlineqrsl)
737 u10_ms = interp_z(10d0, zarray, dataoutlineursl)
738 ELSE
739 ! MOST approach: diagnostics at heights above zdm+z0m to avoid insane values
740 t2_c = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlinetrsl)
741 q2_gkg = interp_z(2d0 + zd_rsl + z0_rsl, zarray, dataoutlineqrsl)
742 u10_ms = interp_z(10d0 + zd_rsl + z0_rsl, zarray, dataoutlineursl)
743 END IF
744 ! get relative humidity:
745 rh2 = qa2rh(q2_gkg, press_hpa, t2_c)
746

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_dts().

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
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0