Skip to content
Snippets Groups Projects
pbldia_new_sfx.f90 8.19 KiB
Newer Older
module pbldia_new_sfx

use sfx_phys_const

implicit none 
real,parameter,private::undef=-999.0,R=287.,appa=0.287

contains

! in this subroutine interpolation is performed using integration of universal functions from z1 to z2, 
! where z1 is measurement or lowest model level height and z2 is required height
! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx)
!C*    DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND),
!C*    OF POTENTIAL TEMPERATURE DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE),
!C*    OF SPECIFIC HUMIDITY DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE)
!C*    INPUT DATA USED:
!C*    ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE)
!C*    AR2(1)    - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH)
!C*                HEIGHT OF CONSTANT FLUX LAYER TOP
!C*    AR2(5)    - ROUGHNESS LENGTH FOR MOMENTUM
!C*    AR2(6)    - ROUGHNESS LENGTH FOR TEMPERATURE AND HUMIDITY
!C*    AR2(8)    - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM
!C*    AR2(9)    - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY
!C*    ARRAY ARDIN
!C*    ARDIN(1)  - MODULE OF WIND VELOCITY AT THE TOP OF CONSTANT FLUX LAYER
!C*    ARDIN(2)  - POTENTIAL TEMPERATURE DEFICIT IN CONSTANT FLUX LAYER
!C*    ARDIN(3)  - SPECIFIC HUMIDITY DEFICIT IN CONSTANT FLUX LAYER
!C*    ARDIN(4)  - DIMENSIONAL HEIGHT OF CONSTANT FLUX LAYER TOP
!C*    ARDIN(5)  - = HWIND
!C*    ARDIN(6)  - = HTEMP
!C*    OUTPUT DATA (ARRAY ARDOUT):
!C*    ARDOUT(1)  - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT
!C*    ARDOUT(2)  - POTENTIAL TEMPERATURE at REQUIRED HEIGHT
!C*    ARDOUT(3)  - SPECIFIC HUMIDITY at REQUIRED HEIGHT
!C*    UFWIND     - UNIVERSAL FUNCTION FOR WIND VELOCITY
!C*    UFTEMP     - UNIVERSAL FUNCTION FOR SCALARS
!     AR2(1) - NON-DIMENSIONAL CFL HIGHT                            =
!     AR2(2) - RICHARDSON NUMBER                                    =
!     AR2(3) - REYNODS NUMBER                                       =
!     AR2(4) - LN(ZU/ZT)                                            =
!     AR2(5) - DYNAMICAL ROUGHNESS ZU (M)                           =
!     AR2(6) - THERMAL   ROUGHNESS ZT (M)                           =
!     AR2(7) - CRITICAL RICHARDSON NUMBER                           =
!     AR2(8) - TRANSFER COEFFICIENT FOR MOMENTUM                    =
!     AR2(9) - TRANSFER COEFFICIENT FR HEAT                         =
!     AR2(10)- COEFFICIENT OF TURBULENCE (KM) AT CFL HIGHT (M**2/S) =
!     AR2(11)- ALFT=KT/KM ( KT-COEFFICIENT OF TURBULENCE FOR HEAT)  =


subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
use sfx_sheba, only: &
                get_psi_sheba => get_psi

      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
      real,intent(out)::ARDOUT(3) !ARDOUT(4)
!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
      real,parameter::kappa = 0.4 
      real psi_m,psi_h,psi_m_hs,psi_h_hs
      real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta

      HWIND = ARDIN(5)
      HTEMP = ARDIN(6)
      HIN = ardin(4)
      zeta = AR2(1)

      USTAR = AR2(8) * ARDIN(1)
      TSTAR = AR2(9) * ARDIN(2)
      QSTAR = AR2(9) * ARDIN(3)


     if(zeta.gt.zetalim) zeta=zetalim
!     if(L.lt.Llim.and.L.gt.0) L=Llim
      L = HIN/zeta  ! MO length

     call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
     call get_psi_sheba(psi_m,psi_h,HWIND/L)
     UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
     WIND = (USTAR/kappa) * UFWIND
     call get_psi_sheba(psi_m,psi_h,HTEMP/L)
      UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)

      DTETA = (TSTAR/kappa) * UFTEMP
      DQ = (QSTAR/kappa) * UFTEMP


      ARDOUT(1) = ARDIN(1)+WIND
      ARDOUT(2) = DTETA+ardin(2)
      ARDOUT(3) = DQ+ardin(3)


      RETURN
  END subroutine pbldia_new_sheba

subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
use sfx_most, only: &
                get_psi_most => get_psi

      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
      real,intent(out)::ARDOUT(3) !ARDOUT(4)
!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
      real,parameter::kappa = 0.4 
      real psi_m,psi_h,psi_m_hs,psi_h_hs
      real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta

      HWIND = ARDIN(5)
      HTEMP = ARDIN(6)
      HIN = ardin(4)
      zeta = AR2(1)

      USTAR = AR2(8) * ARDIN(1)
      TSTAR = AR2(9) * ARDIN(2)
      QSTAR = AR2(9) * ARDIN(3)

     if(zeta.gt.zetalim) zeta=zetalim
!     if(L.lt.Llim.and.L.gt.0) L=Llim
      L = HIN/zeta  ! MO length

     call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
     call get_psi_most(psi_m,psi_h,HWIND/L)
     UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
     WIND = (USTAR/kappa) * UFWIND
     call get_psi_most(psi_m,psi_h,HTEMP/L)
      UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
!  if(L>0)then
!   !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
!      UFTEMP = ALOG(HTEMP/HIN) 
!  endif

      DTETA = (TSTAR/kappa) * UFTEMP
      DQ = (QSTAR/kappa) * UFTEMP

      ARDOUT(1) = ARDIN(1)+WIND
      ARDOUT(2) = DTETA+ardin(2)
      ARDOUT(3) = DQ+ardin(3)


      RETURN
  END subroutine pbldia_new_most

subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
!use sfx_esm, only: &
!                get_psi_esm => get_psi
      real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10)
      real,intent(out)::ARDOUT(3) !ARDOUT(4)
!      real,parameter::Llim = 50.0 !minimum value of L for stable SL (25)
      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
      real,parameter::kappa = 0.4 
      real psi_m,psi_h,psi_m_hs,psi_h_hs
      real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta

      HWIND = ARDIN(5)
      HTEMP = ARDIN(6)
      HIN = ardin(4)
      zeta = AR2(1)

      USTAR = AR2(8) * ARDIN(1)
      TSTAR = AR2(9) * ARDIN(2)
      QSTAR = AR2(9) * ARDIN(3)

     if(zeta.gt.zetalim) zeta=zetalim
!     if(L.lt.Llim.and.L.gt.0) L=Llim
      L = HIN/zeta  ! MO length


     call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
     UFWIND = psi_m
     WIND = (USTAR/kappa) * UFWIND
     call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
      UFTEMP = psi_h
!  if(L>0)then
!   !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T
!      UFTEMP = ALOG(HTEMP/HIN) 
!  endif

      DTETA = (TSTAR/kappa) * UFTEMP
      DQ = (QSTAR/kappa) * UFTEMP

      ARDOUT(1) = ARDIN(1)+WIND
      ARDOUT(2) = DTETA+ardin(2)
      ARDOUT(3) = DQ+ardin(3)


      RETURN
  END subroutine pbldia_new_esm


subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
use sfx_esm_param
    !< @brief universal functions (momentum) & (heat): neutral case
    ! ----------------------------------------------------------------------------
    real, intent(out) :: psi_m, psi_h   !< universal functions
    real, intent(in) :: z1,z2,L            !

  real an5,D1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2

      an5 = (Pr_t_inf_inv/Pr_t_0_inv)**4
      D1=(2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*&
&        (alpha_h-alpha_m)))/(2.0E0*alpha_h**2)


      zeta1 = z1 / L
      zeta2 = z2 / L

      IF(zeta2 .GE. 0.) THEN

      psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
      psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)

      ELSE

      GMZ2 = SQRT(SQRT(1.E0 - alpha_m * zeta2))
      GMZ1 = SQRT(SQRT(1.E0 - alpha_m * zeta1))
      GMST = SQRT(SQRT(1.E0 - alpha_m * D1))
      GTZ2 = SQRT(1.E0 - alpha_h * zeta2)
      GTZ1 = SQRT(1.E0 - alpha_h * zeta1)
      GTST = SQRT(1.E0 - alpha_h * D1)

      IF(zeta2 .GE. D1) THEN
      psi_m = FIM(GMZ2) - FIM(GMZ1)
      psi_h = FIT(GTZ2) - FIT(GTZ1)
      ELSE
      psi_m = (3.E0/GMST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIM(GMST) - FIM(GMZ1)
      psi_h = (3.E0/GTST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIT(GTST) - FIT(GTZ1)
      ENDIF
      ENDIF

end subroutine

  REAL FUNCTION FIM(XX)
      IMPLICIT NONE
      REAL X, XX
      X = AMAX1(XX,1.000001E0)
      FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X)
      RETURN
  END function

  REAL FUNCTION FIT(XX)
      IMPLICIT NONE
      REAL X, XX
      X = AMAX1(XX,1.000001E0)
      FIT = ALOG((X-1.E0)/(X+1.E0))
      RETURN
  END function



end module pbldia_new_sfx