Newer
Older
module pbldia_new_sfx
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
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
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
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
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
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
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
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
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
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
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