Skip to content
Snippets Groups Projects
Commit ab80b72e authored by Anna Shestakova's avatar Anna Shestakova
Browse files

Delete pbldia_new_sfx.f90

parent d1fc905f
Branches
No related tags found
No related merge requests found
module pbldia_new_sfx
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
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment