-
Notifications
You must be signed in to change notification settings - Fork 12
/
surface_varnce_module.F90
507 lines (397 loc) · 20.3 KB
/
surface_varnce_module.F90
1
2
3
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
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
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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
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
266
267
268
269
270
271
272
273
274
275
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
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
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
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
!-------------------------------------------------------------------------
! $Id$
!===============================================================================
module surface_varnce_module
implicit none
private ! Default to private
public :: calc_surface_varnce
contains
!=============================================================================
subroutine calc_surface_varnce( upwp_sfc, vpwp_sfc, wpthlp_sfc, wprtp_sfc, &
um_sfc, vm_sfc, Lscale_up_sfc, wpsclrp_sfc, &
wp2_splat_sfc, tau_zm_sfc, &
wp2_sfc, up2_sfc, vp2_sfc, &
thlp2_sfc, rtp2_sfc, rtpthlp_sfc, &
sclrp2_sfc, &
sclrprtp_sfc, &
sclrpthlp_sfc )
! Description:
! This subroutine computes estimate of the surface thermodynamic and wind
! component second order moments.
! References:
! Andre, J. C., G. De Moor, P. Lacarrere, G. Therry, and R. Du Vachat, 1978.
! Modeling the 24-Hour Evolution of the Mean and Turbulent Structures of
! the Planetary Boundary Layer. J. Atmos. Sci., 35, 1861 -- 1883.
!-----------------------------------------------------------------------
use parameters_model, only: &
T0 ! Variable(s)
use constants_clubb, only: &
four, & ! Variable(s)
two, &
one, &
two_thirds, &
one_third, &
one_fourth, &
zero, &
ten, &
grav, &
eps, &
w_tol_sqd, &
thl_tol, &
rt_tol, &
max_mag_correlation_flux, &
fstderr
use parameters_model, only: &
sclr_dim ! Variable(s)
use numerical_check, only: &
surface_varnce_check ! Procedure
use error_code, only: &
clubb_at_least_debug_level, & ! Procedure
err_code, & ! Error Indicator
clubb_fatal_error ! Constant
use array_index, only: &
iisclr_rt, & ! Index for a scalar emulating rt
iisclr_thl ! Index for a scalar emulating thetal
use stats_type_utilities, only: &
stat_end_update_pt, & ! Procedure(s)
stat_update_var_pt
use clubb_precision, only: &
core_rknd ! Variable(s)
use parameters_tunable, only: &
up2_vp2_factor ! Variable
implicit none
! External
intrinsic :: sqrt, max
! Constant Parameters
! Logical for Andre et al., 1978 parameterization.
logical, parameter :: l_andre_1978 = .false.
real( kind = core_rknd ), parameter :: &
a_const = 1.8_core_rknd, &
z_const = one, & ! Defined height of 1 meter [m]
! Vince Larson increased ufmin to stabilize arm_97. 24 Jul 2007
! ufmin = 0.0001_core_rknd, &
ufmin = 0.01_core_rknd, & ! Minimum allowable value of u* [m/s]
! End Vince Larson's change.
! Vince Larson changed in order to make correlations between [-1,1]. 31 Jan 2008.
! sclr_var_coef = 0.25_core_rknd, & ! This value is made up! - Vince Larson 12 Jul 2005
sclr_var_coef = 0.4_core_rknd, & ! This value is made up! - Vince Larson 12 Jul 2005
! End Vince Larson's change
! Vince Larson reduced surface spike in scalar variances associated
! w/ Andre et al. 1978 scheme
reduce_coef = 0.2_core_rknd
! Input Variables
real( kind = core_rknd ), intent(in) :: &
upwp_sfc, & ! Surface u momentum flux, <u'w'>|_sfc [m^2/s^2]
vpwp_sfc, & ! Surface v momentum flux, <v'w'>|_sfc [m^2/s^2]
wpthlp_sfc, & ! Surface thetal flux, <w'thl'>|_sfc [K m/s]
wprtp_sfc, & ! Surface moisture flux, <w'rt'>|_sfc [kg/kg m/s]
um_sfc, & ! Surface u wind component, <u> [m/s]
vm_sfc, & ! Surface v wind component, <v> [m/s]
Lscale_up_sfc, & ! Upward component of Lscale at surface [m]
wp2_splat_sfc, & ! Tendency of <w'^2> due to splatting of eddies at zm(1) [m^2/s^3]
tau_zm_sfc ! Turbulent dissipation time at level zm(1) [s]
real( kind = core_rknd ), intent(in), dimension(sclr_dim) :: &
wpsclrp_sfc ! Passive scalar flux, <w'sclr'>|_sfc [units m/s]
! Output Variables
real( kind = core_rknd ), intent(out) :: &
wp2_sfc, & ! Surface variance of w, <w'^2>|_sfc [m^2/s^2]
up2_sfc, & ! Surface variance of u, <u'^2>|_sfc [m^2/s^2]
vp2_sfc, & ! Surface variance of v, <v'^2>|_sfc [m^2/s^2]
thlp2_sfc, & ! Surface variance of theta-l, <thl'^2>|_sfc [K^2]
rtp2_sfc, & ! Surface variance of rt, <rt'^2>|_sfc [(kg/kg)^2]
rtpthlp_sfc ! Surface covariance of rt and theta-l [kg K/kg]
real( kind = core_rknd ), intent(out), dimension(sclr_dim) :: &
sclrp2_sfc, & ! Surface variance of passive scalar [units^2]
sclrprtp_sfc, & ! Surface covariance of pssv scalar and rt [units kg/kg]
sclrpthlp_sfc ! Surface covariance of pssv scalar and theta-l [units K]
! Local Variables
real( kind = core_rknd ) :: &
ustar2, & ! Square of surface friction velocity, u* [m^2/s^2]
wstar, & ! Convective velocity, w* [m/s]
uf, & ! Surface friction vel., u*, in older version [m/s]
min_wp2_sfc_val ! Minimum value of wp2_sfc that guarantees [m^2/s^2]
! correlation of (w,rt) and (w,thl) is within (-1,1)
! Variables for Andre et al., 1978 parameterization.
real( kind = core_rknd ) :: &
um_sfc_sqd, & ! Surface value of <u>^2 [m^2/s^2]
vm_sfc_sqd, & ! Surface value of <v>^2 [m^2/s^2]
usp2_sfc, & ! u_s (vector oriented w/ mean sfc. wind) variance [m^2/s^2]
vsp2_sfc, & ! v_s (vector perpen. to mean sfc. wind) variance [m^2/s^2]
wp2_splat_sfc_correction ! Reduction in wp2_sfc due to splatting [m^2/s^2]
real( kind = core_rknd ) :: &
ustar, & ! Surface friction velocity, u* [m/s]
Lngth, & ! Monin-Obukhov length [m]
zeta ! Dimensionless height z_const/Lngth, where z_const = 1 m. [-]
integer :: i ! Loop index
if ( l_andre_1978 ) then
! Calculate <u>^2 and <v>^2.
um_sfc_sqd = um_sfc**2
vm_sfc_sqd = vm_sfc**2
! Calculate surface friction velocity, u*.
ustar = max( ( upwp_sfc**2 + vpwp_sfc**2 )**(one_fourth), ufmin )
if ( abs(wpthlp_sfc) > eps) then
! Find Monin-Obukhov Length (Andre et al., 1978, p. 1866).
Lngth = - ( ustar**3 ) &
/ ( 0.35_core_rknd * (one/T0) * grav * wpthlp_sfc )
! Find the value of dimensionless height zeta
! (Andre et al., 1978, p. 1866).
zeta = z_const / Lngth
else ! wpthlp_sfc = 0
! The value of Monin-Obukhov length is +inf when ustar < 0 and -inf
! when ustar > 0. Either way, zeta = 0.
zeta = zero
endif ! wpthlp_sfc /= 0
! Andre et al, 1978, Eq. 29.
! Notes: 1) "reduce_coef" is a reduction coefficient intended to make
! the values of rtp2, thlp2, and rtpthlp smaller at the
! surface.
! 2) With the reduction coefficient having a value of 0.2, the
! surface correlations of both w & rt and w & thl have a value
! of about 0.845. These correlations are greater if zeta < 0.
! The correlations have a value greater than 1 if
! zeta <= -0.212.
! 3) The surface correlation of rt & thl is 1.
! Brian Griffin; February 2, 2008.
if ( zeta < zero ) then
thlp2_sfc = reduce_coef &
* ( wpthlp_sfc**2 / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
rtp2_sfc = reduce_coef &
* ( wprtp_sfc**2 / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
rtpthlp_sfc = reduce_coef &
* ( wprtp_sfc * wpthlp_sfc / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
wp2_sfc = ( ustar**2 ) &
* ( 1.75_core_rknd + two * (-zeta)**(two_thirds) )
else
thlp2_sfc = reduce_coef &
* four * ( wpthlp_sfc**2 / ustar**2 )
rtp2_sfc = reduce_coef &
* four * ( wprtp_sfc**2 / ustar**2 )
rtpthlp_sfc = reduce_coef &
* four * ( wprtp_sfc * wpthlp_sfc / ustar**2 )
wp2_sfc = 1.75_core_rknd * ustar**2
endif
thlp2_sfc = max( thl_tol**2, thlp2_sfc )
rtp2_sfc = max( rt_tol**2, rtp2_sfc )
! Calculate wstar following Andre et al., 1978, p. 1866.
! w* = ( ( 1 / T0 ) * g * <w'thl'>|_sfc * z_i )^(1/3);
! where z_i is the height of the mixed layer. The value of CLUBB's
! upward component of mixing length, Lscale_up, at the surface will be
! used as z_i.
wstar = ( (one/T0) * grav * wpthlp_sfc * Lscale_up_sfc )**(one_third)
! Andre et al., 1978, Eq. 29.
! Andre et al. (1978) defines horizontal wind surface variances in terms
! of orientation with the mean surface wind. The vector u_s is the wind
! vector oriented with the mean surface wind. The vector v_s is the wind
! vector oriented perpendicular to the mean surface wind. Thus, <u_s> is
! equal to the mean surface wind (both in speed and direction), and <v_s>
! is 0. Equation 29 gives the formula for the variance of u_s, which is
! <u_s'^2> (usp2_sfc in the code), and the formula for the variance of
! v_s, which is <v_s'^2> (vsp2_sfc in the code).
if ( wpthlp_sfc > zero ) then
usp2_sfc = four * ustar**2 + 0.3_core_rknd * wstar**2
vsp2_sfc = 1.75_core_rknd * ustar**2 + 0.3_core_rknd * wstar**2
else
usp2_sfc = four * ustar**2
vsp2_sfc = 1.75_core_rknd * ustar**2
endif
! Add effect of vertical compression of eddies on horizontal gustiness.
! First, ensure that wp2_sfc does not make the correlation
! of (w,thl) or (w,rt) outside (-1,1).
! Perhaps in the future we should also ensure that the correlations
! of (w,u) and (w,v) are not outside (-1,1).
min_wp2_sfc_val &
= max( w_tol_sqd, &
wprtp_sfc**2 / ( rtp2_sfc * max_mag_correlation_flux**2 ), &
wpthlp_sfc**2 / ( thlp2_sfc * max_mag_correlation_flux**2 ) )
if ( wp2_sfc + tau_zm_sfc * wp2_splat_sfc < min_wp2_sfc_val ) then
! splatting correction drives wp2_sfc to overly small value
wp2_splat_sfc_correction = -wp2_sfc + min_wp2_sfc_val
wp2_sfc = min_wp2_sfc_val
else
wp2_splat_sfc_correction = tau_zm_sfc * wp2_splat_sfc
wp2_sfc = wp2_sfc + wp2_splat_sfc_correction
end if
usp2_sfc = usp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
vsp2_sfc = vsp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
! Variance of u, <u'^2>, at the surface can be found from <u_s'^2>,
! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
! <u'^2>|_sfc = <u_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
! + <v_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
! where <u>^2 + <v>^2 /= 0.
up2_sfc &
= usp2_sfc * ( um_sfc_sqd / max( um_sfc_sqd + vm_sfc_sqd , eps ) ) &
+ vsp2_sfc * ( vm_sfc_sqd / max( um_sfc_sqd + vm_sfc_sqd , eps ) )
! Variance of v, <v'^2>, at the surface can be found from <u_s'^2>,
! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
! <v'^2>|_sfc = <v_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
! + <u_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
! where <u>^2 + <v>^2 /= 0.
vp2_sfc &
= vsp2_sfc * ( um_sfc_sqd / max( um_sfc_sqd + vm_sfc_sqd , eps ) ) &
+ usp2_sfc * ( vm_sfc_sqd / max( um_sfc_sqd + vm_sfc_sqd , eps ) )
! Passive scalars
if ( sclr_dim > 0 ) then
do i = 1, sclr_dim
! Notes: 1) "reduce_coef" is a reduction coefficient intended to
! make the values of sclrprtp, sclrpthlp, and sclrp2
! smaller at the surface.
! 2) With the reduction coefficient having a value of 0.2,
! the surface correlation of w & sclr has a value of
! about 0.845. The correlation is greater if zeta < 0.
! The correlation has a value greater than 1 if
! zeta <= -0.212.
! 3) The surface correlations of both rt & sclr and
! thl & sclr are 1.
! Brian Griffin; February 2, 2008.
if ( zeta < zero ) then
sclrprtp_sfc(i) &
= reduce_coef &
* ( wpsclrp_sfc(i) * wprtp_sfc / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
sclrpthlp_sfc(i) &
= reduce_coef &
* ( wpsclrp_sfc(i) * wpthlp_sfc / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
sclrp2_sfc(i) &
= reduce_coef &
* ( wpsclrp_sfc(i)**2 / ustar**2 ) &
* four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)
else
sclrprtp_sfc(i) &
= reduce_coef &
* four * ( wpsclrp_sfc(i) * wprtp_sfc / ustar**2 )
sclrpthlp_sfc(i) &
= reduce_coef &
* four * ( wpsclrp_sfc(i) * wpthlp_sfc / ustar**2 )
sclrp2_sfc(i) &
= reduce_coef &
* four * ( wpsclrp_sfc(i)**2 / ustar**2 )
endif
enddo ! i = 1, sclr_dim
endif
else ! Previous code.
! Compute ustar^2
ustar2 = sqrt( upwp_sfc * upwp_sfc + vpwp_sfc * vpwp_sfc )
! Compute wstar following Andre et al., 1976
if ( wpthlp_sfc > zero ) then
wstar = ( one/T0 * grav * wpthlp_sfc * z_const )**(one_third)
else
wstar = zero
endif
! Surface friction velocity following Andre et al. 1978
uf = sqrt( ustar2 + 0.3_core_rknd * wstar * wstar )
uf = max( ufmin, uf )
! Compute estimate for surface second order moments
wp2_sfc = a_const * uf**2
up2_sfc = up2_vp2_factor * a_const * uf**2 ! From Andre, et al. 1978
vp2_sfc = up2_vp2_factor * a_const * uf**2 ! " "
! Notes: 1) With "a" having a value of 1.8, the surface correlations of
! both w & rt and w & thl have a value of about 0.878.
! 2) The surface correlation of rt & thl is 0.5.
! Brian Griffin; February 2, 2008.
thlp2_sfc = 0.4_core_rknd * a_const * ( wpthlp_sfc / uf )**2
thlp2_sfc = max( thl_tol**2, thlp2_sfc )
rtp2_sfc = 0.4_core_rknd * a_const * ( wprtp_sfc / uf )**2
rtp2_sfc = max( rt_tol**2, rtp2_sfc )
rtpthlp_sfc = 0.2_core_rknd * a_const &
* ( wpthlp_sfc / uf ) * ( wprtp_sfc / uf )
! Add effect of vertical compression of eddies on horizontal gustiness.
! First, ensure that wp2_sfc does not make the correlation
! of (w,thl) or (w,rt) outside (-1,1).
! Perhaps in the future we should also ensure that the correlations
! of (w,u) and (w,v) are not outside (-1,1).
min_wp2_sfc_val &
= max( w_tol_sqd, &
wprtp_sfc**2 / ( rtp2_sfc * max_mag_correlation_flux**2 ), &
wpthlp_sfc**2 / ( thlp2_sfc * max_mag_correlation_flux**2 ) )
if ( wp2_sfc + tau_zm_sfc * wp2_splat_sfc < min_wp2_sfc_val ) then
! splatting correction drives wp2_sfc to overly small values
wp2_splat_sfc_correction = -wp2_sfc + min_wp2_sfc_val
wp2_sfc = min_wp2_sfc_val
else
wp2_splat_sfc_correction = tau_zm_sfc * wp2_splat_sfc
wp2_sfc = wp2_sfc + wp2_splat_sfc_correction
end if
up2_sfc = up2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
vp2_sfc = vp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
! Passive scalars
if ( sclr_dim > 0 ) then
do i = 1, sclr_dim
! Vince Larson changed coeffs to make correlations between [-1,1].
! 31 Jan 2008
! sclrprtp_sfc(i) &
! = a * (wprtp_sfc / uf) * (wpsclrp_sfc(i) / uf)
! sclrpthlp_sfc(i) &
! = a * (wpthlp_sfc / uf) * (wpsclrp_sfc(i) / uf)
! sclrp2_sfc(i) &
! = sclr_var_coef * a * ( wpsclrp_sfc(i) / uf )**2
! Notes: 1) With "a" having a value of 1.8 and "sclr_var_coef"
! having a value of 0.4, the surface correlation of
! w & sclr has a value of about 0.878.
! 2) With "sclr_var_coef" having a value of 0.4, the
! surface correlations of both rt & sclr and
! thl & sclr are 0.5.
! Brian Griffin; February 2, 2008.
! We use the following if..then's to make sclr_rt and sclr_thl
! close to the actual thlp2/rtp2 at the surface.
! -dschanen 25 Sep 08
if ( i == iisclr_rt ) then
! If we are trying to emulate rt with the scalar, then we
! use the variance coefficient from above
sclrprtp_sfc(i) = 0.4_core_rknd * a_const &
* ( wprtp_sfc / uf ) * ( wpsclrp_sfc(i) / uf )
else
sclrprtp_sfc(i) = 0.2_core_rknd * a_const &
* ( wprtp_sfc / uf ) * ( wpsclrp_sfc(i) / uf )
endif
if ( i == iisclr_thl ) then
! As above, but for thetal
sclrpthlp_sfc(i) = 0.4_core_rknd * a_const &
* ( wpthlp_sfc / uf ) &
* ( wpsclrp_sfc(i) / uf )
else
sclrpthlp_sfc(i) = 0.2_core_rknd * a_const &
* ( wpthlp_sfc / uf ) &
* ( wpsclrp_sfc(i) / uf )
endif
sclrp2_sfc(i) = sclr_var_coef * a_const &
* ( wpsclrp_sfc(i) / uf )**2
! End Vince Larson's change.
enddo ! 1,...sclr_dim
endif ! sclr_dim > 0
endif ! l_andre_1978
if ( clubb_at_least_debug_level( 2 ) ) then
call surface_varnce_check( wp2_sfc, up2_sfc, vp2_sfc, &
thlp2_sfc, rtp2_sfc, rtpthlp_sfc, &
sclrp2_sfc, sclrprtp_sfc, sclrpthlp_sfc )
if ( err_code == clubb_fatal_error ) then
write(fstderr,*) "Error in calc_surface_varnce"
write(fstderr,*) "Intent(in)"
write(fstderr,*) "upwp_sfc = ", upwp_sfc
write(fstderr,*) "vpwp_sfc = ", vpwp_sfc
write(fstderr,*) "wpthlp_sfc = ", wpthlp_sfc
write(fstderr,*) "wprtp_sfc = ", wprtp_sfc
if ( sclr_dim > 0 ) then
write(fstderr,*) "wpsclrp_sfc = ", wpsclrp_sfc
endif
write(fstderr,*) "Intent(out)"
write(fstderr,*) "wp2_sfc = ", wp2_sfc
write(fstderr,*) "up2_sfc = ", up2_sfc
write(fstderr,*) "vp2_sfc = ", vp2_sfc
write(fstderr,*) "thlp2_sfc = ", thlp2_sfc
write(fstderr,*) "rtp2_sfc = ", rtp2_sfc
write(fstderr,*) "rtpthlp_sfc = ", rtpthlp_sfc
if ( sclr_dim > 0 ) then
write(fstderr,*) "sclrp2_sfc = ", sclrp2_sfc
write(fstderr,*) "sclrprtp_sfc = ", sclrprtp_sfc
write(fstderr,*) "sclrpthlp_sfc = ", sclrpthlp_sfc
endif
endif ! err_code == clubb_fatal_error
endif ! clubb_at_least_debug_level ( 2 )
return
end subroutine calc_surface_varnce
!===============================================================================
end module surface_varnce_module