From 6c50cbf4bd1253dc852aabe9246f71e83c574dd0 Mon Sep 17 00:00:00 2001 From: sliuclimate Date: Wed, 24 Jan 2024 16:01:01 +0800 Subject: [PATCH 1/3] format MOD_TurbulenceLEddy.F90 tuned --- main/MOD_TurbulenceLEddy.F90 | 516 +++++++++++++++++------------------ 1 file changed, 258 insertions(+), 258 deletions(-) diff --git a/main/MOD_TurbulenceLEddy.F90 b/main/MOD_TurbulenceLEddy.F90 index 7db4fefe..9b50d486 100644 --- a/main/MOD_TurbulenceLEddy.F90 +++ b/main/MOD_TurbulenceLEddy.F90 @@ -1,32 +1,32 @@ MODULE MOD_TurbulenceLEddy !----------------------------------------------------------------------- - use MOD_Precision - IMPLICIT NONE - SAVE + USE MOD_Precision + IMPLICIT NONE + SAVE ! PUBLIC MEMBER FUNCTIONS: - public :: moninobuk_leddy - public :: moninobukm_leddy + PUBLIC :: moninobuk_leddy + PUBLIC :: moninobukm_leddy ! PRIVATE MEMBER FUNCTIONS: - private :: psi + PRIVATE :: psi !----------------------------------------------------------------------- - CONTAINS +CONTAINS !----------------------------------------------------------------------- - subroutine moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & - ustar,fh2m,fq2m,fm10m,fm,fh,fq) + SUBROUTINE moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & + ustar,fh2m,fq2m,fm10m,fm,fh,fq) ! ====================================================================== ! ! Implement the LZD2022 scheme (Liu et al., 2022), which accounts for large -! eddy effects by inlcuding the boundary layer height in the phim function, +! eddy effects by inlcuding the boundary layer height in the phim FUNCTION, ! to compute friction velocity, relation for potential temperature and ! humidity profiles of surface boundary layer. ! @@ -41,167 +41,167 @@ subroutine moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & ! ! ====================================================================== - use MOD_Precision - use MOD_Const_Physical, only : vonkar - implicit none + USE MOD_Precision + USE MOD_Const_Physical, only : vonkar + IMPLICIT NONE ! ---------------------- dummy argument -------------------------------- - real(r8), INTENT(in) :: hu ! observational height of wind [m] - real(r8), INTENT(in) :: ht ! observational height of temperature [m] - real(r8), INTENT(in) :: hq ! observational height of humidity [m] - real(r8), INTENT(in) :: displa ! displacement height [m] - real(r8), INTENT(in) :: z0m ! roughness length, momentum [m] - real(r8), INTENT(in) :: z0h ! roughness length, sensible heat [m] - real(r8), INTENT(in) :: z0q ! roughness length, latent heat [m] - real(r8), INTENT(in) :: obu ! monin-obukhov length (m) - real(r8), INTENT(in) :: um ! wind speed including the stablity effect [m/s] - real(r8), INTENT(in) :: hpbl ! atmospheric boundary layer height [m] - - real(r8), INTENT(out) :: ustar ! friction velocity [m/s] - real(r8), INTENT(out) :: fh2m ! relation for temperature at 2m - real(r8), INTENT(out) :: fq2m ! relation for specific humidity at 2m - real(r8), INTENT(out) :: fm10m ! integral of profile function for momentum at 10m - real(r8), INTENT(out) :: fm ! integral of profile function for momentum - real(r8), INTENT(out) :: fh ! integral of profile function for heat - real(r8), INTENT(out) :: fq ! integral of profile function for moisture - + real(r8), intent(in) :: hu ! observational height of wind [m] + real(r8), intent(in) :: ht ! observational height of temperature [m] + real(r8), intent(in) :: hq ! observational height of humidity [m] + real(r8), intent(in) :: displa ! displacement height [m] + real(r8), intent(in) :: z0m ! roughness length, momentum [m] + real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] + real(r8), intent(in) :: z0q ! roughness length, latent heat [m] + real(r8), intent(in) :: obu ! monin-obukhov length (m) + real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] + real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] + + real(r8), intent(out) :: ustar ! friction velocity [m/s] + real(r8), intent(out) :: fh2m ! relation for temperature at 2m + real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m + real(r8), intent(out) :: fm10m ! integral of profile FUNCTION for momentum at 10m + real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum + real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat + real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture + !------------------------ local variables ------------------------------ - real(r8) zldis ! reference height "minus" zero displacement heght [m] - real(r8) zetam, & - zetam2 ! transition point of flux-gradient relation (wind profile) - real(r8) zetat ! transition point of flux-gradient relation (temp. profile) - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme - real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 - real(r8) Bm2 ! max(Bm, 0.2722) + real(r8) zldis ! reference height "minus" zero displacement heght [m] + real(r8) zetam, & + zetam2 ! transition point of flux-gradient relation (wind profile) + real(r8) zetat ! transition point of flux-gradient relation (temp. profile) + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme + real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 + real(r8) Bm2 ! max(Bm, 0.2722) -! real(r8), external :: psi ! stability function for unstable case +! real(r8), external :: psi ! stability FUNCTION for unstable CASE !----------------------------------------------------------------------- ! adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions. ! wind profile - zldis=hu-displa - zeta=zldis/obu + zldis=hu-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.05 ! zetazi = max(5.*hu, hpbl)/obu - if(zetazi >= 0.) then !stable - zetazi = min(200.,max(zetazi,1.e-5)) - else !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - endif + IF(zetazi >= 0.) THEN !stable + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF Bm = 0.0047 * (-zetazi) + 0.1854 zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) Bm2 = max(Bm, 0.2722) zetam2 = min(zetam, -0.13) - if(zeta < zetam2)then ! zeta < zetam2 - fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) - ustar = vonkar*um / fm + ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - else if(zeta < 0.)then ! zetam2 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) ustar = vonkar*um / fm - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu ustar = vonkar*um / fm - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) ustar = vonkar*um / fm - endif + ENDIF - ! for 10 meter wind-velocity - zldis=10.+z0m - zeta=zldis/obu +! for 10 meter wind-velocity + zldis=10.+z0m + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! - if(zeta < zetam2)then ! zeta < zetam2 + IF(zeta < zetam2)THEN ! zeta < zetam2 fm10m = log(zetam2*obu/z0m) - psi(1,zetam2) & + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - else if(zeta < 0.)then ! zetam2 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 fm10m = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fm10m = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fm10m = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF - ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 +! for 2 meter screen temperature + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - endif - - ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - elseif (zeta < 0.) then ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - else if (zeta <= 1.) then ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - else ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - endif + ENDIF + +! for 2 meter screen humidity + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF - end subroutine moninobuk_leddy + END SUBROUTINE moninobuk_leddy - subroutine moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpbl, & - ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) + SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpbl, & + ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) ! ====================================================================== ! @@ -221,243 +221,243 @@ subroutine moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpb ! fm, fq and phih for roughness sublayer u/k profile calculation ! Shaofeng Liu, 05/2023: implement the LZD2022 scheme (Liu et al., 2022), which ! accounts for large eddy effects by including the -! boundary leyer height in the phim function. +! boundary leyer height in the phim FUNCTION. ! ====================================================================== - use MOD_Precision - use MOD_Const_Physical, only : vonkar - implicit none + USE MOD_Precision + USE MOD_Const_Physical, only : vonkar + IMPLICIT NONE ! ---------------------- dummy argument -------------------------------- - real(r8), INTENT(in) :: hu ! observational height of wind [m] - real(r8), INTENT(in) :: ht ! observational height of temperature [m] - real(r8), INTENT(in) :: hq ! observational height of humidity [m] - real(r8), INTENT(in) :: displa ! displacement height [m] - real(r8), INTENT(in) :: displat ! displacement height of the top layer [m] - real(r8), INTENT(in) :: z0m ! roughness length, momentum [m] - real(r8), INTENT(in) :: z0h ! roughness length, sensible heat [m] - real(r8), INTENT(in) :: z0q ! roughness length, latent heat [m] - real(r8), INTENT(in) :: z0mt ! roughness length of the top layer, latent heat [m] - real(r8), INTENT(in) :: htop ! canopy top height of the top layer [m] - real(r8), INTENT(in) :: obu ! monin-obukhov length (m) - real(r8), INTENT(in) :: um ! wind speed including the stablity effect [m/s] - real(r8), INTENT(in) :: hpbl ! atmospheric boundary layer height [m] - - real(r8), INTENT(out) :: ustar ! friction velocity [m/s] - real(r8), INTENT(out) :: fh2m ! relation for temperature at 2m - real(r8), INTENT(out) :: fq2m ! relation for specific humidity at 2m - real(r8), INTENT(out) :: fmtop ! integral of profile function for momentum at 10m - real(r8), INTENT(out) :: fm ! integral of profile function for momentum - real(r8), INTENT(out) :: fh ! integral of profile function for heat - real(r8), INTENT(out) :: fq ! integral of profile function for moisture - real(r8), INTENT(out) :: fht ! integral of profile function for heat at the top layer - real(r8), INTENT(out) :: fqt ! integral of profile function for moisture at the top layer - real(r8), INTENT(out) :: phih ! phi(h), similarity function for sensible heat + real(r8), intent(in) :: hu ! observational height of wind [m] + real(r8), intent(in) :: ht ! observational height of temperature [m] + real(r8), intent(in) :: hq ! observational height of humidity [m] + real(r8), intent(in) :: displa ! displacement height [m] + real(r8), intent(in) :: displat ! displacement height of the top layer [m] + real(r8), intent(in) :: z0m ! roughness length, momentum [m] + real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] + real(r8), intent(in) :: z0q ! roughness length, latent heat [m] + real(r8), intent(in) :: z0mt ! roughness length of the top layer, latent heat [m] + real(r8), intent(in) :: htop ! canopy top height of the top layer [m] + real(r8), intent(in) :: obu ! monin-obukhov length (m) + real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] + real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] + + real(r8), intent(out) :: ustar ! friction velocity [m/s] + real(r8), intent(out) :: fh2m ! relation for temperature at 2m + real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m + real(r8), intent(out) :: fmtop ! integral of profile FUNCTION for momentum at 10m + real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum + real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat + real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture + real(r8), intent(out) :: fht ! integral of profile FUNCTION for heat at the top layer + real(r8), intent(out) :: fqt ! integral of profile FUNCTION for moisture at the top layer + real(r8), intent(out) :: phih ! phi(h), similarity FUNCTION for sensible heat !------------------------ local variables ------------------------------ - real(r8) zldis ! reference height "minus" zero displacement heght [m] - real(r8) zetam, & - zetam2 ! transition point of flux-gradient relation (wind profile) - real(r8) zetat ! transition point of flux-gradient relation (temp. profile) - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme - real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 - real(r8) Bm2 ! max(Bm, 0.2722) + real(r8) zldis ! reference height "minus" zero displacement heght [m] + real(r8) zetam, & + zetam2 ! transition point of flux-gradient relation (wind profile) + real(r8) zetat ! transition point of flux-gradient relation (temp. profile) + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme + real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 + real(r8) Bm2 ! max(Bm, 0.2722) -! real(r8), external :: psi ! stability function for unstable case +! real(r8), external :: psi ! stability FUNCTION for unstable CASE !----------------------------------------------------------------------- ! adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions. ! wind profile - zldis=hu-displa - zeta=zldis/obu + zldis=hu-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.05 ! ! zetazi = hpbl/obu zetazi = max(5.*hu, hpbl)/obu - if(zetazi >= 0.) then !stable - zetazi = min(200.,max(zetazi,1.e-5)) - else !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - endif + IF(zetazi >= 0.) THEN !stable + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF Bm = 0.0047 * (-zetazi) + 0.1854 zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) Bm2 = max(Bm, 0.2722) zetam2 = min(zetam, -0.13) - if(zeta < zetam2)then ! zeta < zetam2 + IF(zeta < zetam2)THEN ! zeta < zetam2 fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - else if(zeta < 0.)then ! zetam2 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) ustar = vonkar*um / fm - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu ustar = vonkar*um / fm - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) ustar = vonkar*um / fm - endif + ENDIF - ! for canopy top wind-velocity - !NOTE: changed for canopy top wind-velocity (no wake assumed) - zldis=htop-displa - zeta=zldis/obu +! for canopy top wind-velocity +!NOTE: changed for canopy top wind-velocity (no wake assumed) + zldis=htop-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! ! zetam=1.574 - if(zeta < zetam2)then ! zeta < zetam2 + IF(zeta < zetam2)THEN ! zeta < zetam2 fmtop = log(zetam2*obu/z0m) - psi(1,zetam2) & + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - else if(zeta < 0.)then ! zetam2 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 fmtop = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fmtop = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fmtop = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF - ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 +! for 2 meter screen temperature + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - endif + ENDIF - ! for top layer temperature - zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 +! for top layer temperature + zldis=displat+z0mt-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fht = log(-zetat*obu/z0h)-psi(2,-zetat) & + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fht = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fht = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fht = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - endif - - ! for canopy top phi(h) - ! CESM TECH NOTE EQ. (5.31) - zldis=htop-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 - phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) - else if(zeta < 0.)then ! -1 <= zeta < 0 - phih = (1. - 16.*zeta)**(-0.5) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 - phih = 1. + 5.*zeta - else ! 1 < zeta, phi=5+zeta - phih = 5. + zeta - endif + ENDIF + +! for canopy top phi(h) +! CESM TECH NOTE eq. (5.31) + zldis=htop-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + phih = (1. - 16.*zeta)**(-0.5) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + phih = 1. + 5.*zeta + ELSE ! 1 < zeta, phi=5+zeta + phih = 5. + zeta + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - else if(zeta < 0.)then ! -1 <= zeta < 0 + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - else if(zeta <= 1.)then ! 0 <= ztea <= 1 + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - else ! 1 < zeta, phi=5+zeta + ELSE ! 1 < zeta, phi=5+zeta fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - endif - - ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - elseif (zeta < 0.) then ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - else if (zeta <= 1.) then ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - else ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - endif - - ! for top layer humidity - zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - if(zeta < -zetat)then ! zeta < -1 - fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - elseif (zeta < 0.) then ! -1 <= zeta < 0 - fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - else if (zeta <= 1.) then ! 0 <= zeta <= 1 - fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu - else ! 1 < zeta, phi=5+zeta - fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - endif - - end subroutine moninobukm_leddy + ENDIF + +! for 2 meter screen humidity + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF + +! for top layer humidity + zldis=displat+z0mt-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF + + END SUBROUTINE moninobukm_leddy - real(r8) function psi(k,zeta) + real(r8) FUNCTION psi(k,zeta) !======================================================================= -! stability function for unstable case (rib < 0) +! stability FUNCTION for unstable CASE (rib < 0) - use MOD_Precision - implicit none + USE MOD_Precision + IMPLICIT NONE - integer k - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) chik ! + integer k + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) chik ! chik = (1.-16.*zeta)**0.25 - if(k == 1)then - psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) - else - psi = 2.*log((1.+chik*chik)*0.5) - endif + IF(k == 1)THEN + psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) + ELSE + psi = 2.*log((1.+chik*chik)*0.5) + ENDIF - end function psi + END FUNCTION psi END MODULE MOD_TurbulenceLEddy From 5b7f4dcad99c54a0abd1a43332ad371e72e66d47 Mon Sep 17 00:00:00 2001 From: sliuclimate Date: Thu, 25 Jan 2024 13:29:21 +0800 Subject: [PATCH 2/3] MOD_TurbulenceLEddy.F90 format tuned --- main/MOD_TurbulenceLEddy.F90 | 412 +++++++++++++++++------------------ 1 file changed, 206 insertions(+), 206 deletions(-) diff --git a/main/MOD_TurbulenceLEddy.F90 b/main/MOD_TurbulenceLEddy.F90 index 9b50d486..7cdb5e9b 100644 --- a/main/MOD_TurbulenceLEddy.F90 +++ b/main/MOD_TurbulenceLEddy.F90 @@ -89,119 +89,119 @@ SUBROUTINE moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & ! zetazi = max(5.*hu, hpbl)/obu IF(zetazi >= 0.) THEN !stable - zetazi = min(200.,max(zetazi,1.e-5)) - ELSE !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - ENDIF + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF Bm = 0.0047 * (-zetazi) + 0.1854 zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) Bm2 = max(Bm, 0.2722) zetam2 = min(zetam, -0.13) - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ustar = vonkar*um / fm - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ustar = vonkar*um / fm + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ustar = vonkar*um / fm + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ustar = vonkar*um / fm ELSE ! 1 < zeta, phi=5+zeta - fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ustar = vonkar*um / fm - ENDIF + fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ustar = vonkar*um / fm + ENDIF ! for 10 meter wind-velocity - zldis=10.+z0m - zeta=zldis/obu + zldis=10.+z0m + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm10m = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm10m = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm10m = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm10m = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ELSE ! 1 < zeta, phi=5+zeta - fm10m = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm10m = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm10m = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ELSE ! 1 < zeta, phi=5+zeta + fm10m = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq = log(-zetat*obu/z0q) - psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF END SUBROUTINE moninobuk_leddy SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpbl, & - ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) + ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) ! ====================================================================== ! @@ -278,161 +278,161 @@ SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpb ! ! zetazi = hpbl/obu zetazi = max(5.*hu, hpbl)/obu - IF(zetazi >= 0.) THEN !stable - zetazi = min(200.,max(zetazi,1.e-5)) - ELSE !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - ENDIF + IF(zetazi >= 0.) THEN !stable + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF Bm = 0.0047 * (-zetazi) + 0.1854 zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) Bm2 = max(Bm, 0.2722) zetam2 = min(zetam, -0.13) - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) - ustar = vonkar*um / fm + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ustar = vonkar*um / fm - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ustar = vonkar*um / fm - ELSE ! 1 < zeta, phi=5+zeta - fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ustar = vonkar*um / fm - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ustar = vonkar*um / fm + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ustar = vonkar*um / fm + ELSE ! 1 < zeta, phi=5+zeta + fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ustar = vonkar*um / fm + ENDIF ! for canopy top wind-velocity !NOTE: changed for canopy top wind-velocity (no wake assumed) - zldis=htop-displa - zeta=zldis/obu + zldis=htop-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! ! zetam=1.574 IF(zeta < zetam2)THEN ! zeta < zetam2 - fmtop = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + fmtop = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fmtop = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fmtop = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ELSE ! 1 < zeta, phi=5+zeta - fmtop = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fmtop = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fmtop = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ELSE ! 1 < zeta, phi=5+zeta + fmtop = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for top layer temperature zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fht = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fht = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fht = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fht = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fht = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fht = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fht = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fht = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for canopy top phi(h) ! CESM TECH NOTE eq. (5.31) - zldis=htop-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - phih = (1. - 16.*zeta)**(-0.5) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - phih = 1. + 5.*zeta - ELSE ! 1 < zeta, phi=5+zeta - phih = 5. + zeta - ENDIF + zldis=htop-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + phih = (1. - 16.*zeta)**(-0.5) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + phih = 1. + 5.*zeta + ELSE ! 1 < zeta, phi=5+zeta + phih = 5. + zeta + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq = log(-zetat*obu/z0q) - psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF ! for top layer humidity - zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF + zldis=displat+z0mt-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF END SUBROUTINE moninobukm_leddy @@ -450,12 +450,12 @@ real(r8) FUNCTION psi(k,zeta) real(r8) zeta ! dimensionless height used in Monin-Obukhov theory real(r8) chik ! - chik = (1.-16.*zeta)**0.25 - IF(k == 1)THEN - psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) - ELSE - psi = 2.*log((1.+chik*chik)*0.5) - ENDIF + chik = (1.-16.*zeta)**0.25 + IF(k == 1)THEN + psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) + ELSE + psi = 2.*log((1.+chik*chik)*0.5) + ENDIF END FUNCTION psi From a18ce7dbf90838c590bb011ee153dad374658237 Mon Sep 17 00:00:00 2001 From: sliuclimate Date: Thu, 25 Jan 2024 16:14:47 +0800 Subject: [PATCH 3/3] retuned format MOD_TurbulenceLEddy.F90 --- main/MOD_TurbulenceLEddy.F90 | 618 +++++++++++++++++------------------ 1 file changed, 309 insertions(+), 309 deletions(-) diff --git a/main/MOD_TurbulenceLEddy.F90 b/main/MOD_TurbulenceLEddy.F90 index 7cdb5e9b..d79f88b9 100644 --- a/main/MOD_TurbulenceLEddy.F90 +++ b/main/MOD_TurbulenceLEddy.F90 @@ -20,7 +20,7 @@ MODULE MOD_TurbulenceLEddy !----------------------------------------------------------------------- - SUBROUTINE moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & + SUBROUTINE moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & ustar,fh2m,fq2m,fm10m,fm,fh,fq) ! ====================================================================== @@ -41,167 +41,167 @@ SUBROUTINE moninobuk_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um, hpbl, & ! ! ====================================================================== - USE MOD_Precision - USE MOD_Const_Physical, only : vonkar - IMPLICIT NONE + USE MOD_Precision + USE MOD_Const_Physical, only : vonkar + IMPLICIT NONE ! ---------------------- dummy argument -------------------------------- - real(r8), intent(in) :: hu ! observational height of wind [m] - real(r8), intent(in) :: ht ! observational height of temperature [m] - real(r8), intent(in) :: hq ! observational height of humidity [m] - real(r8), intent(in) :: displa ! displacement height [m] - real(r8), intent(in) :: z0m ! roughness length, momentum [m] - real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] - real(r8), intent(in) :: z0q ! roughness length, latent heat [m] - real(r8), intent(in) :: obu ! monin-obukhov length (m) - real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] - real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] - - real(r8), intent(out) :: ustar ! friction velocity [m/s] - real(r8), intent(out) :: fh2m ! relation for temperature at 2m - real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m - real(r8), intent(out) :: fm10m ! integral of profile FUNCTION for momentum at 10m - real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum - real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat - real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture + real(r8), intent(in) :: hu ! observational height of wind [m] + real(r8), intent(in) :: ht ! observational height of temperature [m] + real(r8), intent(in) :: hq ! observational height of humidity [m] + real(r8), intent(in) :: displa ! displacement height [m] + real(r8), intent(in) :: z0m ! roughness length, momentum [m] + real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] + real(r8), intent(in) :: z0q ! roughness length, latent heat [m] + real(r8), intent(in) :: obu ! monin-obukhov length (m) + real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] + real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] + + real(r8), intent(out) :: ustar ! friction velocity [m/s] + real(r8), intent(out) :: fh2m ! relation for temperature at 2m + real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m + real(r8), intent(out) :: fm10m ! integral of profile FUNCTION for momentum at 10m + real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum + real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat + real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture !------------------------ local variables ------------------------------ - real(r8) zldis ! reference height "minus" zero displacement heght [m] - real(r8) zetam, & - zetam2 ! transition point of flux-gradient relation (wind profile) - real(r8) zetat ! transition point of flux-gradient relation (temp. profile) - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme - real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 - real(r8) Bm2 ! max(Bm, 0.2722) + real(r8) zldis ! reference height "minus" zero displacement heght [m] + real(r8) zetam, & + zetam2 ! transition point of flux-gradient relation (wind profile) + real(r8) zetat ! transition point of flux-gradient relation (temp. profile) + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme + real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 + real(r8) Bm2 ! max(Bm, 0.2722) ! real(r8), external :: psi ! stability FUNCTION for unstable CASE !----------------------------------------------------------------------- ! adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions. ! wind profile - zldis=hu-displa - zeta=zldis/obu + zldis=hu-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.05 ! - zetazi = max(5.*hu, hpbl)/obu - IF(zetazi >= 0.) THEN !stable - zetazi = min(200.,max(zetazi,1.e-5)) - ELSE !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - ENDIF - - Bm = 0.0047 * (-zetazi) + 0.1854 - zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) - Bm2 = max(Bm, 0.2722) - zetam2 = min(zetam, -0.13) - - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + zetazi = max(5.*hu, hpbl)/obu + IF(zetazi >= 0.) THEN !stable + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF + + Bm = 0.0047 * (-zetazi) + 0.1854 + zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) + Bm2 = max(Bm, 0.2722) + zetam2 = min(zetam, -0.13) + + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ustar = vonkar*um / fm - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ustar = vonkar*um / fm + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ustar = vonkar*um / fm + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ustar = vonkar*um / fm ELSE ! 1 < zeta, phi=5+zeta - fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ustar = vonkar*um / fm - ENDIF + fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ustar = vonkar*um / fm + ENDIF ! for 10 meter wind-velocity - zldis=10.+z0m - zeta=zldis/obu + zldis=10.+z0m + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm10m = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm10m = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm10m = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm10m = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ELSE ! 1 < zeta, phi=5+zeta - fm10m = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm10m = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm10m = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ELSE ! 1 < zeta, phi=5+zeta + fm10m = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq = log(-zetat*obu/z0q) - psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF - - END SUBROUTINE moninobuk_leddy - - - SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpbl, & - ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF + + END SUBROUTINE moninobuk_leddy + + + SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpbl, & + ustar,fh2m,fq2m,htop,fmtop,fm,fh,fq,fht,fqt,phih) ! ====================================================================== ! @@ -224,240 +224,240 @@ SUBROUTINE moninobukm_leddy(hu,ht,hq,displa,z0m,z0h,z0q,obu,um,displat,z0mt, hpb ! boundary leyer height in the phim FUNCTION. ! ====================================================================== - USE MOD_Precision - USE MOD_Const_Physical, only : vonkar - IMPLICIT NONE + USE MOD_Precision + USE MOD_Const_Physical, only : vonkar + IMPLICIT NONE ! ---------------------- dummy argument -------------------------------- - real(r8), intent(in) :: hu ! observational height of wind [m] - real(r8), intent(in) :: ht ! observational height of temperature [m] - real(r8), intent(in) :: hq ! observational height of humidity [m] - real(r8), intent(in) :: displa ! displacement height [m] - real(r8), intent(in) :: displat ! displacement height of the top layer [m] - real(r8), intent(in) :: z0m ! roughness length, momentum [m] - real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] - real(r8), intent(in) :: z0q ! roughness length, latent heat [m] - real(r8), intent(in) :: z0mt ! roughness length of the top layer, latent heat [m] - real(r8), intent(in) :: htop ! canopy top height of the top layer [m] - real(r8), intent(in) :: obu ! monin-obukhov length (m) - real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] - real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] - - real(r8), intent(out) :: ustar ! friction velocity [m/s] - real(r8), intent(out) :: fh2m ! relation for temperature at 2m - real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m - real(r8), intent(out) :: fmtop ! integral of profile FUNCTION for momentum at 10m - real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum - real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat - real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture - real(r8), intent(out) :: fht ! integral of profile FUNCTION for heat at the top layer - real(r8), intent(out) :: fqt ! integral of profile FUNCTION for moisture at the top layer - real(r8), intent(out) :: phih ! phi(h), similarity FUNCTION for sensible heat + real(r8), intent(in) :: hu ! observational height of wind [m] + real(r8), intent(in) :: ht ! observational height of temperature [m] + real(r8), intent(in) :: hq ! observational height of humidity [m] + real(r8), intent(in) :: displa ! displacement height [m] + real(r8), intent(in) :: displat ! displacement height of the top layer [m] + real(r8), intent(in) :: z0m ! roughness length, momentum [m] + real(r8), intent(in) :: z0h ! roughness length, sensible heat [m] + real(r8), intent(in) :: z0q ! roughness length, latent heat [m] + real(r8), intent(in) :: z0mt ! roughness length of the top layer, latent heat [m] + real(r8), intent(in) :: htop ! canopy top height of the top layer [m] + real(r8), intent(in) :: obu ! monin-obukhov length (m) + real(r8), intent(in) :: um ! wind speed including the stablity effect [m/s] + real(r8), intent(in) :: hpbl ! atmospheric boundary layer height [m] + + real(r8), intent(out) :: ustar ! friction velocity [m/s] + real(r8), intent(out) :: fh2m ! relation for temperature at 2m + real(r8), intent(out) :: fq2m ! relation for specific humidity at 2m + real(r8), intent(out) :: fmtop ! integral of profile FUNCTION for momentum at 10m + real(r8), intent(out) :: fm ! integral of profile FUNCTION for momentum + real(r8), intent(out) :: fh ! integral of profile FUNCTION for heat + real(r8), intent(out) :: fq ! integral of profile FUNCTION for moisture + real(r8), intent(out) :: fht ! integral of profile FUNCTION for heat at the top layer + real(r8), intent(out) :: fqt ! integral of profile FUNCTION for moisture at the top layer + real(r8), intent(out) :: phih ! phi(h), similarity FUNCTION for sensible heat !------------------------ local variables ------------------------------ - real(r8) zldis ! reference height "minus" zero displacement heght [m] - real(r8) zetam, & - zetam2 ! transition point of flux-gradient relation (wind profile) - real(r8) zetat ! transition point of flux-gradient relation (temp. profile) - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme - real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 - real(r8) Bm2 ! max(Bm, 0.2722) + real(r8) zldis ! reference height "minus" zero displacement heght [m] + real(r8) zetam, & + zetam2 ! transition point of flux-gradient relation (wind profile) + real(r8) zetat ! transition point of flux-gradient relation (temp. profile) + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) zetazi ! hpbl/obu, dimensionless height used in the LZD2022 scheme + real(r8) Bm ! Coefficient of the LZD2022 scheme: Bm = 0.0047*(-hpbl/L) + 0.1854 + real(r8) Bm2 ! max(Bm, 0.2722) ! real(r8), external :: psi ! stability FUNCTION for unstable CASE !----------------------------------------------------------------------- ! adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions. ! wind profile - zldis=hu-displa - zeta=zldis/obu + zldis=hu-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.05 ! ! zetazi = hpbl/obu - zetazi = max(5.*hu, hpbl)/obu - IF(zetazi >= 0.) THEN !stable - zetazi = min(200.,max(zetazi,1.e-5)) - ELSE !unstable - zetazi = max(-1.e4,min(zetazi,-1.e-5)) - ENDIF - - Bm = 0.0047 * (-zetazi) + 0.1854 - zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) - Bm2 = max(Bm, 0.2722) - zetam2 = min(zetam, -0.13) - - IF(zeta < zetam2)THEN ! zeta < zetam2 - fm = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) - ustar = vonkar*um / fm + zetazi = max(5.*hu, hpbl)/obu + IF(zetazi >= 0.) THEN !stable + zetazi = min(200.,max(zetazi,1.e-5)) + ELSE !unstable + zetazi = max(-1.e4,min(zetazi,-1.e-5)) + ENDIF + + Bm = 0.0047 * (-zetazi) + 0.1854 + zetam = 0.5*Bm**4 * ( -16. - sqrt(256. + 4./Bm**4) ) + Bm2 = max(Bm, 0.2722) + zetam2 = min(zetam, -0.13) + + IF(zeta < zetam2)THEN ! zeta < zetam2 + fm = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + ustar = vonkar*um / fm ! ! End: Shaofeng Liu, 2023.05.05 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ustar = vonkar*um / fm - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ustar = vonkar*um / fm - ELSE ! 1 < zeta, phi=5+zeta - fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ustar = vonkar*um / fm - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fm = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ustar = vonkar*um / fm + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fm = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ustar = vonkar*um / fm + ELSE ! 1 < zeta, phi=5+zeta + fm = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ustar = vonkar*um / fm + ENDIF ! for canopy top wind-velocity !NOTE: changed for canopy top wind-velocity (no wake assumed) - zldis=htop-displa - zeta=zldis/obu + zldis=htop-displa + zeta=zldis/obu ! ! Begin: Shaofeng Liu, 2023.05.18 ! ! zetam=1.574 - IF(zeta < zetam2)THEN ! zeta < zetam2 - fmtop = log(zetam2*obu/z0m) - psi(1,zetam2) & - + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) + IF(zeta < zetam2)THEN ! zeta < zetam2 + fmtop = log(zetam2*obu/z0m) - psi(1,zetam2) & + + psi(1,z0m/obu) - 2.*Bm2 * ( (-zeta)**(-0.5)-(-zetam2)**(-0.5) ) ! ! End: Shaofeng Liu, 2023.05.18 ! - ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 - fmtop = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fmtop = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu - ELSE ! 1 < zeta, phi=5+zeta - fmtop = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) - ENDIF + ELSE IF(zeta < 0.)THEN ! zetam2 <= zeta < 0 + fmtop = log(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fmtop = log(zldis/z0m) + 5.*zeta - 5.*z0m/obu + ELSE ! 1 < zeta, phi=5+zeta + fmtop = log(obu/z0m) + 5. - 5.*z0m/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! temperature profile - zldis=ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen temperature - zldis=2.+z0h ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fh2m = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fh2m = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fh2m = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fh2m = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for top layer temperature - zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fht = log(-zetat*obu/z0h)-psi(2,-zetat) & - + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fht = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fht = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu - ELSE ! 1 < zeta, phi=5+zeta - fht = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=displat+z0mt-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fht = log(-zetat*obu/z0h)-psi(2,-zetat) & + + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fht = log(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fht = log(zldis/z0h) + 5.*zeta - 5.*z0h/obu + ELSE ! 1 < zeta, phi=5+zeta + fht = log(obu/z0h) + 5. - 5.*z0h/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for canopy top phi(h) ! CESM TECH NOTE eq. (5.31) - zldis=htop-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - phih = (1. - 16.*zeta)**(-0.5) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - phih = 1. + 5.*zeta - ELSE ! 1 < zeta, phi=5+zeta - phih = 5. + zeta - ENDIF + zldis=htop-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + phih = 0.9*vonkar**(1.333)*(-zeta)**(-0.333) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + phih = (1. - 16.*zeta)**(-0.5) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + phih = 1. + 5.*zeta + ELSE ! 1 < zeta, phi=5+zeta + phih = 5. + zeta + ENDIF ! humidity profile - zldis=hq-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq = log(-zetat*obu/z0q) - psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 - fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) - ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 - fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) - ENDIF + zldis=hq-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq = log(-zetat*obu/z0q) - psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF(zeta < 0.)THEN ! -1 <= zeta < 0 + fq = log(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) + ELSE IF(zeta <= 1.)THEN ! 0 <= ztea <= 1 + fq = log(zldis/z0q) + 5.*zeta - 5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq = log(obu/z0q) + 5. - 5.*z0q/obu + (5.*log(zeta)+zeta-1.) + ENDIF ! for 2 meter screen humidity - zldis=2.+z0h - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF + zldis=2.+z0h + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fq2m = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fq2m = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fq2m = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fq2m = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF ! for top layer humidity - zldis=displat+z0mt-displa ! ht-displa - zeta=zldis/obu - zetat=0.465 - IF(zeta < -zetat)THEN ! zeta < -1 - fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & - + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) - ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 - fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) - ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 - fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu - ELSE ! 1 < zeta, phi=5+zeta - fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) - ENDIF - - END SUBROUTINE moninobukm_leddy + zldis=displat+z0mt-displa ! ht-displa + zeta=zldis/obu + zetat=0.465 + IF(zeta < -zetat)THEN ! zeta < -1 + fqt = log(-zetat*obu/z0q)-psi(2,-zetat) & + + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) + ELSE IF (zeta < 0.) THEN ! -1 <= zeta < 0 + fqt = log(zldis/z0q)-psi(2,zeta)+psi(2,z0q/obu) + ELSE IF (zeta <= 1.) THEN ! 0 <= zeta <= 1 + fqt = log(zldis/z0q)+5.*zeta-5.*z0q/obu + ELSE ! 1 < zeta, phi=5+zeta + fqt = log(obu/z0q)+5.-5.*z0q/obu+(5.*log(zeta)+zeta-1.) + ENDIF + + END SUBROUTINE moninobukm_leddy - real(r8) FUNCTION psi(k,zeta) + real(r8) FUNCTION psi(k,zeta) !======================================================================= ! stability FUNCTION for unstable CASE (rib < 0) - USE MOD_Precision - IMPLICIT NONE - - integer k - real(r8) zeta ! dimensionless height used in Monin-Obukhov theory - real(r8) chik ! - - chik = (1.-16.*zeta)**0.25 - IF(k == 1)THEN - psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) - ELSE - psi = 2.*log((1.+chik*chik)*0.5) - ENDIF - - END FUNCTION psi + USE MOD_Precision + IMPLICIT NONE + + integer k + real(r8) zeta ! dimensionless height used in Monin-Obukhov theory + real(r8) chik ! + + chik = (1.-16.*zeta)**0.25 + IF(k == 1)THEN + psi = 2.*log((1.+chik)*0.5)+log((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.) + ELSE + psi = 2.*log((1.+chik*chik)*0.5) + ENDIF + + END FUNCTION psi END MODULE MOD_TurbulenceLEddy