Newer
Older
module pbldia_new_sfx
implicit none
private
public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
pbldia_new_most,pbldia_new_esm
Anna Shestakova
committed
real,parameter,private::kappa=0.4,R=287.,appa=0.287
Anna Shestakova
committed
! 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
Anna Shestakova
committed
!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_sheba(AR2,ARDIN,ARDOUT)
use sfx_sheba, only: &
get_psi_sheba => get_psi
Anna Shestakova
committed
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !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
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
Anna Shestakova
committed
L = HIN/zeta
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)
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
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
Anna Shestakova
committed
return
end subroutine pbldia_new_sheba
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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 = 2. !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
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
if(zeta.gt.zetalim) zeta=zetalim
L = HIN/zeta
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_stable_sheba(psi_m,psi_h,HTEMP/L,HTEMP/L)
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
Anna Shestakova
committed
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !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
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
Anna Shestakova
committed
L = HIN/zeta
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)
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)
Anna Shestakova
committed
return
end subroutine pbldia_new_most
subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
Anna Shestakova
committed
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !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
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
Anna Shestakova
committed
L = HIN/zeta
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
else
UFWIND = ALOG(HWIND/HIN)
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)
Anna Shestakova
committed
return
end subroutine pbldia_new_esm
subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
Anna Shestakova
committed
Anna Shestakova
committed
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
Anna Shestakova
committed
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
Anna Shestakova
committed
if(zeta2 .ge. 0.) then
psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)
Anna Shestakova
committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
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