module pbldia_new_sfx
implicit none
private
public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
pbldia_new_most,pbldia_new_esm,pbldia_new_log,pbldia_new_sheba_coare
real,parameter,private::kappa=0.4,R=287.,appa=0.287
contains
! in this module 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(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) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
!C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
!C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
!C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
subroutine pbldia_new_log(AR2,ARDIN,ARDOUT)
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real hwind,htemp,ustar,tstar,qstar
real dteta,dq,wind,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)
WIND = (USTAR/kappa)*ALOG(hwind/HIN)
DTETA = (TSTAR/kappa)*ALOG(htemp/HIN)
DQ = (QSTAR/kappa)*ALOG(htemp/HIN)
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_log
subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
use sfx_sheba, only: get_psi_sheba => get_psi
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for
! stable SL for wind
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
real Rib,ksi,wg,f,L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if (zeta.ne.0) then
L = HIN/zeta
if(zeta.gt.zetalim) then
L_lim = HIN/zetalim
else
L_lim = L
end if
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L_lim)
call get_psi_sheba(psi_m,psi_h,HWIND/L_lim)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_sheba(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
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_sheba_coare(AR2,ARDIN,ARDOUT)
use sfx_sheba_coare, only: &
get_psi_coare => get_psi_a, &
get_psi_stable_sheba => get_psi_stable
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
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,L_lim
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.ne.0)then
L = HIN/zeta
if(zeta.gt.zetalim) then
L_lim = HIN/zetalim
else
L_lim = L
end if
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
else
call get_psi_coare(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
call get_psi_coare(psi_m,psi_h,HWIND/L,HTEMP/L)
endif
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
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_coare
subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT)
use sfx_sheba_noniterative, only: &
get_psi_stable_sheba => get_psi_stable
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
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
real Rib,wg,ksi,f, L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
L = HIN/zeta
if(zeta.gt.zetalim) then
L_lim = HIN/zetalim
else
L_lim = L
end if
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else ! unstable stratification (functions from ESM)
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
endif
else
! neutral stratification (z/L=0)
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
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_noit
subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
use sfx_most, only: &
get_psi_most => get_psi
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
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
real Rib,ksi,wg,f,L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
L = HIN/zeta
if(zeta.gt.zetalim) then
L_lim = HIN/zetalim
else
L_lim = L
end if
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L_lim)
call get_psi_most(psi_m,psi_h,HWIND/L_lim)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_most(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
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)
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L_lim,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.ne.0)then
L = HIN/zeta
if(zeta.gt.zetalim) then
L_lim = HIN/zetalim
else
L_lim = L
end if
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L_lim)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
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
real, intent(out) :: psi_m, psi_h
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 get_psi_esm1
! functions fim,fit were copied from the original pbldia
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