Line data Source code
1 : !-------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module sfc_varnce_module
5 :
6 : implicit none
7 :
8 : private ! Default to private
9 :
10 : public :: calc_sfc_varnce
11 :
12 : contains
13 :
14 : !=============================================================================
15 8935056 : subroutine calc_sfc_varnce( nz, ngrdcol, sclr_dim, sclr_idx, &
16 8935056 : gr, dt, sfc_elevation, &
17 8935056 : upwp_sfc, vpwp_sfc, wpthlp, wprtp_sfc, &
18 8935056 : um, vm, Lscale_up, wpsclrp_sfc, &
19 8935056 : lhs_splat_wp2, tau_zm, &
20 : !wp2_splat_sfc, tau_zm_sfc, &
21 : l_vary_convect_depth, &
22 8935056 : up2_sfc_coef, &
23 8935056 : a_const, &
24 : stats_metadata, &
25 8935056 : stats_zm, &
26 8935056 : wp2, up2, vp2, &
27 8935056 : thlp2, rtp2, rtpthlp, &
28 8935056 : sclrp2, sclrprtp, sclrpthlp )
29 :
30 : ! Description:
31 : ! This subroutine computes estimate of the surface thermodynamic and wind
32 : ! component second order moments.
33 :
34 : ! References:
35 : ! Andre, J. C., G. De Moor, P. Lacarrere, G. Therry, and R. Du Vachat, 1978.
36 : ! Modeling the 24-Hour Evolution of the Mean and Turbulent Structures of
37 : ! the Planetary Boundary Layer. J. Atmos. Sci., 35, 1861 -- 1883.
38 : !-----------------------------------------------------------------------
39 :
40 : use grid_class, only: &
41 : grid
42 :
43 : use parameters_model, only: &
44 : T0 ! Variable(s)
45 :
46 : use constants_clubb, only: &
47 : four, & ! Variable(s)
48 : two, &
49 : one, &
50 : two_thirds, &
51 : one_third, &
52 : one_fourth, &
53 : zero, &
54 : grav, &
55 : eps, &
56 : w_tol_sqd, &
57 : thl_tol, &
58 : rt_tol, &
59 : max_mag_correlation_flux, &
60 : fstderr, &
61 : wp2_max
62 :
63 : use numerical_check, only: &
64 : sfc_varnce_check ! Procedure
65 :
66 : use error_code, only: &
67 : clubb_at_least_debug_level, & ! Procedure
68 : err_code, & ! Error Indicator
69 : clubb_fatal_error ! Constant
70 :
71 : use array_index, only: &
72 : sclr_idx_type
73 :
74 : use stats_type, only: &
75 : stats ! Type
76 :
77 : use stats_type_utilities, only: &
78 : stat_end_update_pt, & ! Procedure(s)
79 : stat_begin_update_pt, &
80 : stat_update_var_pt
81 :
82 : use stats_variables, only: &
83 : stats_metadata_type
84 :
85 : use clubb_precision, only: &
86 : core_rknd ! Variable(s)
87 :
88 : implicit none
89 :
90 : ! Constant Parameters
91 :
92 : ! Logical for Andre et al., 1978 parameterization.
93 : logical, parameter :: l_andre_1978 = .false.
94 :
95 : real( kind = core_rknd ), parameter :: &
96 : z_const = one, & ! Defined height of 1 meter [m]
97 : ! Vince Larson increased ufmin to stabilize arm_97. 24 Jul 2007
98 : ! ufmin = 0.0001_core_rknd, &
99 : ufmin = 0.01_core_rknd, & ! Minimum allowable value of u* [m/s]
100 : ! End Vince Larson's change.
101 : ! Vince Larson changed in order to make correlations between [-1,1]. 31 Jan 2008.
102 : ! sclr_var_coef = 0.25_core_rknd, & ! This value is made up! - Vince Larson 12 Jul 2005
103 : sclr_var_coef = 0.4_core_rknd, & ! This value is made up! - Vince Larson 12 Jul 2005
104 : ! End Vince Larson's change
105 : ! Vince Larson reduced surface spike in scalar variances associated
106 : ! w/ Andre et al. 1978 scheme
107 : reduce_coef = 0.2_core_rknd
108 :
109 : !-------------------------- Input Variables --------------------------
110 : integer, intent(in) :: &
111 : nz, &
112 : ngrdcol, &
113 : sclr_dim
114 :
115 : type (sclr_idx_type), intent(in) :: &
116 : sclr_idx
117 :
118 : type(grid), intent(in) :: &
119 : gr
120 :
121 : real( kind = core_rknd ) :: &
122 : dt
123 :
124 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
125 : sfc_elevation, &
126 : upwp_sfc, & ! Surface u momentum flux, <u'w'>|_sfc [m^2/s^2]
127 : vpwp_sfc, & ! Surface v momentum flux, <v'w'>|_sfc [m^2/s^2]
128 : wprtp_sfc ! Surface moisture flux, <w'rt'>|_sfc [kg/kg m/s]
129 :
130 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
131 : wpthlp, & ! Surface thetal flux, <w'thl'>|_sfc [K m/s]
132 : um, & ! Surface u wind component, <u> [m/s]
133 : vm, & ! Surface v wind component, <v> [m/s]
134 : Lscale_up, & ! Upward component of Lscale at surface [m]
135 : lhs_splat_wp2, &
136 : !wp2_splat, & ! Tendency of <w'^2> due to splatting of eddies at zm(1) [m^2/s^3]
137 : tau_zm ! Turbulent dissipation time at level zm(1) [s]
138 :
139 :
140 : real( kind = core_rknd ), dimension(ngrdcol,sclr_dim), intent(in) :: &
141 : wpsclrp_sfc ! Passive scalar flux, <w'sclr'>|_sfc [units m/s]
142 :
143 : logical, intent(in) :: &
144 : l_vary_convect_depth
145 :
146 : real( kind = core_rknd ), dimension(ngrdcol) :: &
147 : up2_sfc_coef, & ! CLUBB tunable parameter up2_sfc_coef [-]
148 : a_const ! Coefficient in front of wp2, up2, and vp2
149 :
150 : type (stats_metadata_type), intent(in) :: &
151 : stats_metadata
152 :
153 : !-------------------------- InOut Variables --------------------------
154 : type (stats), target, intent(inout), dimension(ngrdcol) :: &
155 : stats_zm
156 :
157 : !-------------------------- Output Variables --------------------------
158 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
159 : wp2, & ! Surface variance of w, <w'^2>|_sfc [m^2/s^2]
160 : up2, & ! Surface variance of u, <u'^2>|_sfc [m^2/s^2]
161 : vp2, & ! Surface variance of v, <v'^2>|_sfc [m^2/s^2]
162 : thlp2, & ! Surface variance of theta-l, <thl'^2>|_sfc [K^2]
163 : rtp2, & ! Surface variance of rt, <rt'^2>|_sfc [(kg/kg)^2]
164 : rtpthlp ! Surface covariance of rt and theta-l [kg K/kg]
165 :
166 : real( kind = core_rknd ), dimension(ngrdcol,nz,sclr_dim), intent(inout) :: &
167 : sclrp2, & ! Surface variance of passive scalar [units^2]
168 : sclrprtp, & ! Surface covariance of pssv scalar and rt [units kg/kg]
169 : sclrpthlp ! Surface covariance of pssv scalar and theta-l [units K]
170 :
171 : !-------------------------- Local Variables --------------------------
172 : real( kind = core_rknd ), dimension(ngrdcol) :: &
173 17870112 : uf, & ! Surface friction vel., u*, in older version [m/s]
174 17870112 : depth_pos_wpthlp, & ! Thickness of the layer near the surface with wpthlp > 0 [m]
175 17870112 : min_wp2_sfc_val ! Minimum value of wp2_sfc that guarantees [m^2/s^2]
176 : ! correlation of (w,rt) and (w,thl) is within (-1,1)
177 :
178 : ! Variables for Andre et al., 1978 parameterization.
179 : real( kind = core_rknd ), dimension(ngrdcol) :: &
180 8935056 : um_sfc_sqd, & ! Surface value of <u>^2 [m^2/s^2]
181 8935056 : vm_sfc_sqd, & ! Surface value of <v>^2 [m^2/s^2]
182 8935056 : usp2_sfc, & ! u_s (vector oriented w/ mean sfc. wind) variance [m^2/s^2]
183 8935056 : vsp2_sfc, & ! v_s (vector perpen. to mean sfc. wind) variance [m^2/s^2]
184 8935056 : ustar, & ! Surface friction velocity, u* [m/s]
185 8935056 : zeta, & ! Dimensionless height z_const/Lngth, where z_const = 1 m. [-]
186 17870112 : wp2_splat_sfc_correction ! Reduction in wp2_sfc due to splatting [m^2/s^2]
187 :
188 : real( kind = core_rknd ) :: &
189 : ustar2, & ! Square of surface friction velocity, u* [m^2/s^2]
190 : wstar ! Convective velocity, w* [m/s]
191 :
192 : real( kind = core_rknd ) :: &
193 : Lngth ! Monin-Obukhov length [m]
194 :
195 : integer :: i, k, sclr ! Loop index
196 :
197 : !-------------------------- Begin Code --------------------------
198 :
199 : !$acc enter data create( uf, depth_pos_wpthlp, min_wp2_sfc_val, &
200 : !$acc um_sfc_sqd, vm_sfc_sqd, usp2_sfc, vsp2_sfc, &
201 : !$acc ustar, zeta, wp2_splat_sfc_correction )
202 :
203 : ! Reflect surface varnce changes in budget
204 8935056 : if ( stats_metadata%l_stats_samp ) then
205 :
206 : !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp )
207 :
208 0 : do i = 1, ngrdcol
209 : call stat_begin_update_pt( stats_metadata%ithlp2_sf, 1, & ! intent(in)
210 0 : thlp2(i,1) / dt, & ! intent(in)
211 0 : stats_zm(i) ) ! intent(inout)
212 :
213 : call stat_begin_update_pt( stats_metadata%irtp2_sf, 1, & ! intent(in)
214 : rtp2(i,1) / dt, & ! intent(in)
215 0 : stats_zm(i) ) ! intent(inout)
216 :
217 : call stat_begin_update_pt( stats_metadata%irtpthlp_sf, 1, & ! intent(in)
218 : rtpthlp(i,1) / dt, & ! intent(in)
219 0 : stats_zm(i) ) ! intent(inout)
220 :
221 : call stat_begin_update_pt( stats_metadata%iup2_sf, 1, & ! intent(in)
222 : up2(i,1) / dt, & ! intent(in)
223 0 : stats_zm(i) ) ! intent(inout)
224 :
225 : call stat_begin_update_pt( stats_metadata%ivp2_sf, 1, & ! intent(in)
226 : vp2(i,1) / dt, & ! intent(in)
227 0 : stats_zm(i) ) ! intent(inout)
228 :
229 : call stat_begin_update_pt( stats_metadata%iwp2_sf, 1, & ! intent(in)
230 : wp2(i,1) / dt, & ! intent(in)
231 0 : stats_zm(i) ) ! intent(inout)
232 : end do
233 : end if
234 :
235 : !$acc parallel loop gang vector default(present)
236 149194656 : do i = 1, ngrdcol
237 : ! Find thickness of layer near surface with positive heat flux.
238 : ! This is used when l_vary_convect_depth=.true. in order to determine wp2.
239 149194656 : if ( wpthlp(i,1) <= zero ) then
240 46511616 : depth_pos_wpthlp(i) = one ! When sfc heat flux is negative, set depth to 1 m.
241 : else ! When sfc heat flux is positive, march up sounding until wpthlp 1st becomes negative.
242 : k = 1
243 571111866 : do while ( wpthlp(i,k) > zero .and. (gr%zm(i,k)-sfc_elevation(i)) < 1000._core_rknd )
244 477363882 : k = k + 1
245 : end do
246 93747984 : depth_pos_wpthlp(i) = max( one, gr%zm(i,k)-sfc_elevation(i) )
247 : end if
248 : end do
249 : !$acc end parallel loop
250 :
251 : ! a_const used to be set here; now it is a tunable parameter
252 : !if ( .not. l_vary_convect_depth ) then
253 : ! a_const = 1.8_core_rknd
254 : !else
255 : ! a_const = 0.6_core_rknd
256 : !end if
257 :
258 : if ( l_andre_1978 ) then
259 :
260 : ! Calculate <u>^2 and <v>^2.
261 : !$acc parallel loop gang vector default(present)
262 : do i = 1, ngrdcol
263 : um_sfc_sqd(i) = um(i,2)**2
264 : vm_sfc_sqd(i) = vm(i,2)**2
265 : end do
266 : !$acc end parallel loop
267 :
268 : ! Calculate surface friction velocity, u*.
269 : !$acc parallel loop gang vector default(present)
270 : do i = 1, ngrdcol
271 : ustar(i) = max( ( upwp_sfc(i)**2 + vpwp_sfc(i)**2 )**(one_fourth), ufmin )
272 : end do
273 : !$acc end parallel loop
274 :
275 : !$acc parallel loop gang vector default(present)
276 : do i = 1, ngrdcol
277 : if ( abs(wpthlp(i,1)) > eps) then
278 :
279 : ! Find Monin-Obukhov Length (Andre et al., 1978, p. 1866).
280 : Lngth = - ( ustar(i)**3 ) &
281 : / ( 0.35_core_rknd * (one/T0) * grav * wpthlp(i,1) )
282 :
283 : ! Find the value of dimensionless height zeta
284 : ! (Andre et al., 1978, p. 1866).
285 : zeta(i) = z_const / Lngth
286 :
287 : else ! wpthlp = 0
288 :
289 : ! The value of Monin-Obukhov length is +inf when ustar < 0 and -inf
290 : ! when ustar > 0. Either way, zeta = 0.
291 : zeta(i) = zero
292 :
293 : endif ! wpthlp /= 0
294 : end do
295 : !$acc end parallel loop
296 :
297 : ! Andre et al, 1978, Eq. 29.
298 : ! Notes: 1) "reduce_coef" is a reduction coefficient intended to make
299 : ! the values of rtp2, thlp2, and rtpthlp smaller at the
300 : ! surface.
301 : ! 2) With the reduction coefficient having a value of 0.2, the
302 : ! surface correlations of both w & rt and w & thl have a value
303 : ! of about 0.845. These correlations are greater if zeta < 0.
304 : ! The correlations have a value greater than 1 if
305 : ! zeta <= -0.212.
306 : ! 3) The surface correlation of rt & thl is 1.
307 : ! Brian Griffin; February 2, 2008.
308 : !$acc parallel loop gang vector default(present)
309 : do i = 1, ngrdcol
310 : if ( zeta(i) < zero ) then
311 :
312 : thlp2(i,1) = reduce_coef &
313 : * ( wpthlp(i,1)**2 / ustar(i)**2 ) &
314 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
315 :
316 : rtp2(i,1) = reduce_coef &
317 : * ( wprtp_sfc(i)**2 / ustar(i)**2 ) &
318 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
319 :
320 : rtpthlp(i,1) = reduce_coef &
321 : * ( wprtp_sfc(i) * wpthlp(i,1) / ustar(i)**2 ) &
322 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
323 :
324 : wp2(i,1) = ( ustar(i)**2 ) &
325 : * ( 1.75_core_rknd + two * (-zeta(i))**(two_thirds) )
326 :
327 : else
328 :
329 : thlp2(i,1) = reduce_coef &
330 : * four * ( wpthlp(i,1)**2 / ustar(i)**2 )
331 :
332 : rtp2(i,1) = reduce_coef &
333 : * four * ( wprtp_sfc(i)**2 / ustar(i)**2 )
334 :
335 : rtpthlp(i,1) = reduce_coef &
336 : * four * ( wprtp_sfc(i) * wpthlp(i,1) / ustar(i)**2 )
337 :
338 : wp2(i,1) = 1.75_core_rknd * ustar(i)**2
339 :
340 : endif
341 : end do
342 : !$acc end parallel loop
343 :
344 : !$acc parallel loop gang vector default(present)
345 : do i = 1, ngrdcol
346 : thlp2(i,1) = max( thl_tol**2, thlp2(i,1) )
347 : rtp2(i,1) = max( rt_tol**2, rtp2(i,1) )
348 : end do
349 : !$acc end parallel loop
350 :
351 : ! Calculate wstar following Andre et al., 1978, p. 1866.
352 : ! w* = ( ( 1 / T0 ) * g * <w'thl'>|_sfc * z_i )^(1/3);
353 : ! where z_i is the height of the mixed layer. The value of CLUBB's
354 : ! upward component of mixing length, Lscale_up, at the surface will be
355 : ! used as z_i.
356 : !$acc parallel loop gang vector default(present)
357 : do i = 1, ngrdcol
358 : wstar = ( (one/T0) * grav * wpthlp(i,1) * Lscale_up(i,2) )**(one_third)
359 :
360 : ! Andre et al., 1978, Eq. 29.
361 : ! Andre et al. (1978) defines horizontal wind surface variances in terms
362 : ! of orientation with the mean surface wind. The vector u_s is the wind
363 : ! vector oriented with the mean surface wind. The vector v_s is the wind
364 : ! vector oriented perpendicular to the mean surface wind. Thus, <u_s> is
365 : ! equal to the mean surface wind (both in speed and direction), and <v_s>
366 : ! is 0. Equation 29 gives the formula for the variance of u_s, which is
367 : ! <u_s'^2> (usp2_sfc in the code), and the formula for the variance of
368 : ! v_s, which is <v_s'^2> (vsp2_sfc in the code).
369 : if ( wpthlp(i,1) > zero ) then
370 :
371 : usp2_sfc(i) = four * ustar(i)**2 + 0.3_core_rknd * wstar**2
372 :
373 : vsp2_sfc(i) = 1.75_core_rknd * ustar(i)**2 + 0.3_core_rknd * wstar**2
374 :
375 : else
376 :
377 : usp2_sfc(i) = four * ustar(i)**2
378 :
379 : vsp2_sfc(i) = 1.75_core_rknd * ustar(i)**2
380 :
381 : endif
382 : end do
383 : !$acc end parallel loop
384 :
385 : ! Add effect of vertical compression of eddies on horizontal gustiness.
386 : ! First, ensure that wp2 does not make the correlation
387 : ! of (w,thl) or (w,rt) outside (-1,1).
388 : ! Perhaps in the future we should also ensure that the correlations
389 : ! of (w,u) and (w,v) are not outside (-1,1).
390 : !$acc parallel loop gang vector default(present)
391 : do i = 1, ngrdcol
392 : min_wp2_sfc_val(i) &
393 : = max( w_tol_sqd, &
394 : wprtp_sfc(i)**2 / ( rtp2(i,1) * max_mag_correlation_flux**2 ), &
395 : wpthlp(i,1)**2 / ( thlp2(i,1) * max_mag_correlation_flux**2 ) )
396 : end do
397 : !$acc end parallel loop
398 :
399 : !$acc parallel loop gang vector default(present)
400 : do i = 1, ngrdcol
401 : if ( wp2(i,1) - tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1) &
402 : < min_wp2_sfc_val(i) ) then
403 : !if ( wp2(i,1) + tau_zm(i,1) * wp2_splat(i,1) < min_wp2_sfc_val(i) ) then
404 : ! splatting correction drives wp2 to overly small value
405 : wp2_splat_sfc_correction(i) = -wp2(i,1) + min_wp2_sfc_val(i)
406 : wp2(i,1) = min_wp2_sfc_val(i)
407 : else
408 : wp2_splat_sfc_correction(i) = tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1)
409 : !wp2_splat_sfc_correction(i) = tau_zm(i,1) * wp2_splat(i,1)
410 : wp2(i,1) = wp2(i,1) + wp2_splat_sfc_correction(i)
411 : end if
412 : end do
413 : !$acc end parallel loop
414 :
415 : !$acc parallel loop gang vector default(present)
416 : do i = 1, ngrdcol
417 : usp2_sfc(i) = usp2_sfc(i) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
418 : vsp2_sfc(i) = vsp2_sfc(i) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
419 : end do
420 : !$acc end parallel loop
421 :
422 : ! Variance of u, <u'^2>, at the surface can be found from <u_s'^2>,
423 : ! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
424 : ! <u'^2>|_sfc = <u_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
425 : ! + <v_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
426 : ! where <u>^2 + <v>^2 /= 0.
427 : !$acc parallel loop gang vector default(present)
428 : do i = 1, ngrdcol
429 : up2(i,1) = usp2_sfc(i) * ( um_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) ) &
430 : + vsp2_sfc(i) * ( vm_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )
431 : end do
432 : !$acc end parallel loop
433 :
434 : ! Variance of v, <v'^2>, at the surface can be found from <u_s'^2>,
435 : ! <v_s'^2>, and mean winds (at the surface) <u> and <v>, such that:
436 : ! <v'^2>|_sfc = <v_s'^2> * [ <u>^2 / ( <u>^2 + <v>^2 ) ]
437 : ! + <u_s'^2> * [ <v>^2 / ( <u>^2 + <v>^2 ) ];
438 : ! where <u>^2 + <v>^2 /= 0.
439 : !$acc parallel loop gang vector default(present)
440 : do i = 1, ngrdcol
441 : vp2(i,1) = vsp2_sfc(i) * ( um_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) ) &
442 : + usp2_sfc(i) * ( vm_sfc_sqd(i) / max( um_sfc_sqd(i) + vm_sfc_sqd(i) , eps ) )
443 : end do
444 : !$acc end parallel loop
445 :
446 : ! Passive scalars
447 : if ( sclr_dim > 0 ) then
448 :
449 : !$acc parallel loop gang vector collapse(2) default(present)
450 : do sclr = 1, sclr_dim
451 : do i = 1, ngrdcol
452 :
453 : ! Notes: 1) "reduce_coef" is a reduction coefficient intended to
454 : ! make the values of sclrprtp, sclrpthlp, and sclrp2
455 : ! smaller at the surface.
456 : ! 2) With the reduction coefficient having a value of 0.2,
457 : ! the surface correlation of w & sclr has a value of
458 : ! about 0.845. The correlation is greater if zeta < 0.
459 : ! The correlation has a value greater than 1 if
460 : ! zeta <= -0.212.
461 : ! 3) The surface correlations of both rt & sclr and
462 : ! thl & sclr are 1.
463 : ! Brian Griffin; February 2, 2008.
464 : if ( zeta(i) < zero ) then
465 :
466 : sclrprtp(i,1,sclr) &
467 : = reduce_coef &
468 : * ( wpsclrp_sfc(i,sclr) * wprtp_sfc(i) / ustar(i)**2 ) &
469 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
470 :
471 : sclrpthlp(i,1,sclr) &
472 : = reduce_coef &
473 : * ( wpsclrp_sfc(i,sclr) * wpthlp(i,1) / ustar(i)**2 ) &
474 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
475 :
476 : sclrp2(i,1,sclr) &
477 : = reduce_coef &
478 : * ( wpsclrp_sfc(i,sclr)**2 / ustar(i)**2 ) &
479 : * four * ( one - 8.3_core_rknd * zeta(i) )**(-two_thirds)
480 :
481 : else
482 :
483 : sclrprtp(i,1,sclr) &
484 : = reduce_coef &
485 : * four * ( wpsclrp_sfc(i,sclr) * wprtp_sfc(i) / ustar(i)**2 )
486 :
487 : sclrpthlp(i,1,sclr) &
488 : = reduce_coef &
489 : * four * ( wpsclrp_sfc(i,sclr) * wpthlp(i,1) / ustar(i)**2 )
490 :
491 : sclrp2(i,1,sclr) &
492 : = reduce_coef &
493 : * four * ( wpsclrp_sfc(i,sclr)**2 / ustar(i)**2 )
494 :
495 : endif
496 :
497 : end do
498 : enddo ! i = 1, sclr_dim
499 : !$acc end parallel loop
500 :
501 : endif
502 :
503 :
504 : else ! Previous code.
505 :
506 : ! Compute ustar^2
507 : !$acc parallel loop gang vector default(present)
508 149194656 : do i = 1, ngrdcol
509 140259600 : ustar2 = sqrt( upwp_sfc(i) * upwp_sfc(i) + vpwp_sfc(i) * vpwp_sfc(i) )
510 :
511 : ! Compute wstar following Andre et al., 1976
512 140259600 : if ( wpthlp(i,1) > zero ) then
513 93747984 : if ( .not. l_vary_convect_depth ) then
514 93747984 : wstar = ( one/T0 * grav * wpthlp(i,1) * z_const )**(one_third)
515 : else
516 0 : wstar = ( one/T0 * grav * wpthlp(i,1) * 0.2_core_rknd * depth_pos_wpthlp(i) )**(one_third)
517 : end if
518 : else
519 : wstar = zero
520 : endif
521 :
522 : ! Surface friction velocity following Andre et al. 1978
523 140259600 : if ( .not. l_vary_convect_depth ) then
524 140259600 : uf(i) = sqrt( ustar2 + 0.3_core_rknd * wstar * wstar )
525 : else
526 0 : uf(i) = sqrt( ustar2 + wstar * wstar )
527 : end if
528 :
529 149194656 : uf(i) = max( ufmin, uf(i) )
530 : end do
531 : !$acc end parallel loop
532 :
533 : ! Compute estimate for surface second order moments
534 : !$acc parallel loop gang vector default(present)
535 149194656 : do i = 1, ngrdcol
536 140259600 : wp2(i,1) = a_const(i) * uf(i)**2
537 140259600 : up2(i,1) = up2_sfc_coef(i) * a_const(i) * uf(i)**2 ! From Andre, et al. 1978
538 149194656 : vp2(i,1) = up2_sfc_coef(i) * a_const(i) * uf(i)**2 ! " "
539 : end do
540 : !$acc end parallel loop
541 :
542 : ! Notes: 1) With "a" having a value of 1.8, the surface correlations of
543 : ! both w & rt and w & thl have a value of about 0.878.
544 : ! 2) The surface correlation of rt & thl is 0.5.
545 : ! Brian Griffin; February 2, 2008.
546 :
547 8935056 : if ( .not. l_vary_convect_depth ) then
548 : !$acc parallel loop gang vector default(present)
549 149194656 : do i = 1, ngrdcol
550 140259600 : thlp2(i,1) = 0.4_core_rknd * a_const(i) * ( wpthlp(i,1) / uf(i) )**2
551 140259600 : rtp2(i,1) = 0.4_core_rknd * a_const(i) * ( wprtp_sfc(i) / uf(i) )**2
552 : rtpthlp(i,1) = 0.2_core_rknd * a_const(i) * ( wpthlp(i,1) / uf(i) ) &
553 149194656 : * ( wprtp_sfc(i) / uf(i) )
554 : end do
555 : !$acc end parallel loop
556 : else
557 : !$acc parallel loop gang vector default(present)
558 0 : do i = 1, ngrdcol
559 0 : thlp2(i,1) = ( wpthlp(i,1) / uf(i) )**2 / ( max_mag_correlation_flux**2 * a_const(i) )
560 0 : rtp2(i,1) = ( wprtp_sfc(i) / uf(i) )**2 / ( max_mag_correlation_flux**2 * a_const(i) )
561 0 : rtpthlp(i,1) = max_mag_correlation_flux * sqrt( thlp2(i,1) * rtp2(i,1) )
562 : end do
563 : !$acc end parallel loop
564 : end if
565 :
566 : !$acc parallel loop gang vector default(present)
567 149194656 : do i = 1, ngrdcol
568 140259600 : thlp2(i,1) = max( thl_tol**2, thlp2(i,1) )
569 149194656 : rtp2(i,1) = max( rt_tol**2, rtp2(i,1) )
570 : end do
571 : !$acc end parallel loop
572 :
573 : ! Add effect of vertical compression of eddies on horizontal gustiness.
574 : ! First, ensure that wp2 does not make the correlation
575 : ! of (w,thl) or (w,rt) outside (-1,1).
576 : ! Perhaps in the future we should also ensure that the correlations
577 : ! of (w,u) and (w,v) are not outside (-1,1).
578 : !$acc parallel loop gang vector default(present)
579 149194656 : do i = 1, ngrdcol
580 140259600 : min_wp2_sfc_val(i) &
581 : = max( w_tol_sqd, &
582 0 : wprtp_sfc(i)**2 / ( rtp2(i,1) * max_mag_correlation_flux**2 ), &
583 289454256 : wpthlp(i,1)**2 / ( thlp2(i,1) * max_mag_correlation_flux**2 ) )
584 : end do
585 : !$acc end parallel loop
586 :
587 : !$acc parallel loop gang vector default(present)
588 149194656 : do i = 1, ngrdcol
589 140259600 : if ( wp2(i,1) - tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1) &
590 : < min_wp2_sfc_val(i) ) then
591 : !if ( wp2(i,1) + tau_zm(i,1) * wp2_splat(i,1) < min_wp2_sfc_val(i) ) then
592 : ! splatting correction drives wp2 to overly small values
593 98548 : wp2_splat_sfc_correction(i) = -wp2(i,1) + min_wp2_sfc_val(i)
594 98548 : wp2(i,1) = min_wp2_sfc_val(i)
595 : else
596 140161052 : wp2_splat_sfc_correction(i) = tau_zm(i,1) * lhs_splat_wp2(i,1) * wp2(i,1)
597 : !wp2_splat_sfc_correction(i) = tau_zm(i,1) * wp2_splat(i,1)
598 140161052 : wp2(i,1) = wp2(i,1) + wp2_splat_sfc_correction(i)
599 : end if
600 :
601 140259600 : up2(i,1) = up2(i,1) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
602 149194656 : vp2(i,1) = vp2(i,1) - 0.5_core_rknd * wp2_splat_sfc_correction(i)
603 : end do
604 : !$acc end parallel loop
605 :
606 : ! Passive scalars
607 8935056 : if ( sclr_dim > 0 ) then
608 : !$acc parallel loop gang vector collapse(2) default(present)
609 0 : do sclr = 1, sclr_dim
610 0 : do i = 1, ngrdcol
611 : ! Vince Larson changed coeffs to make correlations between [-1,1].
612 : ! 31 Jan 2008
613 : ! sclrprtp(i) &
614 : ! = a * (wprtp_sfc / uf) * (wpsclrp(i) / uf)
615 : ! sclrpthlp(i) &
616 : ! = a * (wpthlp / uf) * (wpsclrp(i) / uf)
617 : ! sclrp2(i) &
618 : ! = sclr_var_coef * a * ( wpsclrp(i) / uf )**2
619 : ! Notes: 1) With "a" having a value of 1.8 and "sclr_var_coef"
620 : ! having a value of 0.4, the surface correlation of
621 : ! w & sclr has a value of about 0.878.
622 : ! 2) With "sclr_var_coef" having a value of 0.4, the
623 : ! surface correlations of both rt & sclr and
624 : ! thl & sclr are 0.5.
625 : ! Brian Griffin; February 2, 2008.
626 :
627 : ! We use the following if..then's to make sclr_rt and sclr_thl
628 : ! close to the actual thlp2/rtp2 at the surface.
629 : ! -dschanen 25 Sep 08
630 0 : if ( sclr == sclr_idx%iisclr_rt ) then
631 : ! If we are trying to emulate rt with the scalar, then we
632 : ! use the variance coefficient from above
633 0 : sclrprtp(i,1,sclr) = 0.4_core_rknd * a_const(i) &
634 0 : * ( wprtp_sfc(i) / uf(i) ) * ( wpsclrp_sfc(i,sclr) / uf(i) )
635 : else
636 0 : sclrprtp(i,1,sclr) = 0.2_core_rknd * a_const(i) &
637 0 : * ( wprtp_sfc(i) / uf(i) ) * ( wpsclrp_sfc(i,sclr) / uf(i) )
638 : endif
639 :
640 0 : if ( sclr == sclr_idx%iisclr_thl ) then
641 : ! As above, but for thetal
642 0 : sclrpthlp(i,1,sclr) = 0.4_core_rknd * a_const(i) &
643 : * ( wpthlp(i,1) / uf(i) ) &
644 0 : * ( wpsclrp_sfc(i,sclr) / uf(i) )
645 : else
646 0 : sclrpthlp(i,1,sclr) = 0.2_core_rknd * a_const(i) &
647 : * ( wpthlp(i,1) / uf(i) ) &
648 0 : * ( wpsclrp_sfc(i,sclr) / uf(i) )
649 : endif
650 :
651 0 : sclrp2(i,1,sclr) = sclr_var_coef * a_const(i) &
652 0 : * ( wpsclrp_sfc(i,sclr) / uf(i) )**2
653 :
654 : ! End Vince Larson's change.
655 :
656 : end do
657 : enddo ! 1,...sclr_dim
658 : !$acc end parallel loop
659 :
660 : endif ! sclr_dim > 0
661 :
662 : endif ! l_andre_1978
663 :
664 : ! Clip wp2 at wp2_max, same as in advance_wp2_wp3
665 : !$acc parallel loop gang vector default(present)
666 149194656 : do i = 1, ngrdcol
667 149194656 : wp2(i,1) = min( wp2(i,1), wp2_max )
668 : end do
669 : !$acc end parallel loop
670 :
671 : !$acc parallel loop gang vector default(present)
672 149194656 : do i = 1, ngrdcol
673 149194656 : if ( abs(gr%zm(i,1)-sfc_elevation(i)) > abs(gr%zm(i,1)+sfc_elevation(i))*eps/2 ) then
674 :
675 : ! Variances for cases where the lowest level is not at the surface.
676 : ! Eliminate surface effects on lowest level variances.
677 0 : wp2(i,1) = w_tol_sqd
678 0 : up2(i,1) = w_tol_sqd
679 0 : vp2(i,1) = w_tol_sqd
680 0 : thlp2(i,1) = thl_tol**2
681 0 : rtp2(i,1) = rt_tol**2
682 0 : rtpthlp(i,1) = 0.0_core_rknd
683 :
684 : end if ! gr%zm(1,1) == sfc_elevation
685 : end do
686 : !$acc end parallel loop
687 :
688 8935056 : if ( sclr_dim > 0 ) then
689 :
690 : !$acc parallel loop gang vector collapse(2) default(present)
691 0 : do sclr = 1, sclr_dim
692 0 : do i = 1, ngrdcol
693 0 : if ( abs(gr%zm(i,1)-sfc_elevation(i)) > abs(gr%zm(i,1)+sfc_elevation(i))*eps/2 ) then
694 :
695 : ! Variances for cases where the lowest level is not at the surface.
696 : ! Eliminate surface effects on lowest level variances.
697 0 : sclrp2(i,1,sclr) = 0.0_core_rknd
698 0 : sclrprtp(i,1,sclr) = 0.0_core_rknd
699 0 : sclrpthlp(i,1,sclr) = 0.0_core_rknd
700 : end if
701 : end do
702 : end do
703 : end if
704 :
705 8935056 : if ( stats_metadata%l_stats_samp ) then
706 :
707 : !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp )
708 :
709 0 : do i = 1, ngrdcol
710 : call stat_end_update_pt( stats_metadata%ithlp2_sf, 1, & ! intent(in)
711 0 : thlp2(i,1) / dt, & ! intent(in)
712 0 : stats_zm(i) ) ! intent(inout)
713 :
714 : call stat_end_update_pt( stats_metadata%irtp2_sf, 1, & ! intent(in)
715 : rtp2(i,1) / dt, & ! intent(in)
716 0 : stats_zm(i) ) ! intent(inout)
717 :
718 : call stat_end_update_pt( stats_metadata%irtpthlp_sf, 1, & ! intent(in)
719 : rtpthlp(i,1) / dt, & ! intent(in)
720 0 : stats_zm(i) ) ! intent(inout)
721 :
722 : call stat_end_update_pt( stats_metadata%iup2_sf, 1, & ! intent(in)
723 : up2(i,1) / dt, & ! intent(in)
724 0 : stats_zm(i) ) ! intent(inout)
725 :
726 : call stat_end_update_pt( stats_metadata%ivp2_sf, 1, & ! intent(in)
727 : vp2(i,1) / dt, & ! intent(in)
728 0 : stats_zm(i) ) ! intent(inout)
729 :
730 : call stat_end_update_pt( stats_metadata%iwp2_sf, 1, & ! intent(in)
731 : wp2(i,1) / dt, & ! intent(in)
732 0 : stats_zm(i) ) ! intent(inout)
733 : end do
734 : end if
735 :
736 8935056 : if ( clubb_at_least_debug_level( 2 ) ) then
737 :
738 : !$acc update host( wp2, up2, vp2, thlp2, rtp2, rtpthlp, &
739 : !$acc upwp_sfc, vpwp_sfc, wpthlp, wprtp_sfc )
740 :
741 : !$acc update host( sclrp2, sclrprtp, sclrpthlp ) if ( sclr_dim > 0 )
742 :
743 0 : do i = 1, ngrdcol
744 0 : call sfc_varnce_check( sclr_dim, wp2(i,1), up2(i,1), vp2(i,1), & ! intent(in)
745 : thlp2(i,1), rtp2(i,1), rtpthlp(i,1), & ! intent(in)
746 0 : sclrp2(i,1,:), sclrprtp(i,1,:), sclrpthlp(i,1,:) ) ! intent(in)
747 : end do
748 :
749 0 : if ( err_code == clubb_fatal_error ) then
750 :
751 0 : write(fstderr,*) "Error in calc_sfc_varnce"
752 0 : write(fstderr,*) "Intent(in)"
753 :
754 0 : write(fstderr,*) "upwp_sfc = ", upwp_sfc
755 0 : write(fstderr,*) "vpwp_sfc = ", vpwp_sfc
756 0 : write(fstderr,*) "wpthlp = ", wpthlp
757 0 : write(fstderr,*) "wprtp_sfc = ", wprtp_sfc
758 :
759 0 : if ( sclr_dim > 0 ) then
760 0 : write(fstderr,*) "wpsclrp = ", wpsclrp_sfc
761 : endif
762 :
763 0 : write(fstderr,*) "Intent(out)"
764 :
765 0 : write(fstderr,*) "wp2 = ", wp2
766 0 : write(fstderr,*) "up2 = ", up2
767 0 : write(fstderr,*) "vp2 = ", vp2
768 0 : write(fstderr,*) "thlp2 = ", thlp2
769 0 : write(fstderr,*) "rtp2 = ", rtp2
770 0 : write(fstderr,*) "rtpthlp = ", rtpthlp
771 :
772 0 : if ( sclr_dim > 0 ) then
773 0 : write(fstderr,*) "sclrp2 = ", sclrp2
774 0 : write(fstderr,*) "sclrprtp = ", sclrprtp
775 0 : write(fstderr,*) "sclrpthlp = ", sclrpthlp
776 : endif
777 :
778 : endif ! err_code == clubb_fatal_error
779 :
780 : endif ! clubb_at_least_debug_level ( 2 )! Update surface stats
781 :
782 : !$acc exit data delete( uf, depth_pos_wpthlp, min_wp2_sfc_val, &
783 : !$acc um_sfc_sqd, vm_sfc_sqd, usp2_sfc, vsp2_sfc, &
784 : !$acc ustar, zeta, wp2_splat_sfc_correction )
785 :
786 8935056 : return
787 :
788 : end subroutine calc_sfc_varnce
789 :
790 : !===============================================================================
791 :
792 : end module sfc_varnce_module
|