Line data Source code
1 :
2 : module advance_windm_edsclrm_module
3 :
4 : implicit none
5 :
6 : private ! Set Default Scope
7 :
8 : public :: advance_windm_edsclrm
9 :
10 : private :: windm_edsclrm_solve, &
11 : compute_uv_tndcy, &
12 : windm_edsclrm_lhs, &
13 : windm_edsclrm_rhs
14 :
15 :
16 : ! Private named constants to avoid string comparisons
17 : integer, parameter, private :: &
18 : windm_edsclrm_um = 1, & ! Named constant to handle um solves
19 : windm_edsclrm_vm = 2, & ! Named constant to handle vm solves
20 : windm_edsclrm_scalar = 3, & ! Named constant to handle scalar solves
21 : clip_upwp = 10, & ! Named constant for upwp clipping
22 : ! NOTE: This must be the same as the clip_upwp
23 : ! declared in clip_explicit!
24 : clip_vpwp = 11 ! Named constant for vpwp clipping
25 : ! NOTE: This must be the same as the clip_vpwp
26 : ! declared in clip_explicit!
27 :
28 : contains
29 :
30 : !=============================================================================
31 8935056 : subroutine advance_windm_edsclrm( nz, ngrdcol, edsclr_dim, gr, dt, &
32 8935056 : wm_zt, Km_zm, Kmh_zm, &
33 8935056 : ug, vg, um_ref, vm_ref, &
34 8935056 : wp2, up2, vp2, um_forcing, vm_forcing, &
35 8935056 : edsclrm_forcing, &
36 8935056 : rho_ds_zm, invrs_rho_ds_zt, &
37 8935056 : fcor, l_implemented, &
38 : nu_vert_res_dep, &
39 : tridiag_solve_method, &
40 : l_predict_upwp_vpwp, &
41 : l_upwind_xm_ma, &
42 : l_uv_nudge, &
43 : l_tke_aniso, &
44 : l_lmm_stepping, &
45 : l_linearize_pbl_winds, &
46 : order_xp2_xpyp, order_wp2_wp3, order_windm, &
47 : stats_metadata, &
48 8935056 : stats_zt, stats_zm, stats_sfc, &
49 8935056 : um, vm, edsclrm, &
50 8935056 : upwp, vpwp, wpedsclrp, &
51 8935056 : um_pert, vm_pert, upwp_pert, vpwp_pert )
52 :
53 : ! Description:
54 : ! Solves for both mean horizontal wind components, um and vm, and for the
55 : ! eddy-scalars (passive scalars that don't use the high-order closure).
56 :
57 : ! Uses the LAPACK tridiagonal solver subroutine with 2 + # of scalar(s)
58 : ! back substitutions (since the left hand side matrix is the same for all
59 : ! input variables).
60 :
61 : ! References:
62 : ! Eqn. 8 & 9 on p. 3545 of
63 : ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
64 : ! Method and Model Description'' Golaz, et al. (2002)
65 : ! JAS, Vol. 59, pp. 3540--3551.
66 : !-----------------------------------------------------------------------
67 :
68 : use grid_class, only: &
69 : grid, & ! Type
70 : zm2zt
71 :
72 : use parameters_model, only: &
73 : ts_nudge ! Variable(s)
74 :
75 : use parameters_tunable, only: &
76 : nu_vertical_res_dep ! Type(s)
77 :
78 : use clubb_precision, only: &
79 : core_rknd ! Variable(s)
80 :
81 : use stats_type_utilities, only: &
82 : stat_begin_update, & ! Subroutines
83 : stat_end_update, &
84 : stat_update_var
85 :
86 : use stats_variables, only: &
87 : stats_metadata_type
88 :
89 : use clip_explicit, only: &
90 : clip_covar ! Procedure(s)
91 :
92 : use error_code, only: &
93 : clubb_at_least_debug_level, & ! Procedure
94 : err_code, & ! Error Indicator
95 : clubb_fatal_error ! Constant
96 :
97 : use constants_clubb, only: &
98 : one_half, & ! Constant(s)
99 : zero, &
100 : fstderr, &
101 : eps
102 :
103 : use sponge_layer_damping, only: &
104 : uv_sponge_damp_settings, &
105 : uv_sponge_damp_profile, &
106 : sponge_damp_xm ! Procedure(s)
107 :
108 : use stats_type, only: stats ! Type
109 :
110 : use diffusion, only: &
111 : diffusion_zt_lhs ! Procedure(s)
112 :
113 : use mean_adv, only: &
114 : term_ma_zt_lhs ! Procedures
115 :
116 : use model_flags, only: &
117 : l_upwind_Kh_dp_term
118 :
119 : use advance_helper_module, only: &
120 : calc_xpwp
121 :
122 : implicit none
123 :
124 : ! ------------------------ Input Variables ------------------------
125 : integer, intent(in) :: &
126 : nz, &
127 : ngrdcol, &
128 : edsclr_dim
129 :
130 : type (grid), target, intent(in) :: &
131 : gr
132 :
133 : real( kind = core_rknd ), intent(in) :: &
134 : dt ! Model timestep [s]
135 :
136 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
137 : wm_zt, & ! w wind component on thermodynamic levels [m/s]
138 : Km_zm, & ! Eddy diffusivity of winds on momentum levs. [m^2/s]
139 : Kmh_zm, & ! Eddy diffusivity of themo on momentum levs. [m^s/s]
140 : ug, & ! u (west-to-east) geostrophic wind comp. [m/s]
141 : vg, & ! v (south-to-north) geostrophic wind comp. [m/s]
142 : um_ref, & ! Reference u wind component for nudging [m/s]
143 : vm_ref, & ! Reference v wind component for nudging [m/s]
144 : wp2, & ! w'^2 (momentum levels) [m^2/s^2]
145 : up2, & ! u'^2 (momentum levels) [m^2/s^2]
146 : vp2, & ! v'^2 (momentum levels) [m^2/s^2]
147 : um_forcing, & ! u forcing [m/s/s]
148 : vm_forcing, & ! v forcing [m/s/s]
149 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
150 : invrs_rho_ds_zt ! Inv. dry, static density at thermo. levels [m^3/kg]
151 :
152 : real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(in) :: &
153 : edsclrm_forcing ! Eddy scalar large-scale forcing [{units vary}/s]
154 :
155 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
156 : fcor ! Coriolis parameter [s^-1]
157 :
158 : logical, intent(in) :: &
159 : l_implemented ! Flag for CLUBB being implemented in a larger model.
160 :
161 : type(nu_vertical_res_dep), intent(in) :: &
162 : nu_vert_res_dep ! Vertical resolution dependent nu values
163 :
164 : integer, intent(in) :: &
165 : tridiag_solve_method ! Specifier for method to solve tridiagonal systems
166 :
167 : logical, intent(in) :: &
168 : l_predict_upwp_vpwp, & ! Flag to predict <u'w'> and <v'w'> along with <u> and <v> alongside
169 : ! the advancement of <rt>, <w'rt'>, <thl>, <wpthlp>, <sclr>, and
170 : ! <w'sclr'> in subroutine advance_xm_wpxp. Otherwise, <u'w'> and
171 : ! <v'w'> are still approximated by eddy diffusivity when <u> and <v>
172 : ! are advanced in subroutine advance_windm_edsclrm.
173 : l_upwind_xm_ma, & ! This flag determines whether we want to use an upwind differencing
174 : ! approximation rather than a centered differencing for turbulent or
175 : ! mean advection terms. It affects rtm, thlm, sclrm, um and vm.
176 : l_uv_nudge, & ! For wind speed nudging
177 : l_tke_aniso, & ! For anisotropic turbulent kinetic energy, i.e. TKE = 1/2
178 : ! (u'^2 + v'^2 + w'^2)
179 : l_lmm_stepping, & ! Apply Linear Multistep Method (LMM) Stepping
180 : l_linearize_pbl_winds ! Flag (used by E3SM) to linearize PBL winds
181 :
182 : integer, intent(in) :: &
183 : order_xp2_xpyp, &
184 : order_wp2_wp3, &
185 : order_windm
186 :
187 : type (stats_metadata_type), intent(in) :: &
188 : stats_metadata
189 :
190 : ! ------------------------ Input/Output Variables ------------------------
191 : type (stats), dimension(ngrdcol), target, intent(inout) :: &
192 : stats_zt, &
193 : stats_zm, &
194 : stats_sfc
195 :
196 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
197 : um, & ! Mean u (west-to-east) wind component [m/s]
198 : vm, & ! Mean v (south-to-north) wind component [m/s]
199 : upwp, & ! <u'w'> (momentum levels) [m^2/s^2]
200 : vpwp ! <v'w'> (momentum levels) [m^2/s^2]
201 :
202 : ! Input/Output Variable for eddy-scalars
203 : real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(inout) :: &
204 : edsclrm ! Mean eddy scalar quantity [units vary]
205 :
206 : ! Output Variable for eddy-scalars
207 : real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim), intent(inout) :: &
208 : wpedsclrp ! w'edsclr' (momentum levels) [m/s {units vary}]
209 :
210 : ! Variables used to track perturbed version of winds.
211 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
212 : um_pert, & ! perturbed <u> [m/s]
213 : vm_pert, & ! perturbed <v> [m/s]
214 : upwp_pert, & ! perturbed <u'w'> [m^2/s^2]
215 : vpwp_pert ! perturbed <v'w'> [m^2/s^2]
216 :
217 : ! ------------------------ Local Variables ------------------------
218 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
219 17870112 : um_old, & ! Saved value of mean u (west-to-east) wind component [m/s]
220 17870112 : vm_old ! Saved value of Mean v (south-to-north) wind component [m/s]
221 :
222 : real( kind = core_rknd ), dimension(ngrdcol,nz,edsclr_dim) :: &
223 17870112 : edsclrm_old ! Saved value of mean eddy scalar quantity [units vary]
224 :
225 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
226 17870112 : um_tndcy, & ! u wind component tendency [m/s^2]
227 17870112 : vm_tndcy ! v wind component tendency [m/s^2]
228 :
229 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
230 17870112 : upwp_chnge, & ! Net change of u'w' due to clipping [m^2/s^2]
231 17870112 : vpwp_chnge ! Net change of v'w' due to clipping [m^2/s^2]
232 :
233 : real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
234 17870112 : lhs ! The implicit part of the tridiagonal matrix [units vary]
235 :
236 : real( kind = core_rknd ), dimension(ngrdcol,nz,max(2,edsclr_dim)) :: &
237 17870112 : rhs, &! The explicit part of the tridiagonal matrix [units vary]
238 17870112 : solution ! The solution to the tridiagonal matrix [units vary]
239 :
240 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
241 17870112 : wind_speed, & ! wind speed; sqrt(u^2 + v^2) [m/s]
242 17870112 : wind_speed_pert ! perturbed wind speed; sprt(u^2 + v^2) [m/s]
243 :
244 : real( kind = core_rknd ), dimension(ngrdcol) :: &
245 17870112 : u_star_sqd, & ! Surface friction velocity, u_star, squared [m/s]
246 17870112 : u_star_sqd_pert ! perturbed u_star, squared [m/s]
247 :
248 : logical :: &
249 : l_imp_sfc_momentum_flux ! Flag for implicit momentum surface fluxes.
250 :
251 : integer :: nrhs ! Number of right hand side terms
252 :
253 : integer :: i, k, edsclr, j ! Array index
254 :
255 : logical :: l_first_clip_ts, l_last_clip_ts ! flags for clip_covar
256 :
257 : real( kind = core_rknd ), dimension(ngrdcol) :: &
258 17870112 : nu_zero
259 :
260 : real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
261 17870112 : lhs_diff, & ! LHS diffustion terms
262 17870112 : lhs_ma_zt ! LHS mean advection terms
263 :
264 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
265 17870112 : Km_zt, & ! Eddy diffusivity of winds on momentum levs. [m^2/s]
266 17870112 : Kmh_zt, & ! Eddy diffusivity of themo on momentum levs. [m^s/s]
267 17870112 : Km_zm_p_nu10, & ! Km_zm plus nu_vert_res_dep%nu10
268 17870112 : xpwp ! x'w' for arbitrary x
269 :
270 : integer, parameter :: &
271 : kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
272 : k_tdiag = 2, & ! Thermodynamic main diagonal index.
273 : km1_tdiag = 3 ! Thermodynamic subdiagonal index.
274 :
275 : ! Whether perturbed winds are being solved.
276 : logical :: l_perturbed_wind
277 :
278 : real( kind = core_rknd ), dimension(nz) :: &
279 17870112 : tmp_in
280 :
281 : ! ------------------------ Begin Code ------------------------
282 :
283 : !$acc enter data create( um_old, vm_old, um_tndcy, vm_tndcy, &
284 : !$acc upwp_chnge, vpwp_chnge, lhs, rhs, solution, wind_speed, &
285 : !$acc wind_speed_pert, u_star_sqd, u_star_sqd_pert, &
286 : !$acc nu_zero, lhs_diff, lhs_ma_zt, Km_zt, Kmh_zt, &
287 : !$acc Km_zm_p_nu10, xpwp )
288 :
289 : !$acc enter data if( edsclr_dim > 0) create( edsclrm_old )
290 :
291 : !$acc parallel loop gang vector default(present)
292 149194656 : do i = 1, ngrdcol
293 149194656 : nu_zero(i) = zero
294 : end do
295 : !$acc end parallel loop
296 :
297 : !$acc parallel loop gang vector collapse(2) default(present)
298 768414816 : do k = 1, nz
299 12690480816 : do i = 1, ngrdcol
300 12681545760 : Km_zm_p_nu10(i,k) = Km_zm(i,k) + nu_vert_res_dep%nu10(i)
301 : end do
302 : end do
303 : !$acc end parallel loop
304 :
305 8935056 : l_perturbed_wind = ( .not. l_predict_upwp_vpwp ) .and. l_linearize_pbl_winds
306 :
307 8935056 : if ( .not. l_implemented ) then
308 : call term_ma_zt_lhs( nz, ngrdcol, wm_zt, gr%weights_zt2zm, & ! intent(in)
309 : gr%invrs_dzt, gr%invrs_dzm, & ! intent(in)
310 : l_upwind_xm_ma, & ! intent(in)
311 0 : lhs_ma_zt ) ! intent(out)
312 : else
313 : !$acc parallel loop gang vector collapse(3) default(present)
314 768414816 : do k = 1, nz
315 12690480816 : do i = 1, ngrdcol
316 48447743760 : do j = 1, 3
317 47688264000 : lhs_ma_zt(j,i,k) = zero
318 : end do
319 : end do
320 : end do
321 : !$acc end parallel loop
322 : end if
323 :
324 8935056 : if ( .not. l_predict_upwp_vpwp ) then
325 :
326 0 : Km_zt(:,:) = zm2zt( nz, ngrdcol, gr, Km_zm(:,:), zero )
327 :
328 : ! Calculate diffusion terms
329 : call diffusion_zt_lhs( nz, ngrdcol, gr, Km_zm, Km_zt, nu_vert_res_dep%nu10, & ! In
330 : invrs_rho_ds_zt, rho_ds_zm, & ! In
331 0 : lhs_diff ) ! Out
332 :
333 : ! The lower boundary condition needs to be applied here at level 2.
334 : if ( .not. l_upwind_Kh_dp_term ) then
335 : !$acc parallel loop gang vector default(present)
336 0 : do i = 1, ngrdcol
337 : ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
338 0 : lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
339 0 : * ( Km_zm(i,2) + nu_vert_res_dep%nu10(i) ) &
340 0 : * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
341 :
342 : ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
343 0 : lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
344 : * ( Km_zm(i,2) + nu_vert_res_dep%nu10(i) ) &
345 0 : * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
346 :
347 : ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
348 0 : lhs_diff(km1_tdiag,i,2) = zero
349 : end do
350 : !$acc end parallel loop
351 : else
352 : !$acc parallel loop gang vector default(present)
353 : do i = 1, ngrdcol
354 : ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
355 : lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) &
356 : * ( Km_zt(i,2) + nu_vert_res_dep%nu10(i) ) &
357 : * gr%invrs_dzm(i,2)
358 :
359 : ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
360 : lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) &
361 : * ( Km_zt(i,2) + nu_vert_res_dep%nu10(i) ) &
362 : * gr%invrs_dzm(i,2)
363 :
364 : ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
365 : lhs_diff(km1_tdiag,i,2) = zero
366 : end do
367 : !$acc end parallel loop
368 : end if
369 :
370 0 : if ( l_lmm_stepping ) then
371 : !$acc parallel loop gang vector collapse(2) default(present)
372 0 : do k = 1, nz
373 0 : do i = 1, ngrdcol
374 0 : um_old(i,k) = um(i,k)
375 0 : vm_old(i,k) = vm(i,k)
376 : end do
377 : end do
378 : !$acc end parallel loop
379 : end if ! l_lmm_stepping
380 :
381 : !----------------------------------------------------------------
382 : ! Prepare tridiagonal system for horizontal winds, um and vm
383 : !----------------------------------------------------------------
384 :
385 : ! Compute Coriolis, geostrophic, and other prescribed wind forcings for um.
386 : call compute_uv_tndcy( nz, ngrdcol, windm_edsclrm_um, & ! intent(in)
387 : fcor, vm, vg, & ! intent(in)
388 : um_forcing, l_implemented, & ! intent(in)
389 : stats_metadata, & ! intent(in)
390 : stats_zt, & ! intent(inout)
391 0 : um_tndcy ) ! intent(out)
392 :
393 : ! Compute Coriolis, geostrophic, and other prescribed wind forcings for vm.
394 : call compute_uv_tndcy( nz, ngrdcol, windm_edsclrm_vm, & ! intent(in)
395 : fcor, um, ug, & ! intent(in)
396 : vm_forcing, l_implemented, & ! intent(in)
397 : stats_metadata, & ! intent(in)
398 : stats_zt, & ! intent(inout)
399 0 : vm_tndcy ) ! intent(out)
400 :
401 : ! Momentum surface fluxes, u'w'|_sfc and v'w'|_sfc, are applied through
402 : ! an implicit method, such that:
403 : ! x'w'|_sfc = - ( u_star(t)^2 / wind_speed(t) ) * xm(t+1).
404 0 : l_imp_sfc_momentum_flux = .true.
405 :
406 : ! Compute wind speed (use threshold "eps" to prevent divide-by-zero
407 : ! error).
408 : !$acc parallel loop gang vector collapse(2) default(present)
409 0 : do k = 1, nz
410 0 : do i = 1, ngrdcol
411 0 : wind_speed(i,k) = max( sqrt( um(i,k)**2 + vm(i,k)**2 ), eps )
412 : end do
413 : end do
414 : !$acc end parallel loop
415 :
416 : ! Compute u_star_sqd according to the definition of u_star.
417 : !$acc parallel loop gang vector default(present)
418 0 : do i = 1, ngrdcol
419 0 : u_star_sqd(i) = sqrt( upwp(i,1)**2 + vpwp(i,1)**2 )
420 : end do
421 : !$acc end parallel loop
422 :
423 : ! Compute the explicit portion of the um equation.
424 : ! Build the right-hand side vector.
425 : call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_um, dt, & ! intent(in)
426 : lhs_diff, um, um_tndcy, & ! intent(in)
427 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
428 : l_imp_sfc_momentum_flux, upwp(:,1), & ! intent(in)
429 : stats_metadata, & ! intent(in)
430 : stats_zt, & ! intent(inout)
431 0 : rhs(:,:,windm_edsclrm_um) ) ! intent(out)
432 :
433 : ! Compute the explicit portion of the vm equation.
434 : ! Build the right-hand side vector.
435 : call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_vm, dt, & ! intent(in)
436 : lhs_diff, vm, vm_tndcy, & ! intent(in)
437 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
438 : l_imp_sfc_momentum_flux, vpwp(:,1), & ! intent(in)
439 : stats_metadata, & ! intent(in)
440 : stats_zt, & ! intent(inout)
441 0 : rhs(:,:,windm_edsclrm_vm) ) ! intent(out)
442 :
443 : ! Store momentum flux (explicit component)
444 :
445 : ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
446 : ! upwp(1) = upwp_sfc
447 : ! vpwp(1) = vpwp_sfc
448 :
449 : call calc_xpwp( nz, ngrdcol, gr, &
450 : Km_zm_p_nu10, um, &
451 0 : xpwp )
452 :
453 : ! Solve for x'w' at all intermediate model levels.
454 : ! A Crank-Nicholson timestep is used.
455 : !$acc parallel loop gang vector collapse(2) default(present)
456 0 : do k = 2, nz-1
457 0 : do i = 1, ngrdcol
458 0 : upwp(i,k) = -one_half * xpwp(i,k)
459 : end do
460 : end do
461 : !$acc end parallel loop
462 :
463 : call calc_xpwp( nz, ngrdcol, gr, &
464 : Km_zm_p_nu10, vm, &
465 0 : xpwp )
466 :
467 : !$acc parallel loop gang vector collapse(2) default(present)
468 0 : do k = 2, nz-1
469 0 : do i = 1, ngrdcol
470 0 : vpwp(i,k) = -one_half * xpwp(i,k)
471 : end do
472 : end do
473 : !$acc end parallel loop
474 :
475 : ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
476 : ! means that x'w' at the top model level is 0,
477 : ! since x'w' = - K_zm * d(xm)/dz.
478 : !$acc parallel loop gang vector default(present)
479 0 : do i = 1, ngrdcol
480 0 : upwp(i,nz) = zero
481 0 : vpwp(i,nz) = zero
482 : end do
483 : !$acc end parallel loop
484 :
485 : ! Compute the implicit portion of the um and vm equations.
486 : ! Build the left-hand side matrix.
487 : call windm_edsclrm_lhs( nz, ngrdcol, gr, dt, & ! intent(in)
488 : lhs_ma_zt, lhs_diff, & ! intent(in)
489 : wind_speed, u_star_sqd, & ! intent(in)
490 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
491 : l_implemented, l_imp_sfc_momentum_flux, & ! intent(in)
492 0 : lhs ) ! intent(out)
493 :
494 : ! Decompose and back substitute for um and vm
495 0 : nrhs = 2
496 : call windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, stats_metadata%iwindm_matrix_condt_num, & ! intent(in)
497 : tridiag_solve_method, & ! intent(in)
498 : stats_metadata, & ! intent(in)
499 : stats_sfc, & ! intent(inout)
500 : lhs, rhs, & ! intent(inout)
501 0 : solution ) ! intent(out)
502 :
503 : ! Check for singular matrices and bad LAPACK arguments
504 0 : if ( clubb_at_least_debug_level( 0 ) ) then
505 0 : if ( err_code == clubb_fatal_error ) then
506 0 : write(fstderr,*) "Fatal error solving for um/vm"
507 0 : return
508 : end if
509 : end if
510 :
511 : !----------------------------------------------------------------
512 : ! Update zonal (west-to-east) component of mean wind, um
513 : !----------------------------------------------------------------
514 : !$acc parallel loop gang vector collapse(2) default(present)
515 0 : do k = 1, nz
516 0 : do i = 1, ngrdcol
517 0 : um(i,k) = solution(i,k,windm_edsclrm_um)
518 : end do
519 : end do
520 : !$acc end parallel loop
521 :
522 : !----------------------------------------------------------------
523 : ! Update meridional (south-to-north) component of mean wind, vm
524 : !----------------------------------------------------------------
525 : !$acc parallel loop gang vector collapse(2) default(present)
526 0 : do k = 1, nz
527 0 : do i = 1, ngrdcol
528 0 : vm(i,k) = solution(i,k,windm_edsclrm_vm)
529 : end do
530 : end do
531 : !$acc end parallel loop
532 :
533 0 : if ( stats_metadata%l_stats_samp ) then
534 :
535 : !$acc update host( um, lhs_diff, lhs_ma_zt, invrs_rho_ds_zt, u_star_sqd, &
536 : !$acc rho_ds_zm, wind_speed, vm )
537 :
538 0 : do i = 1, ngrdcol
539 : ! Implicit contributions to um and vm
540 0 : call windm_edsclrm_implicit_stats( nz, windm_edsclrm_um, um(i,:), & ! intent(in)
541 0 : gr%invrs_dzt(i,:), & ! intent(in)
542 0 : lhs_diff(:,i,:), lhs_ma_zt(:,i,:), & ! intent(in)
543 0 : invrs_rho_ds_zt(i,:), u_star_sqd(i), & ! intent(in)
544 : rho_ds_zm(i,:), wind_speed(i,:), & ! intent(in)
545 : l_imp_sfc_momentum_flux, & ! intent(in)
546 : stats_metadata, & ! intent(in)
547 0 : stats_zt(i) ) ! intent(inout)
548 :
549 : call windm_edsclrm_implicit_stats( nz, windm_edsclrm_vm, vm(i,:), & ! intent(in)
550 0 : gr%invrs_dzt(i,:), & ! intent(in)
551 : lhs_diff(:,i,:), lhs_ma_zt(:,i,:), & ! intent(in)
552 : invrs_rho_ds_zt(i,:), u_star_sqd(i), & ! intent(in)
553 : rho_ds_zm(i,:), wind_speed(i,:), & ! intent(in)
554 : l_imp_sfc_momentum_flux, & ! intent(in)
555 : stats_metadata, & ! intent(in)
556 0 : stats_zt(i) ) ! intent(inout)
557 : end do
558 : end if ! stats_metadata%l_stats_samp
559 :
560 0 : if ( l_lmm_stepping ) then
561 : !$acc parallel loop gang vector collapse(2) default(present)
562 0 : do k = 1, nz
563 0 : do i = 1, ngrdcol
564 0 : um(i,k) = one_half * ( um_old(i,k) + um(i,k) )
565 0 : vm(i,k) = one_half * ( vm_old(i,k) + vm(i,k) )
566 : end do
567 : end do
568 : !$acc end parallel loop
569 : endif ! l_lmm_stepping ) then
570 :
571 0 : if ( uv_sponge_damp_settings%l_sponge_damping ) then
572 :
573 : !$acc update host( um, vm, um_ref, vm_ref )
574 :
575 : ! _sponge_damp_settings and _sponge_damp_profile could potentially vary
576 : ! from column to column, but there is no column index in those variables.
577 : ! Thus this code is potentially unsafe when implemented in a host model,
578 : ! which is indicated by l_implemented = T
579 0 : if ( l_implemented ) then
580 : write(fstderr,*) "l_sponge_damping = T and l_implemented = T &
581 0 : -- this is likely unsafe and considered fatal"
582 0 : err_code = clubb_fatal_error
583 0 : return
584 : end if
585 :
586 0 : if ( stats_metadata%l_stats_samp ) then
587 0 : do i = 1, ngrdcol
588 0 : tmp_in(1) = 0.0_core_rknd
589 0 : tmp_in(2:nz) = um(i,2:nz)
590 : call stat_begin_update( nz, stats_metadata%ium_sdmp, tmp_in / dt, & ! intent(in)
591 0 : stats_zt(i) ) ! intent(inout)
592 0 : tmp_in(2:nz) = vm(i,2:nz)
593 : call stat_begin_update( nz, stats_metadata%ivm_sdmp, tmp_in / dt, & ! intent(in)
594 0 : stats_zt(i) ) ! intent(inout)
595 : end do
596 : end if
597 :
598 0 : do i = 1, ngrdcol
599 0 : um(i,1:nz) = sponge_damp_xm( nz, dt, gr%zt(i,:), gr%zm(i,:), &
600 0 : um_ref(i,1:nz), um(i,1:nz), uv_sponge_damp_profile )
601 : end do
602 :
603 0 : do i = 1, ngrdcol
604 0 : vm(i,1:nz) = sponge_damp_xm( nz, dt, gr%zt(i,:), gr%zm(i,:), &
605 0 : vm_ref(i,1:nz), vm(i,1:nz), uv_sponge_damp_profile )
606 : end do
607 :
608 0 : if ( stats_metadata%l_stats_samp ) then
609 0 : do i = 1, ngrdcol
610 0 : tmp_in(1) = 0.0_core_rknd
611 0 : tmp_in(2:nz) = um(i,2:nz)
612 : call stat_end_update( nz, stats_metadata%ium_sdmp, tmp_in / dt, & ! intent(in)
613 0 : stats_zt(i) ) ! intent(inout)
614 0 : tmp_in(2:nz) = vm(i,2:nz)
615 : call stat_end_update( nz, stats_metadata%ivm_sdmp, tmp_in / dt, & ! intent(in)
616 0 : stats_zt(i) ) ! intent(inout)
617 : end do
618 : end if
619 :
620 : !$acc update device( um, vm )
621 :
622 : end if ! uv_sponge_damp_settings%l_sponge_damping
623 :
624 : ! Second part of momentum (implicit component)
625 :
626 : ! Solve for x'w' at all intermediate model levels.
627 : ! A Crank-Nicholson timestep is used.
628 : call calc_xpwp( nz, ngrdcol, gr, &
629 : Km_zm_p_nu10, um, &
630 0 : xpwp )
631 :
632 : !$acc parallel loop gang vector collapse(2) default(present)
633 0 : do k = 2, nz-1
634 0 : do i = 1, ngrdcol
635 0 : upwp(i,k) = upwp(i,k) - one_half * xpwp(i,k)
636 : end do
637 : end do
638 : !$acc end parallel loop
639 :
640 : call calc_xpwp( nz, ngrdcol, gr, &
641 : Km_zm_p_nu10, vm, &
642 0 : xpwp )
643 :
644 : !$acc parallel loop gang vector collapse(2) default(present)
645 0 : do k = 2, nz-1
646 0 : do i = 1, ngrdcol
647 0 : vpwp(i,k) = vpwp(i,k) - one_half * xpwp(i,k)
648 : end do
649 : end do
650 : !$acc end parallel loop
651 :
652 : ! Adjust um and vm if nudging is turned on.
653 0 : if ( l_uv_nudge ) then
654 :
655 : ! Reflect nudging in budget
656 0 : if ( stats_metadata%l_stats_samp ) then
657 : !$acc update host( um, vm )
658 0 : do i = 1, ngrdcol
659 0 : tmp_in(1) = 0.0_core_rknd
660 0 : tmp_in(2:nz) = um(i,2:nz)
661 : call stat_begin_update( nz, stats_metadata%ium_ndg, tmp_in / dt, & ! intent(in)
662 0 : stats_zt(i) ) ! intent(inout)
663 0 : tmp_in(2:nz) = vm(i,2:nz)
664 : call stat_begin_update( nz, stats_metadata%ivm_ndg, tmp_in / dt, & ! intent(in)
665 0 : stats_zt(i) ) ! intent(inout)
666 : end do
667 : end if
668 :
669 : !$acc parallel loop gang vector collapse(2) default(present)
670 0 : do k = 1, nz
671 0 : do i = 1, ngrdcol
672 0 : um(i,k) = um(i,k) - ( ( um(i,k) - um_ref(i,k) ) * (dt/ts_nudge) )
673 0 : vm(i,k) = vm(i,k) - ( ( vm(i,k) - vm_ref(i,k) ) * (dt/ts_nudge) )
674 : end do
675 : end do
676 : !$acc end parallel loop
677 :
678 0 : if ( stats_metadata%l_stats_samp ) then
679 : !$acc update host( um, vm )
680 0 : do i = 1, ngrdcol
681 0 : tmp_in(1) = 0.0_core_rknd
682 0 : tmp_in(2:nz) = um(i,2:nz)
683 : call stat_end_update( nz, stats_metadata%ium_ndg, tmp_in / dt, & ! intent(in)
684 0 : stats_zt(i) ) ! intent(inout)
685 0 : tmp_in(2:nz) = vm(i,2:nz)
686 : call stat_end_update( nz, stats_metadata%ivm_ndg, tmp_in / dt, & ! intent(in)
687 0 : stats_zt(i) ) ! intent(inout)
688 : end do
689 : end if
690 :
691 : end if ! l_uv_nudge
692 :
693 0 : if ( stats_metadata%l_stats_samp ) then
694 : !$acc update host( um_ref, vm_ref )
695 0 : do i = 1, ngrdcol
696 0 : tmp_in(1) = 0.0_core_rknd
697 0 : tmp_in(2:nz) = um_ref(i,2:nz)
698 : call stat_update_var( stats_metadata%ium_ref, tmp_in, & ! intent(in)
699 0 : stats_zt(i) ) ! intent(inout)
700 0 : tmp_in(2:nz) = vm_ref(i,2:nz)
701 : call stat_update_var( stats_metadata%ivm_ref, tmp_in, & ! intent(in)
702 0 : stats_zt(i) ) ! intent(inout)
703 : end do
704 : end if
705 :
706 : if ( order_windm < order_wp2_wp3 &
707 0 : .and. order_windm < order_xp2_xpyp ) then
708 0 : l_first_clip_ts = .true.
709 0 : l_last_clip_ts = .false.
710 : elseif ( order_windm > order_wp2_wp3 &
711 0 : .and. order_windm > order_xp2_xpyp ) then
712 0 : l_first_clip_ts = .false.
713 0 : l_last_clip_ts = .true.
714 : else
715 0 : l_first_clip_ts = .false.
716 0 : l_last_clip_ts = .false.
717 : endif
718 :
719 0 : if ( l_tke_aniso ) then
720 :
721 : ! Clipping for u'w'
722 : !
723 : ! Clipping u'w' at each vertical level, based on the
724 : ! correlation of u and w at each vertical level, such that:
725 : ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
726 : ! -1 <= corr_(u,w) <= 1.
727 : !
728 : ! Since u'^2, w'^2, and u'w' are each advanced in different
729 : ! subroutines from each other in advance_clubb_core, clipping for u'w'
730 : ! has to be done three times during each timestep (once after each
731 : ! variable has been updated).
732 : ! This is the third instance of u'w' clipping.
733 : !l_first_clip_ts = .false.
734 : !l_last_clip_ts = .true.
735 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
736 : l_last_clip_ts, dt, wp2, up2, & ! intent(in)
737 : l_predict_upwp_vpwp, & ! intent(in)
738 : stats_metadata, & ! intent(in)
739 : stats_zm, & ! intent(inout)
740 0 : upwp, upwp_chnge ) ! intent(inout)
741 :
742 : ! Clipping for v'w'
743 : !
744 : ! Clipping v'w' at each vertical level, based on the
745 : ! correlation of v and w at each vertical level, such that:
746 : ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
747 : ! -1 <= corr_(v,w) <= 1.
748 : !
749 : ! Since v'^2, w'^2, and v'w' are each advanced in different
750 : ! subroutines from each other in advance_clubb_core, clipping for v'w'
751 : ! has to be done three times during each timestep (once after each
752 : ! variable has been updated).
753 : ! This is the third instance of v'w' clipping.
754 : !l_first_clip_ts = .false.
755 : !l_last_clip_ts = .true.
756 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
757 : l_last_clip_ts, dt, wp2, vp2, & ! intent(in)
758 : l_predict_upwp_vpwp, & ! intent(in)
759 : stats_metadata, & ! intent(in)
760 : stats_zm, & ! intent(inout)
761 0 : vpwp, vpwp_chnge ) ! intent(inout)
762 : else
763 :
764 : ! intent(in) this case, it is assumed that
765 : ! u'^2 == v'^2 == w'^2, and the variables `up2' and `vp2' do not
766 : ! interact with any other variables.
767 : !l_first_clip_ts = .false.
768 : !l_last_clip_ts = .true.
769 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
770 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
771 : l_predict_upwp_vpwp, & ! intent(in)
772 : stats_metadata, & ! intent(in)
773 : stats_zm, & ! intent(inout)
774 0 : upwp, upwp_chnge ) ! intent(inout)
775 :
776 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
777 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
778 : l_predict_upwp_vpwp, & ! intent(in)
779 : stats_metadata, & ! intent(in)
780 : stats_zm, & ! intent(inout)
781 0 : vpwp, vpwp_chnge ) ! intent(inout)
782 : endif ! l_tke_aniso
783 :
784 : ! The values of um(1) and vm(1) are located below the model surface and
785 : ! do not affect the rest of the model. The values of um(1) or vm(1) are
786 : ! simply set to the values of um(2) and vm(2), respectively, after the
787 : ! equation matrices has been solved. Even though um and vm would sharply
788 : ! decrease to a value of 0 at the surface, this is done to avoid
789 : ! confusion on plots of the vertical profiles of um and vm.
790 : !$acc parallel loop gang vector default(present)
791 0 : do i = 1, ngrdcol
792 0 : um(i,1) = um(i,2)
793 0 : vm(i,1) = vm(i,2)
794 : end do
795 : !$acc end parallel loop
796 :
797 : endif ! .not. l_predict_upwp_vpwp
798 :
799 :
800 8935056 : if ( l_perturbed_wind ) then
801 :
802 : !----------------------------------------------------------------
803 : ! Prepare tridiagonal system for horizontal winds, um and vm
804 : !----------------------------------------------------------------
805 :
806 : ! Momentum surface fluxes, u'w'|_sfc and v'w'|_sfc, are applied through
807 : ! an implicit method, such that:
808 : ! x'w'|_sfc = - ( u_star(t)^2 / wind_speed(t) ) * xm(t+1).
809 0 : l_imp_sfc_momentum_flux = .true.
810 :
811 : ! Compute wind speed (use threshold "eps" to prevent divide-by-zero
812 : ! error).
813 : !$acc parallel loop gang vector collapse(2) default(present)
814 0 : do k = 1, nz
815 0 : do i = 1, ngrdcol
816 0 : wind_speed_pert(i,k) = max( sqrt( (um_pert(i,k))**2 + (vm_pert(i,k))**2 ), eps )
817 : end do
818 : end do
819 : !$acc end parallel loop
820 :
821 : ! Compute u_star_sqd according to the definition of u_star.
822 : !$acc parallel loop gang vector default(present)
823 0 : do i = 1, ngrdcol
824 0 : u_star_sqd_pert(i) = sqrt( upwp_pert(i,1)**2 + vpwp_pert(i,1)**2 )
825 : end do
826 : !$acc end parallel loop
827 :
828 : ! Compute the explicit portion of the um equation.
829 : ! Build the right-hand side vector.
830 : call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_um, dt, & ! intent(in)
831 : lhs_diff, um_pert, um_tndcy, & ! intent(in)
832 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
833 : l_imp_sfc_momentum_flux, upwp_pert(:,1), & ! intent(in)
834 : stats_metadata, & ! intent(in)
835 : stats_zt, & ! intent(inout)
836 0 : rhs(:,:,windm_edsclrm_um) ) ! intent(out)
837 :
838 : ! Compute the explicit portion of the vm equation.
839 : ! Build the right-hand side vector.
840 : call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_vm, dt, & ! intent(in)
841 : lhs_diff, vm_pert, vm_tndcy, & ! intent(in)
842 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
843 : l_imp_sfc_momentum_flux, vpwp_pert(:,1), & ! intent(in)
844 : stats_metadata, & ! intent(in)
845 : stats_zt, & ! intent(inout)
846 0 : rhs(:,:,windm_edsclrm_vm) ) ! intent(out)
847 :
848 : ! Store momentum flux (explicit component)
849 :
850 : ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
851 : ! upwp(1) = upwp_sfc
852 : ! vpwp(1) = vpwp_sfc
853 :
854 : ! Solve for x'w' at all intermediate model levels.
855 : ! A Crank-Nicholson timestep is used.
856 : call calc_xpwp( nz, ngrdcol, gr, &
857 : Km_zm_p_nu10, um_pert, &
858 0 : xpwp )
859 :
860 : !$acc parallel loop gang vector collapse(2) default(present)
861 0 : do k = 2, nz-1
862 0 : do i = 1, ngrdcol
863 0 : upwp_pert(i,k) = -one_half * xpwp(i,k)
864 : end do
865 : end do
866 : !$acc end parallel loop
867 :
868 : call calc_xpwp( nz, ngrdcol, gr, &
869 : Km_zm_p_nu10, vm_pert, &
870 0 : xpwp )
871 :
872 : !$acc parallel loop gang vector collapse(2) default(present)
873 0 : do k = 2, nz-1
874 0 : do i = 1, ngrdcol
875 0 : vpwp_pert(i,k) = -one_half * xpwp(i,k)
876 : end do
877 : end do
878 : !$acc end parallel loop
879 :
880 : ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
881 : ! means that x'w' at the top model level is 0,
882 : ! since x'w' = - K_zm * d(xm)/dz.
883 : !$acc parallel loop gang vector default(present)
884 0 : do i = 1, ngrdcol
885 0 : upwp_pert(i,nz) = zero
886 0 : vpwp_pert(i,nz) = zero
887 : end do
888 : !$acc end parallel loop
889 :
890 : ! Compute the implicit portion of the um and vm equations.
891 : ! Build the left-hand side matrix.
892 : call windm_edsclrm_lhs( nz, ngrdcol, gr, dt, & ! intent(in)
893 : lhs_ma_zt, lhs_diff, & ! intent(in)
894 : wind_speed_pert, u_star_sqd_pert, & ! intent(in)
895 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
896 : l_implemented, l_imp_sfc_momentum_flux, & ! intent(in)
897 0 : lhs ) ! intent(out)
898 :
899 : ! Decompose and back substitute for um and vm
900 0 : nrhs = 2
901 : call windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, stats_metadata%iwindm_matrix_condt_num, & ! intent(in)
902 : tridiag_solve_method, & ! intent(in)
903 : stats_metadata, & ! intent(in)
904 : stats_sfc, & ! intent(in)
905 : lhs, rhs, & ! intent(inout)
906 0 : solution ) ! intent(out)
907 :
908 : ! Check for singular matrices and bad LAPACK arguments
909 0 : if ( clubb_at_least_debug_level( 0 ) ) then
910 0 : if ( err_code == clubb_fatal_error ) then
911 0 : write(fstderr,*) "Fatal error solving for um_pert/vm_pert"
912 0 : return
913 : endif
914 : endif
915 :
916 : !----------------------------------------------------------------
917 : ! Update zonal (west-to-east) component of mean wind, um
918 : !----------------------------------------------------------------
919 : !$acc parallel loop gang vector collapse(2) default(present)
920 0 : do k = 1, nz
921 0 : do i = 1, ngrdcol
922 0 : um_pert(i,k) = solution(i,k,windm_edsclrm_um)
923 : end do
924 : end do
925 : !$acc end parallel loop
926 :
927 : !----------------------------------------------------------------
928 : ! Update meridional (south-to-north) component of mean wind, vm
929 : !----------------------------------------------------------------
930 : !$acc parallel loop gang vector collapse(2) default(present)
931 0 : do k = 1, nz
932 0 : do i = 1, ngrdcol
933 0 : vm_pert(i,k) = solution(i,k,windm_edsclrm_vm)
934 : end do
935 : end do
936 : !$acc end parallel loop
937 :
938 : ! Second part of momentum (implicit component)
939 :
940 : ! Solve for x'w' at all intermediate model levels.
941 : ! A Crank-Nicholson timestep is used.
942 : call calc_xpwp( nz, ngrdcol, gr, &
943 : Km_zm_p_nu10, um_pert, &
944 0 : xpwp )
945 :
946 : !$acc parallel loop gang vector collapse(2) default(present)
947 0 : do k = 2, nz-1
948 0 : do i = 1, ngrdcol
949 0 : upwp_pert(i,k) = upwp_pert(i,k) - one_half * xpwp(i,k)
950 : end do
951 : end do
952 : !$acc end parallel loop
953 :
954 : call calc_xpwp( nz, ngrdcol, gr, &
955 : Km_zm_p_nu10, vm_pert, &
956 0 : xpwp )
957 :
958 : !$acc parallel loop gang vector collapse(2) default(present)
959 0 : do k = 2, nz-1
960 0 : do i = 1, ngrdcol
961 0 : vpwp_pert(i,k) = vpwp_pert(i,k) - one_half * xpwp(i,k)
962 : end do
963 : end do
964 : !$acc end parallel loop
965 :
966 0 : if ( l_tke_aniso ) then
967 :
968 : ! Clipping for u'w'
969 : !
970 : ! Clipping u'w' at each vertical level, based on the
971 : ! correlation of u and w at each vertical level, such that:
972 : ! corr_(u,w) = u'w' / [ sqrt(u'^2) * sqrt(w'^2) ];
973 : ! -1 <= corr_(u,w) <= 1.
974 : !
975 : ! Since u'^2, w'^2, and u'w' are each advanced in different
976 : ! subroutines from each other in advance_clubb_core, clipping for u'w'
977 : ! has to be done three times during each timestep (once after each
978 : ! variable has been updated).
979 : ! This is the third instance of u'w' clipping.
980 0 : l_first_clip_ts = .false.
981 0 : l_last_clip_ts = .true.
982 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
983 : l_last_clip_ts, dt, wp2, up2, & ! intent(in)
984 : l_predict_upwp_vpwp, & ! intent(in)
985 : stats_metadata, & ! intent(in)
986 : stats_zm, & ! intent(inout)
987 0 : upwp_pert, upwp_chnge ) ! intent(inout)
988 :
989 : ! Clipping for v'w'
990 : !
991 : ! Clipping v'w' at each vertical level, based on the
992 : ! correlation of v and w at each vertical level, such that:
993 : ! corr_(v,w) = v'w' / [ sqrt(v'^2) * sqrt(w'^2) ];
994 : ! -1 <= corr_(v,w) <= 1.
995 : !
996 : ! Since v'^2, w'^2, and v'w' are each advanced in different
997 : ! subroutines from each other in advance_clubb_core, clipping for v'w'
998 : ! has to be done three times during each timestep (once after each
999 : ! variable has been updated).
1000 : ! This is the third instance of v'w' clipping.
1001 0 : l_first_clip_ts = .false.
1002 0 : l_last_clip_ts = .true.
1003 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
1004 : l_last_clip_ts, dt, wp2, vp2, & ! intent(in)
1005 : l_predict_upwp_vpwp, & ! intent(in)
1006 : stats_metadata, & ! intent(in)
1007 : stats_zm, & ! intent(inout)
1008 0 : vpwp_pert, vpwp_chnge ) ! intent(inout)
1009 : else
1010 :
1011 : ! intent(in) this case, it is assumed that
1012 : ! u'^2 == v'^2 == w'^2, and the variables `up2' and `vp2' do not
1013 : ! interact with any other variables.
1014 0 : l_first_clip_ts = .false.
1015 0 : l_last_clip_ts = .true.
1016 : call clip_covar( nz, ngrdcol, gr, clip_upwp, l_first_clip_ts, & ! intent(in)
1017 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
1018 : l_predict_upwp_vpwp, & ! intent(in)
1019 : stats_metadata, & ! intent(in)
1020 : stats_zm, & ! intent(inout)
1021 0 : upwp_pert, upwp_chnge ) ! intent(inout)
1022 :
1023 : call clip_covar( nz, ngrdcol, gr, clip_vpwp, l_first_clip_ts, & ! intent(in)
1024 : l_last_clip_ts, dt, wp2, wp2, & ! intent(in)
1025 : l_predict_upwp_vpwp, & ! intent(in)
1026 : stats_metadata, & ! intent(in)
1027 : stats_zm, & ! intent(inout)
1028 0 : vpwp_pert, vpwp_chnge ) ! intent(inout)
1029 :
1030 : end if ! l_tke_aniso
1031 :
1032 : ! The values of um(1) and vm(1) are located below the model surface and
1033 : ! do not affect the rest of the model. The values of um(1) or vm(1) are
1034 : ! simply set to the values of um(2) and vm(2), respectively, after the
1035 : ! equation matrices has been solved. Even though um and vm would sharply
1036 : ! decrease to a value of 0 at the surface, this is done to avoid
1037 : ! confusion on plots of the vertical profiles of um and vm.
1038 : !$acc parallel loop gang vector default(present)
1039 0 : do i = 1, ngrdcol
1040 0 : um_pert(i,1) = um_pert(i,2)
1041 0 : vm_pert(i,1) = vm_pert(i,2)
1042 : end do
1043 : !$acc end parallel loop
1044 :
1045 : end if ! l_perturbed_wind
1046 :
1047 : !----------------------------------------------------------------
1048 : ! Prepare tridiagonal system for eddy-scalars
1049 : !----------------------------------------------------------------
1050 :
1051 8935056 : if ( edsclr_dim > 0 ) then
1052 :
1053 8935056 : Kmh_zt(:,:) = zm2zt( nz, ngrdcol, gr, Kmh_zm(:,:), zero )
1054 :
1055 : ! Calculate diffusion terms
1056 : call diffusion_zt_lhs( nz, ngrdcol, gr, Kmh_zm, Kmh_zt, nu_zero, & ! intent(in)
1057 : invrs_rho_ds_zt, rho_ds_zm, & ! intent(in)
1058 8935056 : lhs_diff ) ! intent(out)
1059 :
1060 : ! The lower boundary condition needs to be applied here at level 2.
1061 : if ( .not. l_upwind_Kh_dp_term ) then
1062 :
1063 : !$acc parallel loop gang vector default(present)
1064 149194656 : do i = 1, ngrdcol
1065 : ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
1066 280519200 : lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
1067 : * ( Kmh_zm(i,2) + nu_zero(i) ) &
1068 280519200 : * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
1069 :
1070 : ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
1071 0 : lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
1072 : * ( Kmh_zm(i,2) + nu_zero(i) ) &
1073 140259600 : * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
1074 :
1075 : ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
1076 149194656 : lhs_diff(km1_tdiag,i,2) = zero
1077 : end do
1078 : !$acc end parallel loop
1079 :
1080 : else
1081 :
1082 : !$acc parallel loop gang vector default(present)
1083 : do i = 1, ngrdcol
1084 : ! Thermodynamic superdiagonal: [ x xm(k+1,<t+1>) ]
1085 : lhs_diff(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) &
1086 : * ( Kmh_zt(i,2) + nu_zero(i) ) * gr%invrs_dzm(i,2)
1087 :
1088 : ! Thermodynamic main diagonal: [ x xm(k,<t+1>) ]
1089 : lhs_diff(k_tdiag,i,2) = + gr%invrs_dzt(i,2) &
1090 : * ( Kmh_zt(i,2) + nu_zero(i) ) * gr%invrs_dzm(i,2)
1091 :
1092 : ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
1093 : lhs_diff(km1_tdiag,i,2) = zero
1094 : end do
1095 : !$acc end parallel loop
1096 :
1097 : end if
1098 :
1099 8935056 : if ( l_lmm_stepping ) then
1100 : !$acc parallel loop gang vector collapse(3) default(present)
1101 0 : do j = 1, edsclr_dim
1102 0 : do k = 1, nz
1103 0 : do i = 1, ngrdcol
1104 0 : edsclrm_old(i,k,j) = edsclrm(i,k,j)
1105 : end do
1106 : end do
1107 : end do
1108 : !$acc end parallel loop
1109 : endif ! l_lmm_stepping
1110 :
1111 : ! Eddy-scalar surface fluxes, x'w'|_sfc, are applied through an explicit
1112 : ! method.
1113 8935056 : l_imp_sfc_momentum_flux = .false.
1114 :
1115 : ! Compute the explicit portion of eddy scalar equation.
1116 : ! Build the right-hand side vector.
1117 : ! Because of statistics, we have to use a DO rather than a FORALL here
1118 : ! -dschanen 7 Oct 2008
1119 125090784 : do edsclr = 1, edsclr_dim
1120 : call windm_edsclrm_rhs( nz, ngrdcol, gr, windm_edsclrm_scalar, dt, & ! intent(in)
1121 : lhs_diff, edsclrm(:,:,edsclr), & ! intent(in)
1122 : edsclrm_forcing(:,:,edsclr), & ! intent(in)
1123 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
1124 : l_imp_sfc_momentum_flux, wpedsclrp(:,1,edsclr), & ! intent(in)
1125 : stats_metadata, & ! intent(in)
1126 : stats_zt, & ! intent(inout)
1127 125090784 : rhs(:,:,edsclr) ) ! intent(out)
1128 : enddo
1129 :
1130 :
1131 : ! Store momentum flux (explicit component)
1132 :
1133 : ! The surface flux, x'w'(1) = x'w'|_sfc, is set elsewhere in the model.
1134 : ! wpedsclrp(1,1:edsclr_dim) = wpedsclrp_sfc(1:edsclr_dim)
1135 :
1136 : ! Solve for x'w' at all intermediate model levels.
1137 : ! A Crank-Nicholson timestep is used
1138 125090784 : do edsclr = 1, edsclr_dim
1139 :
1140 : call calc_xpwp( nz, ngrdcol, gr, &
1141 : Km_zm_p_nu10, edsclrm(:,:,edsclr), &
1142 116155728 : xpwp )
1143 :
1144 : !$acc parallel loop gang vector collapse(2) default(present)
1145 9766016208 : do k = 2, nz-1
1146 >16109*10^7 : do i = 1, ngrdcol
1147 >16098*10^7 : wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
1148 : end do
1149 : end do
1150 : !$acc end parallel loop
1151 : end do
1152 :
1153 : ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
1154 : ! means that x'w' at the top model level is 0,
1155 : ! since x'w' = - K_zm * d(xm)/dz.
1156 : !$acc parallel loop gang vector collapse(2) default(present)
1157 125090784 : do j = 1, edsclr_dim
1158 1948465584 : do i = 1, ngrdcol
1159 1939530528 : wpedsclrp(i,nz,j) = zero
1160 : end do
1161 : end do
1162 : !$acc end parallel loop
1163 :
1164 : ! Compute the implicit portion of the xm (eddy-scalar) equations.
1165 : ! Build the left-hand side matrix.
1166 : call windm_edsclrm_lhs( nz, ngrdcol, gr, dt, & ! intent(in)
1167 : lhs_ma_zt, lhs_diff, & ! intent(in)
1168 : wind_speed, u_star_sqd, & ! intent(in)
1169 : rho_ds_zm, invrs_rho_ds_zt, & ! intent(in)
1170 : l_implemented, l_imp_sfc_momentum_flux, & ! intent(in)
1171 8935056 : lhs ) ! intent(out)
1172 :
1173 : ! Decompose and back substitute for all eddy-scalar variables
1174 : call windm_edsclrm_solve( nz, ngrdcol, gr, edsclr_dim, 0, & ! intent(in)
1175 : tridiag_solve_method, & ! intent(in)
1176 : stats_metadata, & ! intent(in)
1177 : stats_sfc, & ! intent(inout)
1178 : lhs, rhs, & ! intent(inout)
1179 8935056 : solution ) ! intent(out)
1180 :
1181 8935056 : if ( clubb_at_least_debug_level( 0 ) ) then
1182 8935056 : if ( err_code == clubb_fatal_error ) then
1183 0 : write(fstderr,*) "Fatal error solving for eddsclrm"
1184 : end if
1185 : end if
1186 :
1187 : !----------------------------------------------------------------
1188 : ! Update Eddy-diff. Passive Scalars
1189 : !----------------------------------------------------------------
1190 : !$acc parallel loop gang vector collapse(3) default(present)
1191 125090784 : do j = 1, edsclr_dim
1192 9998327664 : do k = 1, nz
1193 >16497*10^7 : do i = 1, ngrdcol
1194 >16486*10^7 : edsclrm(i,k,j) = solution(i,k,j)
1195 : end do
1196 : end do
1197 : end do
1198 : !$acc end parallel loop
1199 :
1200 8935056 : if ( l_lmm_stepping ) then
1201 : !$acc parallel loop gang vector collapse(3) default(present)
1202 0 : do j = 1, edsclr_dim
1203 0 : do k = 1, nz
1204 0 : do i = 1, ngrdcol
1205 0 : edsclrm(i,k,j) = one_half * ( edsclrm_old(i,k,j) + edsclrm(i,k,j) )
1206 : end do
1207 : end do
1208 : end do
1209 : !$acc end parallel loop
1210 : endif ! l_lmm_stepping
1211 :
1212 : ! Second part of momentum (implicit component)
1213 :
1214 : ! Solve for x'w' at all intermediate model levels.
1215 : ! A Crank-Nicholson timestep is used.
1216 125090784 : do edsclr = 1, edsclr_dim
1217 :
1218 : call calc_xpwp( nz, ngrdcol, gr, &
1219 : Kmh_zm, edsclrm(:,:,edsclr), &
1220 116155728 : xpwp )
1221 :
1222 : !$acc parallel loop gang vector collapse(2) default(present)
1223 9766016208 : do k = 2, nz-1
1224 >16109*10^7 : do i = 1, ngrdcol
1225 >16098*10^7 : wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
1226 : end do
1227 : end do
1228 : !$acc end parallel loop
1229 : end do
1230 :
1231 : ! Note that the w'edsclr' terms are not clipped, since we don't compute
1232 : ! the variance of edsclr anywhere. -dschanen 7 Oct 2008
1233 :
1234 : ! The value of edsclrm(1) is located below the model surface and does not
1235 : ! effect the rest of the model. The value of edsclrm(1) is simply set to
1236 : ! the value of edsclrm(2) after the equation matrix has been solved.
1237 : !$acc parallel loop gang vector collapse(2) default(present)
1238 125090784 : do edsclr = 1, edsclr_dim
1239 1948465584 : do i = 1, ngrdcol
1240 1939530528 : edsclrm(i,1,edsclr) = edsclrm(i,2,edsclr)
1241 : end do
1242 : end do
1243 : !$acc end parallel loop
1244 :
1245 : endif
1246 :
1247 8935056 : if ( clubb_at_least_debug_level( 0 ) ) then
1248 8935056 : if ( err_code == clubb_fatal_error ) then
1249 :
1250 : !$acc update host( wm_zt, Km_zm, ug, vg, um_ref, vm_ref, wp2, &
1251 : !$acc up2, vp2, um_forcing, vm_forcing, edsclrm_forcing, &
1252 : !$acc fcor, um_old, um, vm_old, vm, edsclrm_old, &
1253 : !$acc edsclrm, upwp, vpwp, wpedsclrp )
1254 :
1255 0 : write(fstderr,*) "Error in advance_windm_edsclrm"
1256 :
1257 0 : write(fstderr,*) "intent(in)"
1258 :
1259 0 : write(fstderr,*) "dt = ", dt
1260 0 : write(fstderr,*) "wm_zt = ", wm_zt
1261 0 : write(fstderr,*) "Km_zm = ", Km_zm
1262 0 : write(fstderr,*) "ug = ", ug
1263 0 : write(fstderr,*) "vg = ", vg
1264 0 : write(fstderr,*) "um_ref = ", um_ref
1265 0 : write(fstderr,*) "vm_ref = ", vm_ref
1266 0 : write(fstderr,*) "wp2 = ", wp2
1267 0 : write(fstderr,*) "up2 = ", up2
1268 0 : write(fstderr,*) "vp2 = ", vp2
1269 0 : write(fstderr,*) "um_forcing = ", um_forcing
1270 0 : write(fstderr,*) "vm_forcing = ", vm_forcing
1271 0 : do edsclr = 1, edsclr_dim
1272 0 : write(fstderr,*) "edsclrm_forcing # = ",edsclr, edsclrm_forcing
1273 : end do
1274 0 : write(fstderr,*) "fcor = ", fcor
1275 0 : write(fstderr,*) "l_implemented = ", l_implemented
1276 :
1277 0 : write(fstderr,*) "intent(inout)"
1278 :
1279 0 : if ( l_lmm_stepping ) &
1280 0 : write(fstderr,*) "um (pre-solve) = ", um_old
1281 0 : write(fstderr,*) "um = ", um
1282 0 : if ( l_lmm_stepping ) &
1283 0 : write(fstderr,*) "vm (pre-solve) = ", vm_old
1284 0 : write(fstderr,*) "vm = ", vm
1285 0 : do edsclr = 1, edsclr_dim
1286 0 : if ( l_lmm_stepping ) &
1287 0 : write(fstderr,*) "edsclrm (pre-solve) # ", edsclr, "=", edsclrm_old(:,:,edsclr)
1288 0 : write(fstderr,*) "edsclrm # ", edsclr, "=", edsclrm(:,:,edsclr)
1289 : end do
1290 0 : write(fstderr,*) "upwp = ", upwp
1291 0 : write(fstderr,*) "vpwp = ", vpwp
1292 0 : write(fstderr,*) "wpedsclrp = ", wpedsclrp
1293 :
1294 0 : return
1295 : end if
1296 : end if
1297 :
1298 : !$acc exit data delete( um_old, vm_old, um_tndcy, vm_tndcy, &
1299 : !$acc upwp_chnge, vpwp_chnge, lhs, rhs, solution, wind_speed, &
1300 : !$acc wind_speed_pert, u_star_sqd, u_star_sqd_pert, &
1301 : !$acc nu_zero, lhs_diff, lhs_ma_zt, Km_zt, Kmh_zt, &
1302 : !$acc Km_zm_p_nu10, xpwp )
1303 :
1304 : !$acc exit data if( edsclr_dim > 0) delete( edsclrm_old )
1305 :
1306 : return
1307 :
1308 : end subroutine advance_windm_edsclrm
1309 :
1310 : !=============================================================================
1311 8935056 : subroutine windm_edsclrm_solve( nz, ngrdcol, gr, nrhs, ixm_matrix_condt_num, &
1312 : tridiag_solve_method, &
1313 : stats_metadata, &
1314 8935056 : stats_sfc, &
1315 8935056 : lhs, rhs, solution )
1316 :
1317 : ! Note: In the "Description" section of this subroutine, the variable
1318 : ! "invrs_dzm" will be written as simply "dzm", and the variable
1319 : ! "invrs_dzt" will be written as simply "dzt". This is being done as
1320 : ! as device to save space and to make some parts of the description
1321 : ! more readable. This change does not pertain to the actual code.
1322 :
1323 : ! Description:
1324 : ! Solves the horizontal wind or eddy-scalar time-tendency equation, and
1325 : ! diagnoses the turbulent flux. A Crank-Nicholson time-stepping algorithm
1326 : ! is used in solving the turbulent advection term and in diagnosing the
1327 : ! turbulent flux.
1328 : !
1329 : ! The rate of change of an eddy-scalar quantity, xm, is:
1330 : !
1331 : ! d(xm)/dt = - w * d(xm)/dz - (1/rho_ds) * d( rho_ds * x'w' )/dz
1332 : ! + xm_forcings.
1333 : !
1334 : !
1335 : ! The Turbulent Advection Term
1336 : ! ----------------------------
1337 : !
1338 : ! The above equation contains a turbulent advection term:
1339 : !
1340 : ! - (1/rho_ds) * d( rho_ds * x'w' )/dz;
1341 : !
1342 : ! where the momentum flux, x'w', is closed using a down gradient approach:
1343 : !
1344 : ! x'w' = - K_zm * d(xm)/dz.
1345 : !
1346 : ! The turbulent advection term becomes:
1347 : !
1348 : ! + (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz;
1349 : !
1350 : ! which is the same as a standard eddy-diffusion term (if "rho_ds * K_zm" in
1351 : ! the term above is substituted for "K_zm" in a standard eddy-diffusion
1352 : ! term, and if the standard eddy-diffusion term is multiplied by
1353 : ! "1/rho_ds"). Thus, the turbulent advection term is treated and solved in
1354 : ! the same way that a standard eddy-diffusion term would be solved. The
1355 : ! term is discretized as follows:
1356 : !
1357 : ! The values of xm are found on the thermodynamic levels, while the values
1358 : ! of K_zm are found on the momentum levels. Additionally, the values of
1359 : ! rho_ds_zm are found on the momentum levels, and the values of
1360 : ! invrs_rho_ds_zt are found on the thermodynamic levels. The
1361 : ! derivatives (d/dz) of xm are taken over the intermediate momentum levels.
1362 : ! At the intermediate momentum levels, d(xm)/dz is multiplied by K_zm and by
1363 : ! rho_ds_zm. Then, the derivative of the whole mathematical expression is
1364 : ! taken over the central thermodynamic level, where it is multiplied by
1365 : ! invrs_rho_ds_zt, which yields the desired result.
1366 : !
1367 : ! ---xm(kp1)----------------------------------------------------- t(k+1)
1368 : !
1369 : ! ===========d(xm)/dz===K_zm(k)=====rho_ds_zm(k)================= m(k)
1370 : !
1371 : ! ---xm(k)---invrs_rho_ds_zt---d[rho_ds_zm*K_zm*d(xm)/dz]/dz----- t(k)
1372 : !
1373 : ! ===========d(xm)/dz===K_zm(km1)===rho_ds_zm(km1)=============== m(k-1)
1374 : !
1375 : ! ---xm(km1)----------------------------------------------------- t(k-1)
1376 : !
1377 : ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
1378 : ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
1379 : ! The letter "t" is used for thermodynamic levels and the letter "m" is used
1380 : ! for momentum levels.
1381 : !
1382 : ! dzt(k) = 1 / ( zm(k) - zm(k-1) )
1383 : ! dzm(k) = 1 / ( zt(k+1) - zt(k) )
1384 : ! dzm(k-1) = 1 / ( zt(k) - zt(k-1) )
1385 : !
1386 : ! The vertically discretized form of the turbulent advection term (treated
1387 : ! as an eddy diffusion term) is written out as:
1388 : !
1389 : ! + invrs_rho_ds_zt(k)
1390 : ! * dzt(k)
1391 : ! * [ rho_ds_zm(k) * K_zm(k) * dzm(k) * ( xm(k+1) - xm(k) )
1392 : ! - rho_ds_zm(k-1) * K_zm(k-1) * dzm(k-1) * ( xm(k) - xm(k-1) ) ].
1393 : !
1394 : ! For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme is
1395 : ! used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
1396 : ! eddy-diffusion term. The discretized implicit portion of the term is
1397 : ! written out as:
1398 : !
1399 : ! + (1/2) * invrs_rho_ds_zt(k)
1400 : ! * dzt(k)
1401 : ! * [ rho_ds_zm(k) * K_zm(k)
1402 : ! * dzm(k) * ( xm(k+1,<t+1>) - xm(k,<t+1>) )
1403 : ! - rho_ds_zm(k-1) * K_zm(k-1)
1404 : ! * dzm(k-1) * ( xm(k,<t+1>) - xm(k-1,<t+1>) ) ].
1405 : !
1406 : ! Note: When the implicit term is brought over to the left-hand side,
1407 : ! the sign is reversed and the leading "+" in front of the term
1408 : ! is changed to a "-".
1409 : !
1410 : ! The discretized explicit portion of the term is written out as:
1411 : !
1412 : ! + (1/2) * invrs_rho_ds_zt(k)
1413 : ! * dzt(k)
1414 : ! * [ rho_ds_zm(k) * K_zm(k)
1415 : ! * dzm(k) * ( xm(k+1,<t>) - xm(k,<t>) )
1416 : ! - rho_ds_zm(k-1) * K_zm(k-1)
1417 : ! * dzm(k-1) * ( xm(k,<t>) - xm(k-1,<t>) ) ].
1418 : !
1419 : ! Timestep index (t) stands for the index of the current timestep, while
1420 : ! timestep index (t+1) stands for the index of the next timestep, which is
1421 : ! being advanced to in solving the d(xm)/dt equation.
1422 : !
1423 : !
1424 : ! Boundary Conditions:
1425 : !
1426 : ! An eddy-scalar quantity is not allowed to flux out the upper boundary.
1427 : ! Thus, a zero-flux boundary condition is used for the upper boundary in the
1428 : ! eddy-diffusion equation.
1429 : !
1430 : ! The lower boundary condition is much more complicated. It is neither a
1431 : ! zero-flux nor a fixed-point boundary condition. Rather, it is a
1432 : ! fixed-flux boundary condition. This term is a turbulent advection term,
1433 : ! but with the eddy-scalars, the only value of x'w' relevant in solving the
1434 : ! d(xm)/dt equation is the value of x'w' at the surface (the first momentum
1435 : ! level), which is written as x'w'|_sfc.
1436 : !
1437 : ! 1) x'w' surface flux; generalized explicit form
1438 : !
1439 : ! The x'w' surface flux is applied to the d(xm)/dt equation through the
1440 : ! turbulent advection term, which is:
1441 : !
1442 : ! - (1/rho_ds) * d( rho_ds * x'w' )/dz.
1443 : !
1444 : ! At most vertical levels, a substitution can be made for x'w', such
1445 : ! that:
1446 : !
1447 : ! x'w' = - K_zm * d(xm)/dz.
1448 : !
1449 : ! However, the same substitution cannot be made at the surface (momentum
1450 : ! level 1), as x'w'|_sfc is a surface flux that is explicitly computed
1451 : ! elsewhere in the model code.
1452 : !
1453 : ! The lower boundary condition, which in this case needs to be applied to
1454 : ! the d(xm)/dt equation at level 2, is discretized as follows:
1455 : !
1456 : ! --xm(3)------------------------------------------------------- t(3)
1457 : !
1458 : ! ========[x'w'(2) = -K_zm(2)*d(xm)/dz]===rho_ds_zm(2)========== m(2)
1459 : !
1460 : ! --xm(2)---invrs_rho_ds_zt(2)---d[rho_ds_zm*K_zm*d(xm)/dz]/dz-- t(2)
1461 : !
1462 : ! ========[x'w'|_sfc]=====================rho_ds_zm(1)========== m(1) sfc
1463 : !
1464 : ! --xm(1)-------(below surface; not applicable)----------------- t(1)
1465 : !
1466 : ! where "sfc" is the level of the model surface or lower boundary.
1467 : !
1468 : ! The vertically discretized form of the turbulent advection term
1469 : ! (treated as an eddy diffusion term), with the explicit surface flux,
1470 : ! x'w'|_sfc, in place, is written out as:
1471 : !
1472 : ! - invrs_rho_ds_zt(2)
1473 : ! * dzt(2) * [ rho_ds_zm(2) * x'w'(2) - rho_ds_zm(1) * x'w'|_sfc ];
1474 : !
1475 : ! which can be re-written as:
1476 : !
1477 : ! + invrs_rho_ds_zt(2)
1478 : ! * dzt(2)
1479 : ! * [ rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
1480 : ! + rho_ds_zm(1) * x'w'|_sfc ];
1481 : !
1482 : ! which can be re-written again as:
1483 : !
1484 : ! + invrs_rho_ds_zt(2)
1485 : ! * dzt(2)
1486 : ! * rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
1487 : ! + invrs_rho_ds_zt(2)
1488 : ! * dzt(2)
1489 : ! * rho_ds_zm(1) * x'w'|_sfc.
1490 : !
1491 : ! For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme
1492 : ! is used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
1493 : ! eddy-diffusion term. The discretized implicit portion of the term is
1494 : ! written out as:
1495 : !
1496 : ! + (1/2) * invrs_rho_ds_zt(2)
1497 : ! * dzt(2)
1498 : ! * [ rho_ds_zm(2) * K_zm(2)
1499 : ! * dzm(2) * ( xm(3,<t+1>) - xm(2,<t+1>) ) ].
1500 : !
1501 : ! Note: When the implicit term is brought over to the left-hand side,
1502 : ! the sign is reversed and the leading "+" in front of the term
1503 : ! is changed to a "-".
1504 : !
1505 : ! The discretized explicit portion of the term is written out as:
1506 : !
1507 : ! + (1/2) * invrs_rho_ds_zt(2)
1508 : ! * dzt(2)
1509 : ! * [ rho_ds_zm(2) * K_zm(2)
1510 : ! * dzm(2) * ( xm(3,<t>) - xm(2,<t>) ) ]
1511 : ! + invrs_rho_ds_zt(2)
1512 : ! * dzt(2)
1513 : ! * rho_ds_zm(1) * x'w'|_sfc.
1514 : !
1515 : ! Note: The x'w'|_sfc portion of the term written above has been pulled
1516 : ! away from the rest of the explicit form written above because
1517 : ! the (1/2) factor due to Crank-Nicholson time_stepping does not
1518 : ! apply to it, as there isn't an implicit portion for x'w'|_sfc.
1519 : !
1520 : ! Timestep index (t) stands for the index of the current timestep, while
1521 : ! timestep index (t+1) stands for the index of the next timestep, which
1522 : ! is being advanced to in solving the d(xm)/dt equation.
1523 : !
1524 : ! 2) x'w' surface flux; implicit form for momentum fluxes u'w' and v'w'
1525 : !
1526 : ! The x'w' surface flux is applied to the d(xm)/dt equation through the
1527 : ! turbulent advection term, which is:
1528 : !
1529 : ! - (1/rho_ds) * d( rho_ds * x'w' )/dz.
1530 : !
1531 : ! At most vertical levels, a substitution can be made for x'w', such
1532 : ! that:
1533 : !
1534 : ! x'w' = - K_zm * d(xm)/dz.
1535 : !
1536 : ! However, the same substitution cannot be made at the surface (momentum
1537 : ! level 1), as x'w'|_sfc is a surface momentum flux that is found by the
1538 : ! following equation:
1539 : !
1540 : ! x'w'|_sfc = - [ u_star^2 / sqrt( um^2 + vm^2 ) ] * xm;
1541 : !
1542 : ! where x'w'|_sfc and xm are either u'w'|_sfc and um, respectively, or
1543 : ! v'w'|_sfc and vm, respectively (um and vm are located at the first
1544 : ! thermodynamic level above the surface, which is thermodynamic level 2),
1545 : ! sqrt( um^2 + vm^2 ) is the wind speed (also at thermodynamic level 2),
1546 : ! and u_star is defined as:
1547 : !
1548 : ! u_star = ( u'w'|_sfc^2 + v'w'|_sfc^2 )^(1/4);
1549 : !
1550 : ! and thus u_star^2 is defined as:
1551 : !
1552 : ! u_star^2 = sqrt( u'w'|_sfc^2 + v'w'|_sfc^2 ).
1553 : !
1554 : ! The value of u_star is either set to a constant value or computed
1555 : ! (through function diag_ustar) based on the surface wind speed, the
1556 : ! height above surface of the surface wind speed (as compared to the
1557 : ! roughness height), and the buoyancy flux at the surface. Either way,
1558 : ! u_star is computed elsewhere in the model, and the values of u'w'|_sfc
1559 : ! and v'w'|_sfc are based on it and computed along with it. The values
1560 : ! of u'w'|_sfc and v'w'|_sfc are then passed into advance_clubb_core,
1561 : ! and are eventually passed into advance_windm_edsclrm. In subroutine
1562 : ! advance_windm_edsclrm, the value of u_star_sqd is then recomputed
1563 : ! based on u'w'|_sfc and v'w'|_sfc. The value of sqrt( u_star_sqd ) is
1564 : ! consistent with the value of the original computation of u_star.
1565 : !
1566 : ! The equation listed above is substituted for x'w'|_sfc. The lower
1567 : ! boundary condition, which in this case needs to be applied to the
1568 : ! d(xm)/dt equation at level 2, is discretized as follows:
1569 : !
1570 : ! --xm(3)------------------------------------------------------- t(3)
1571 : !
1572 : ! ===[x'w'(2) = -K_zm(2)*d(xm)/dz]=================rho_ds_zm(2)= m(2)
1573 : !
1574 : ! --xm(2)---invrs_rho_ds_zt(2)---d[rho_ds_zm*K_zm*d(xm)/dz]/dz-- t(2)
1575 : !
1576 : ! ===[x'w'|_sfc = -[u_star^2/sqrt(um^2+vm^2)]*xm]==rho_ds_zm(1)= m(1) sfc
1577 : !
1578 : ! --xm(1)-------(below surface; not applicable)----------------- t(1)
1579 : !
1580 : ! where "sfc" is the level of the model surface or lower boundary.
1581 : !
1582 : ! The vertically discretized form of the turbulent advection term
1583 : ! (treated as an eddy diffusion term), with the implicit surface momentum
1584 : ! flux in place, is written out as:
1585 : !
1586 : ! - invrs_rho_ds_zt(2)
1587 : ! * dzt(2) * [ rho_ds_zm(2) * x'w'(2) - rho_ds_zm(1) * x'w'|_sfc ];
1588 : !
1589 : ! which can be re-written as:
1590 : !
1591 : ! - invrs_rho_ds_zt(2)
1592 : ! * dzt(2)
1593 : ! * [ rho_ds_zm(2)
1594 : ! * { - K_zm(2) * dzm(2) * ( xm(3) - xm(2) ) }
1595 : ! - rho_ds_zm(1)
1596 : ! * { - [ u_star^2 / sqrt( um(2)^2 + vm(2)^2 ) ] * xm(2) } ];
1597 : !
1598 : ! which can be re-written as:
1599 : !
1600 : ! + invrs_rho_ds_zt(2)
1601 : ! * dzt(2)
1602 : ! * rho_ds_zm(2) * K_zm(2) * dzm(2) * ( xm(3) - xm(2) )
1603 : ! - invrs_rho_ds_zt(2)
1604 : ! * dzt(2)
1605 : ! * rho_ds_zm(1) * [ u_star^2 / sqrt( um(2)^2 + vm(2)^2 ) ] * xm(2).
1606 : !
1607 : ! For this equation, a Crank-Nicholson (semi-implicit) diffusion scheme
1608 : ! is used to solve the (1/rho_ds) * d [ rho_ds * K_zm * d(xm)/dz ] / dz
1609 : ! eddy-diffusion term. The discretized implicit portion of the term is
1610 : ! written out as:
1611 : !
1612 : ! + (1/2) * invrs_rho_ds_zt(2)
1613 : ! * dzt(2)
1614 : ! * [ rho_ds_zm(2) * K_zm(2)
1615 : ! * dzm(2) * ( xm(3,<t+1>) - xm(2,<t+1>) ) ]
1616 : ! - invrs_rho_ds_zt(2)
1617 : ! * dzt(2)
1618 : ! * rho_ds_zm(1)
1619 : ! * [u_star^2/sqrt( um(2,<t>)^2 + vm(2,<t>)^2 )] * xm(2,<t+1>).
1620 : !
1621 : ! Note: When the implicit term is brought over to the left-hand side,
1622 : ! the signs are reversed and the leading "+" in front of the first
1623 : ! part of the term is changed to a "-", while the leading "-" in
1624 : ! front of the second part of the term is changed to a "+".
1625 : !
1626 : ! Note: The x'w'|_sfc portion of the term written above has been pulled
1627 : ! away from the rest of the implicit form written above because
1628 : ! the (1/2) factor due to Crank-Nicholson time_stepping does not
1629 : ! apply to it. The x'w'|_sfc portion of the term is treated
1630 : ! completely implicitly in order to enhance numerical stability.
1631 : !
1632 : ! The discretized explicit portion of the term is written out as:
1633 : !
1634 : ! + (1/2) * invrs_rho_ds_zt(2)
1635 : ! * dzt(2)
1636 : ! * [ rho_ds_zm(2) * K_zm(2)
1637 : ! * dzm(2) * ( xm(3,<t>) - xm(2,<t>) ) ].
1638 : !
1639 : ! Timestep index (t) stands for the index of the current timestep, while
1640 : ! timestep index (t+1) stands for the index of the next timestep, which
1641 : ! is being advanced to in solving the d(xm)/dt equation.
1642 : !
1643 : !
1644 : ! The lower boundary condition for the implicit and explicit portions of the
1645 : ! turbulent advection term, without the x'w'|_sfc portion of the term, can
1646 : ! easily be invoked by using the zero-flux boundary conditions found in the
1647 : ! generalized diffusion function (function diffusion_zt_lhs), which is used
1648 : ! for many other equations in this model. Either the generalized explicit
1649 : ! surface flux needs to be added onto the explicit term after the diffusion
1650 : ! function has been called from subroutine windm_edsclrm_rhs, or the
1651 : ! implicit momentum surface flux needs to be added onto the implicit term
1652 : ! after the diffusion function has been called from subroutine
1653 : ! windm_edsclrm_lhs. However, all other equations in this model that use
1654 : ! zero-flux diffusion have level 1 as the level to which the lower boundary
1655 : ! condition needs to be applied. Thus, an adjuster will have to be used at
1656 : ! level 2 to call diffusion_zt_lhs with level 1 as the input level (the last
1657 : ! variable being passed in during the function call). However, the other
1658 : ! variables passed in (rho_ds_zm*K_zm, gr%dzt(1,:), and gr%dzm(1,:) variables) will
1659 : ! have to be passed in as solving for level 2.
1660 : !
1661 : ! The value of xm(1) is located below the model surface and does not effect
1662 : ! the rest of the model. Since xm can be either a horizontal wind component
1663 : ! or a generic eddy scalar quantity, the value of xm(1) is simply set to the
1664 : ! value of xm(2) after the equation matrix has been solved.
1665 : !
1666 : !
1667 : ! Conservation Properties:
1668 : !
1669 : ! When a fixed-flux lower boundary condition is used (combined with a
1670 : ! zero-flux upper boundary condition), this technique of discretizing the
1671 : ! turbulent advection term (treated as an eddy-diffusion term) leads to
1672 : ! conservative differencing. When the implicit momentum surface flux is
1673 : ! either zero or not used, the column totals for each column in the
1674 : ! left-hand side matrix (for the turbulent advection term) should be equal
1675 : ! to 0. Otherwise, the column total for the second column will be equal to
1676 : ! rho_ds_zm(1) * x'w'|_sfc<t+1>. When the generalized explicit surface
1677 : ! flux is either zero or not used, the column total for the right-hand side
1678 : ! vector (for the turbulent advection term) should be equal to 0.
1679 : ! Otherwise, the column total for the right-hand side vector (for the
1680 : ! turbulent advection term) will be equal to rho_ds_zm(1) * x'w'|_sfc<t>.
1681 : ! This ensures that the total amount of quantity xm over the entire vertical
1682 : ! domain is only changed by the surface flux (neglecting any forcing terms).
1683 : ! The total amount of change is equal to rho_ds_zm(1) * x'w'|_sfc.
1684 : !
1685 : ! To see that this conservation law is satisfied by the left-hand side
1686 : ! matrix, compute the turbulent advection (treated as eddy diffusion) of xm,
1687 : ! neglecting any implicit momentum surface flux, multiply by rho_ds_zt, and
1688 : ! integrate vertically. In discretized matrix notation (where "i" stands
1689 : ! for the matrix column and "j" stands for the matrix row):
1690 : !
1691 : ! 0 = Sum_j Sum_i
1692 : ! (rho_ds_zt)_i ( 1/dzt )_i
1693 : ! ( 0.5_core_rknd * (1/rho_ds_zt) * dzt * (rho_ds_zm*K_zm*dzm) )_ij (xm<t+1>)_j.
1694 : !
1695 : ! The left-hand side matrix,
1696 : ! ( 0.5_core_rknd * (1/rho_ds_zt) * dzt * (rho_ds_zm*K_zm*dzm) )_ij, is partially
1697 : ! written below. The sum over i in the above equation removes (1/rho_ds_zt)
1698 : ! and dzt everywhere from the matrix below. The sum over j leaves the
1699 : ! column totals that are desired, which are 0.
1700 : !
1701 : ! Left-hand side matrix contributions from the turbulent advection term
1702 : ! (treated as an eddy-diffusion term using a Crank-Nicholson timestep);
1703 : ! first five vertical levels:
1704 : !
1705 : ! ------------------------------------------------------------------------------->
1706 : !k=1 | 0 0 0 0
1707 : ! |
1708 : !k=2 | 0 +0.5* -0.5* 0
1709 : ! | (1/rho_ds_zt(k))* (1/rho_ds_zt(k))*
1710 : ! | dzt(k)* dzt(k)*
1711 : ! | rho_ds_zm(k)* rho_ds_zm(k)*
1712 : ! | K_zm(k)*dzm(k) K_zm(k)*dzm(k)
1713 : ! |
1714 : !k=3 | 0 -0.5* +0.5* -0.5*
1715 : ! | (1/rho_ds_zt(k))* (1/rho_ds_zt(k))* (1/rho_ds_zt(k))*
1716 : ! | dzt(k)* dzt(k)* dzt(k)*
1717 : ! | rho_ds_zm(k-1)* [ rho_ds_zm(k)* rho_ds_zm(k)*
1718 : ! | K_zm(k-1)*dzm(k-1) K_zm(k)*dzm(k) K_zm(k)*dzm(k)
1719 : ! | +rho_ds_zm(k-1)*
1720 : ! | K_zm(k-1)*dzm(k-1) ]
1721 : ! |
1722 : !k=4 | 0 0 -0.5* +0.5*
1723 : ! | (1/rho_ds_zt(k))* (1/rho_ds_zt(k))*
1724 : ! | dzt(k)* dzt(k)*
1725 : ! | rho_ds_zm(k-1)* [ rho_ds_zm(k)*
1726 : ! | K_zm(k-1)*dzm(k-1) K_zm(k)*dzm(k)
1727 : ! | +rho_ds_zm(k-1)*
1728 : ! | K_zm(k-1)*dzm(k-1) ]
1729 : ! |
1730 : !k=5 | 0 0 0 -0.5*
1731 : ! | (1/rho_ds_zt(k))*
1732 : ! | dzt(k)*
1733 : ! | rho_ds_zm(k-1)*
1734 : ! | K_zm(k-1)*dzm(k-1)
1735 : ! \ /
1736 : !
1737 : ! Note: The superdiagonal term from level 4 and both the main diagonal and
1738 : ! superdiagonal terms from level 5 are not shown on this diagram.
1739 : !
1740 : ! Note: If an implicit momentum surface flux is used, an additional term,
1741 : ! + (1/rho_ds_zt(2)) * dzt(2) * rho_ds_zm(1)
1742 : ! * [ u_star^2 / sqrt( um(2,<t>)^2 + vm(2,<t>)^2 ) ], is added to
1743 : ! row 2 (k=2), column 2.
1744 : !
1745 : ! To see that the above conservation law is satisfied by the right-hand side
1746 : ! vector, compute the turbulent advection (treated as eddy diffusion) of xm,
1747 : ! neglecting any generalized explicit surface flux, multiply by rho_ds_zt,
1748 : ! and integrate vertically. In discretized matrix notation (where "i"
1749 : ! stands for the matrix column and "j" stands for the matrix row):
1750 : !
1751 : ! 0 = Sum_j Sum_i (rho_ds_zt)_i ( 1/dzt )_i ( rhs_vector )_j.
1752 : !
1753 : ! The right-hand side vector, ( rhs_vector )_j, is partially written below.
1754 : ! The sum over i in the above equation removes (1/rho_ds_zt) and dzt
1755 : ! everywhere from the vector below. The sum over j leaves the column total
1756 : ! that is desired, which is 0.
1757 : !
1758 : ! Right-hand side vector contributions from the turbulent advection term
1759 : ! (treated as an eddy-diffusion term using a Crank-Nicholson timestep);
1760 : ! first five vertical levels:
1761 : !
1762 : ! --------------------------------------------
1763 : !k=1 | 0 |
1764 : ! | |
1765 : ! | |
1766 : !k=2 | +0.5*(1/rho_ds_zt(k))* |
1767 : ! | dzt(k)* |
1768 : ! | [ rho_ds_zm(k)*K_zm(k)* |
1769 : ! | dzm(k)*(xm(k+1,<t>)-xm(k,<t>)) ] |
1770 : ! | |
1771 : !k=3 | +0.5*(1/rho_ds_zt(k))* |
1772 : ! | dzt(k)* |
1773 : ! | [ rho_ds_zm(k)*K_zm(k)* |
1774 : ! | dzm(k)*(xm(k+1,<t>)-xm(k,<t>)) |
1775 : ! | -rho_ds_zm(k-1)*K_zm(k-1)* |
1776 : ! | dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
1777 : ! | |
1778 : !k=4 | +0.5*(1/rho_ds_zt(k))* |
1779 : ! | dzt(k)* |
1780 : ! | [ rho_ds_zm(k)*K_zm(k)* |
1781 : ! | dzm(k)*(xm(k+1,<t>)-xm(k,<t>)) |
1782 : ! | -rho_ds_zm(k-1)*K_zm(k-1)* |
1783 : ! | dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
1784 : ! | |
1785 : !k=5 | +0.5*(1/rho_ds_zt(k))* |
1786 : ! | dzt(k)* |
1787 : ! | [ rho_ds_zm(k)*K_zm(k)* |
1788 : ! | dzm(k)*(xm(k+1,<t>)-xm(k,<t>)) |
1789 : ! | -rho_ds_zm(k-1)*K_zm(k-1)* |
1790 : ! | dzm(k-1)*(xm(k,<t>)-xm(k-1,<t>)) ] |
1791 : ! \ / \ /
1792 : !
1793 : ! Note: If a generalized explicit surface flux is used, an additional term,
1794 : ! + (1/rho_ds_zt(2)) * dzt(2) * rho_ds_zm(1) * x'w'|_sfc, is added to
1795 : ! row 2 (k=2).
1796 : !
1797 : ! Note: Only the contributions by the turbulent advection term are shown
1798 : ! for both the left-hand side matrix and the right-hand side vector.
1799 : ! There are more terms in the equation, and thus more factors to be
1800 : ! added to both the left-hand side matrix (such as time tendency and
1801 : ! mean advection) and the right-hand side vector (such as xm
1802 : ! forcings). The left-hand side matrix is set-up so that a singular
1803 : ! matrix is not encountered.
1804 :
1805 : ! References:
1806 : ! Eqn. 8 & 9 on p. 3545 of
1807 : ! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
1808 : ! Method and Model Description'' Golaz, et al. (2002)
1809 : ! JAS, Vol. 59, pp. 3540--3551.
1810 : !-----------------------------------------------------------------------
1811 :
1812 : use grid_class, only: &
1813 : grid ! Type
1814 :
1815 : use matrix_solver_wrapper, only: &
1816 : tridiag_solve ! Procedure(s)
1817 :
1818 : use stats_variables, only: &
1819 : stats_metadata_type
1820 :
1821 : use stats_type_utilities, only: &
1822 : stat_update_var_pt ! Subroutine
1823 :
1824 : use clubb_precision, only: &
1825 : core_rknd ! Variable(s)
1826 :
1827 : use stats_type, only: stats ! Type
1828 :
1829 : implicit none
1830 :
1831 : ! Constant parameters
1832 :
1833 : integer, parameter :: &
1834 : kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
1835 : k_tdiag = 2, & ! Thermodynamic main diagonal index.
1836 : km1_tdiag = 3 ! Thermodynamic subdiagonal index.
1837 :
1838 : ! ------------------------ Input Variables ------------------------
1839 : integer, intent(in) :: &
1840 : nz, &
1841 : ngrdcol
1842 :
1843 : type (grid), target, intent(in) :: &
1844 : gr
1845 :
1846 : integer, intent(in) :: &
1847 : nrhs ! Number of right-hand side (explicit) vectors & Number of solution vectors.
1848 :
1849 : integer, intent(in) :: &
1850 : ixm_matrix_condt_num ! Stats index of the condition numbers
1851 :
1852 : integer, intent(in) :: &
1853 : tridiag_solve_method ! Specifier for method to solve tridiagonal systems
1854 :
1855 : type (stats_metadata_type), intent(in) :: &
1856 : stats_metadata
1857 :
1858 : ! ------------------------ Inout variables ------------------------
1859 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
1860 : stats_sfc
1861 :
1862 : real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(inout) :: &
1863 : lhs ! Implicit contributions to um, vm, and eddy scalars [units vary]
1864 :
1865 : real( kind = core_rknd ), dimension(ngrdcol,nz,nrhs), intent(inout) :: &
1866 : rhs ! Right-hand side (explicit) contributions.
1867 :
1868 : ! ------------------------ Output variables ------------------------
1869 : real( kind = core_rknd ), dimension(ngrdcol,nz,nrhs), intent(out) :: &
1870 : solution ! Solution to the system of equations [units vary]
1871 :
1872 : ! ------------------------ Local variables ------------------------
1873 : real( kind = core_rknd ), dimension(ngrdcol) :: &
1874 17870112 : rcond ! Estimate of the reciprocal of the condition number on the LHS matrix
1875 :
1876 : integer :: i
1877 :
1878 : ! ------------------------ Begin Code ------------------------
1879 :
1880 : ! Solve tridiagonal system for xm.
1881 8935056 : if ( stats_metadata%l_stats_samp .and. ixm_matrix_condt_num > 0 ) then
1882 :
1883 : call tridiag_solve( "windm_edsclrm", tridiag_solve_method, & ! Intent(in)
1884 : ngrdcol, nz, nrhs, & ! Intent(in)
1885 : lhs, rhs, & ! Intent(inout)
1886 0 : solution, rcond ) ! Intent(out)
1887 :
1888 : ! Est. of the condition number of the variance LHS matrix
1889 0 : do i = 1, ngrdcol
1890 0 : call stat_update_var_pt( ixm_matrix_condt_num, 1, 1.0_core_rknd/rcond(i), & ! intent(in)
1891 0 : stats_sfc(i) ) ! intent(inout)
1892 : end do
1893 : else
1894 :
1895 : call tridiag_solve( "windm_edsclrm", tridiag_solve_method, & ! Intent(in)
1896 : ngrdcol, nz, nrhs, & ! Intent(in)
1897 : lhs, rhs, & ! Intent(inout)
1898 8935056 : solution ) ! Intent(out)
1899 : end if
1900 :
1901 8935056 : return
1902 : end subroutine windm_edsclrm_solve
1903 :
1904 : !=============================================================================
1905 0 : subroutine windm_edsclrm_implicit_stats( nz, solve_type, xm, &
1906 0 : invrs_dzt, &
1907 0 : lhs_diff, lhs_ma_zt, &
1908 0 : invrs_rho_ds_zt, u_star_sqd,&
1909 0 : rho_ds_zm, wind_speed, &
1910 : l_imp_sfc_momentum_flux, &
1911 : stats_metadata, &
1912 : stats_zt )
1913 :
1914 : ! Description:
1915 : ! Compute implicit contributions to um and vm
1916 :
1917 : ! References:
1918 : ! None
1919 : !-----------------------------------------------------------------------
1920 :
1921 : use constants_clubb, only: &
1922 : zero
1923 :
1924 : use stats_type_utilities, only: &
1925 : stat_end_update_pt, & ! Subroutines
1926 : stat_update_var_pt
1927 :
1928 : use clubb_precision, only: &
1929 : core_rknd
1930 :
1931 : use stats_type, only: &
1932 : stats ! Type
1933 :
1934 : use stats_variables, only: &
1935 : stats_metadata_type
1936 :
1937 : implicit none
1938 :
1939 : !---------------------- Input Variables ----------------------
1940 : integer, intent(in) :: &
1941 : nz
1942 :
1943 : integer, intent(in) :: &
1944 : solve_type ! Desc. of what is being solved for
1945 :
1946 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
1947 : xm, & ! Computed value um or vm at <t+1> [m/s]
1948 : invrs_dzt ! The inverse spacing between momentum grid levels;
1949 : ! centered over thermodynamic grid levels.
1950 :
1951 : real( kind = core_rknd ), dimension(3,nz), intent(in) :: &
1952 : lhs_diff, & ! LHS diffustion terms
1953 : lhs_ma_zt ! LHS mean advection terms
1954 :
1955 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
1956 : wind_speed, & ! wind speed; sqrt(u^2 + v^2) [m/s]
1957 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
1958 : invrs_rho_ds_zt ! Inv. dry, static density at thermo. levels [m^3/kg]
1959 :
1960 : real( kind = core_rknd ), intent(in) :: &
1961 : u_star_sqd ! Surface friction velocity, u_star, squared [m/s]
1962 :
1963 : logical, intent(in) :: &
1964 : l_imp_sfc_momentum_flux ! Flag for implicit momentum surface fluxes.
1965 :
1966 : type (stats_metadata_type), intent(in) :: &
1967 : stats_metadata
1968 :
1969 : !---------------------- InOut variables ----------------------
1970 : type (stats), target, intent(inout) :: &
1971 : stats_zt
1972 :
1973 : !---------------------- Local variables ----------------------
1974 : integer :: k, kp1, km1 ! Array indices
1975 :
1976 : real( kind = core_rknd ), dimension(nz) :: &
1977 0 : imp_sfc_flux
1978 :
1979 : ! Budget indices
1980 : integer :: ixm_ma, ixm_ta
1981 :
1982 : !---------------------- Begin Code ----------------------
1983 :
1984 0 : select case ( solve_type )
1985 : case ( windm_edsclrm_um )
1986 0 : ixm_ma = stats_metadata%ium_ma
1987 0 : ixm_ta = stats_metadata%ium_ta
1988 :
1989 : case ( windm_edsclrm_vm )
1990 0 : ixm_ma = stats_metadata%ivm_ma
1991 0 : ixm_ta = stats_metadata%ivm_ta
1992 :
1993 : case default
1994 0 : ixm_ma = 0
1995 0 : ixm_ta = 0
1996 :
1997 : end select
1998 :
1999 0 : imp_sfc_flux(:) = zero
2000 :
2001 0 : if ( l_imp_sfc_momentum_flux ) then
2002 :
2003 : ! Statistics: implicit contributions for um or vm.
2004 :
2005 : ! xm term ta is modified at level 2 to include the effects of the
2006 : ! surface flux. In this case, this effects the implicit portion of
2007 : ! the term, which handles the main diagonal for the turbulent advection term
2008 0 : if ( stats_metadata%ium_ta + stats_metadata%ivm_ta > 0 ) then
2009 0 : imp_sfc_flux(2) = - invrs_rho_ds_zt(2) * invrs_dzt(2) &
2010 0 : * rho_ds_zm(1) * ( u_star_sqd / wind_speed(2) )
2011 : endif
2012 :
2013 : endif ! l_imp_sfc_momentum_flux
2014 :
2015 : ! Finalize implicit contributions for xm
2016 :
2017 0 : do k = 2, nz-1, 1
2018 :
2019 0 : km1 = max( k-1, 2 )
2020 0 : kp1 = min( k+1, nz )
2021 :
2022 : ! xm mean advection
2023 : ! xm term ma is completely implicit; call stat_update_var_pt.
2024 : call stat_update_var_pt( ixm_ma, k, & ! intent(in)
2025 0 : - lhs_ma_zt(3,k) * xm(km1) &
2026 : - lhs_ma_zt(2,k) * xm(k) &
2027 0 : - lhs_ma_zt(1,k) * xm(kp1), & ! intent(in)
2028 0 : stats_zt ) ! intent(inout)
2029 :
2030 : ! xm turbulent transport (implicit component)
2031 : ! xm term ta has both implicit and explicit components;
2032 : ! call stat_end_update_pt.
2033 : call stat_end_update_pt( ixm_ta, k, & ! intent(in)
2034 0 : - 0.5_core_rknd * lhs_diff(3,k) * xm(km1) &
2035 : + ( - 0.5_core_rknd * lhs_diff(2,k) &
2036 : + imp_sfc_flux(k) ) * xm(k) &
2037 : - 0.5_core_rknd * lhs_diff(1,k) * xm(kp1), & ! intent(in)
2038 0 : stats_zt ) ! intent(inout)
2039 :
2040 : enddo
2041 :
2042 :
2043 : ! Upper boundary conditions
2044 0 : k = nz
2045 0 : km1 = max( k-1, 1 )
2046 :
2047 : ! xm mean advection
2048 : ! xm term ma is completely implicit; call stat_update_var_pt.
2049 : call stat_update_var_pt( ixm_ma, k, & ! intent(in)
2050 0 : - lhs_ma_zt(3,k) * xm(km1) &
2051 : - lhs_ma_zt(2,k) * xm(k), & ! intent(in)
2052 0 : stats_zt ) ! intent(inout)
2053 :
2054 : ! xm turbulent transport (implicit component)
2055 : ! xm term ta has both implicit and explicit components;
2056 : ! call stat_end_update_pt.
2057 : call stat_end_update_pt( ixm_ta, k, & ! intent(in)
2058 0 : - 0.5_core_rknd * lhs_diff(3,k) * xm(km1) &
2059 : + ( - 0.5_core_rknd * lhs_diff(2,k) &
2060 : + imp_sfc_flux(k) ) * xm(k), & ! intent(in)
2061 0 : stats_zt ) ! intent(inout)
2062 :
2063 :
2064 0 : return
2065 : end subroutine windm_edsclrm_implicit_stats
2066 :
2067 : !=============================================================================
2068 0 : subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
2069 0 : fcor, perp_wind_m, perp_wind_g, &
2070 0 : xm_forcing, l_implemented, &
2071 : stats_metadata, &
2072 0 : stats_zt, &
2073 0 : xm_tndcy )
2074 :
2075 : ! Description:
2076 : ! Computes the explicit tendency for the um and vm wind components.
2077 : !
2078 : ! The only explicit tendency that is involved in the d(um)/dt or d(vm)/dt
2079 : ! equations is the Coriolis tendency.
2080 : !
2081 : ! The d(um)/dt equation contains the term:
2082 : !
2083 : ! - f * ( v_g - vm );
2084 : !
2085 : ! where f is the Coriolis parameter and v_g is the v component of the
2086 : ! geostrophic wind.
2087 : !
2088 : ! Likewise, the d(vm)/dt equation contains the term:
2089 : !
2090 : ! + f * ( u_g - um );
2091 : !
2092 : ! where u_g is the u component of the geostrophic wind.
2093 : !
2094 : ! This term is treated completely explicitly. The values of um, vm, u_g,
2095 : ! and v_g are all found on the thermodynamic levels.
2096 : !
2097 : ! Wind forcing from the GCSS cases is also added here.
2098 : !
2099 : ! References:
2100 : !-----------------------------------------------------------------------
2101 :
2102 : use grid_class, only: &
2103 : grid ! Type
2104 :
2105 : use stats_type_utilities, only: &
2106 : stat_update_var
2107 :
2108 : use stats_variables, only: &
2109 : stats_metadata_type
2110 :
2111 : use clubb_precision, only: &
2112 : core_rknd ! Variable(s)
2113 :
2114 : use stats_type, only: stats ! Type
2115 :
2116 : implicit none
2117 :
2118 : ! -------------------------- Input Variables --------------------------
2119 : integer, intent(in) :: &
2120 : nz, &
2121 : ngrdcol
2122 :
2123 : integer, intent(in) :: &
2124 : solve_type ! Description of what is being solved for
2125 :
2126 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
2127 : fcor ! Coriolis parameter [s^-1]
2128 :
2129 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2130 : perp_wind_m, & ! Perpendicular component of the mean wind (e.g. v, for the u-eqn) [m/s]
2131 : perp_wind_g, & ! Perpendicular component of the geostropic wind (e.g. vg) [m/s]
2132 : xm_forcing ! Prescribed wind forcing [m/s/s]
2133 :
2134 : logical, intent(in) :: &
2135 : l_implemented ! Flag for CLUBB being implemented in a larger model.
2136 :
2137 : type (stats_metadata_type), intent(in) :: &
2138 : stats_metadata
2139 :
2140 : ! -------------------------- InOut Variables --------------------------
2141 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
2142 : stats_zt
2143 :
2144 : ! -------------------------- Output Variables --------------------------
2145 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2146 : xm_tndcy ! xm tendency [m/s^2]
2147 :
2148 : ! -------------------------- Local Variables --------------------------
2149 : integer :: &
2150 : ixm_gf, &
2151 : ixm_cf, &
2152 : ixm_f
2153 :
2154 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
2155 0 : xm_gf, &
2156 0 : xm_cf
2157 :
2158 : integer :: i, k
2159 :
2160 : ! -------------------------- Begin Code --------------------------
2161 :
2162 : !$acc enter data create( xm_gf, xm_cf )
2163 :
2164 : ! Initialize variables to 0.
2165 0 : xm_gf = 0.0_core_rknd
2166 0 : xm_cf = 0.0_core_rknd
2167 0 : xm_tndcy = 0.0_core_rknd
2168 :
2169 0 : if ( .not. l_implemented ) then
2170 : ! Only compute the Coriolis term if the model is running on it's own,
2171 : ! and is not part of a larger, host model.
2172 :
2173 0 : select case ( solve_type )
2174 :
2175 : case ( windm_edsclrm_um )
2176 :
2177 0 : ixm_gf = stats_metadata%ium_gf
2178 0 : ixm_cf = stats_metadata%ium_cf
2179 0 : ixm_f = stats_metadata%ium_f
2180 :
2181 : !$acc parallel loop gang vector collapse(2) default(present)
2182 0 : do k = 2, nz
2183 0 : do i = 1, ngrdcol
2184 0 : xm_gf(i,k) = - fcor(i) * perp_wind_g(i,k)
2185 : end do
2186 : end do
2187 : !$acc end parallel loop
2188 :
2189 : !$acc parallel loop gang vector collapse(2) default(present)
2190 0 : do k = 2, nz
2191 0 : do i = 1, ngrdcol
2192 0 : xm_cf(i,k) = fcor(i) * perp_wind_m(i,k)
2193 : end do
2194 : end do
2195 : !$acc end parallel loop
2196 :
2197 : case ( windm_edsclrm_vm )
2198 :
2199 0 : ixm_gf = stats_metadata%ivm_gf
2200 0 : ixm_cf = stats_metadata%ivm_cf
2201 0 : ixm_f = stats_metadata%ivm_f
2202 :
2203 : !$acc parallel loop gang vector collapse(2) default(present)
2204 0 : do k = 2, nz
2205 0 : do i = 1, ngrdcol
2206 0 : xm_gf(i,k) = fcor(i) * perp_wind_g(i,k)
2207 : end do
2208 : end do
2209 : !$acc end parallel loop
2210 :
2211 : !$acc parallel loop gang vector collapse(2) default(present)
2212 0 : do k = 2, nz
2213 0 : do i = 1, ngrdcol
2214 0 : xm_cf(i,k) = -fcor(i) * perp_wind_m(i,k)
2215 : end do
2216 : end do
2217 : !$acc end parallel loop
2218 :
2219 : case default
2220 :
2221 0 : ixm_gf = 0
2222 0 : ixm_cf = 0
2223 0 : ixm_f = 0
2224 :
2225 : !$acc parallel loop gang vector collapse(2) default(present)
2226 0 : do k = 2, nz
2227 0 : do i = 1, ngrdcol
2228 0 : xm_gf(i,k) = 0._core_rknd
2229 0 : xm_cf(i,k) = 0._core_rknd
2230 : end do
2231 : end do
2232 : !$acc end parallel loop
2233 :
2234 : end select
2235 :
2236 : !$acc parallel loop gang vector collapse(2) default(present)
2237 0 : do k = 2, nz
2238 0 : do i = 1, ngrdcol
2239 0 : xm_tndcy(i,k) = xm_gf(i,k) + xm_cf(i,k) + xm_forcing(i,k)
2240 : end do
2241 : end do
2242 : !$acc end parallel loop
2243 :
2244 0 : if ( stats_metadata%l_stats_samp ) then
2245 :
2246 : !$acc update host( xm_gf, xm_cf, xm_forcing )
2247 :
2248 0 : do i = 1, ngrdcol
2249 : ! xm term gf is completely explicit; call stat_update_var.
2250 0 : call stat_update_var( ixm_gf, xm_gf(i,:), & ! intent(in)
2251 0 : stats_zt(i) ) ! intent(inout)
2252 :
2253 : ! xm term cf is completely explicit; call stat_update_var.
2254 : call stat_update_var( ixm_cf, xm_cf(i,:), & ! intent(in)
2255 0 : stats_zt(i) ) ! intent(inout)
2256 :
2257 : ! xm term F
2258 : call stat_update_var( ixm_f, xm_forcing(i,:), & ! intent(in)
2259 0 : stats_zt(i) ) ! intent(inout)
2260 : end do
2261 : endif
2262 :
2263 : else ! implemented in a host model.
2264 :
2265 : !$acc parallel loop gang vector collapse(2) default(present)
2266 0 : do k = 2, nz
2267 0 : do i = 1, ngrdcol
2268 0 : xm_tndcy(i,k) = 0.0_core_rknd
2269 : end do
2270 : end do
2271 : !$acc end parallel loop
2272 :
2273 : endif
2274 :
2275 : !$acc exit data delete( xm_gf, xm_cf )
2276 :
2277 0 : return
2278 :
2279 : end subroutine compute_uv_tndcy
2280 :
2281 : !======================================================================================
2282 8935056 : subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
2283 8935056 : lhs_ma_zt, lhs_diff, &
2284 8935056 : wind_speed, u_star_sqd, &
2285 8935056 : rho_ds_zm, invrs_rho_ds_zt, &
2286 : l_implemented, l_imp_sfc_momentum_flux, &
2287 8935056 : lhs )
2288 : ! Description:
2289 : ! Calculate the implicit portion of the horizontal wind or eddy-scalar
2290 : ! time-tendency equation. See the description in subroutine
2291 : ! windm_edsclrm_solve for more details.
2292 : !
2293 : ! Notes:
2294 : ! Lower Boundary:
2295 : ! The lower boundary condition is a fixed-flux boundary condition, which
2296 : ! gets added into the time-tendency equation at level 2.
2297 : ! The value of xm(1) is located below the model surface and does not effect
2298 : ! the rest of the model. Since xm can be either a horizontal wind component
2299 : ! or a generic eddy scalar quantity, the value of xm(1) is simply set to the
2300 : ! value of xm(2) after the equation matrix has been solved.
2301 : !
2302 : ! --- THIS SUBROUTINE HAS BEEN OPTIMIZED ---
2303 : ! Simple changes to this procedure may adversely affect computational speed
2304 : ! - Gunther Huebler, Aug. 2018, clubb:ticket:834
2305 : !----------------------------------------------------------------------------------
2306 :
2307 : use grid_class, only: &
2308 : grid ! Type
2309 :
2310 : use clubb_precision, only: &
2311 : core_rknd ! Variable(s)
2312 :
2313 : implicit none
2314 :
2315 : ! ----------------------- Input Variables -----------------------
2316 : integer, intent(in) :: &
2317 : nz, &
2318 : ngrdcol
2319 :
2320 : type (grid), target, intent(in) :: gr
2321 :
2322 : real( kind = core_rknd ), intent(in) :: &
2323 : dt ! Model timestep [s]
2324 :
2325 : real( kind = core_rknd ), intent(in), dimension(3,ngrdcol,nz) :: &
2326 : lhs_diff, & ! LHS diffustion terms
2327 : lhs_ma_zt ! LHS mean advection terms
2328 :
2329 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2330 : wind_speed, & ! wind speed; sqrt( u^2 + v^2 ) [m/s]
2331 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
2332 : invrs_rho_ds_zt ! Inv. dry, static density at thermo. levels [m^3/kg]
2333 :
2334 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
2335 : u_star_sqd ! Surface friction velocity, u_*, squared [m/s]
2336 :
2337 : logical, intent(in) :: &
2338 : l_implemented, & ! Flag for CLUBB being implemented in a larger model.
2339 : l_imp_sfc_momentum_flux ! Flag for implicit momentum surface fluxes.
2340 :
2341 : ! ----------------------- Output Variable -----------------------
2342 : real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
2343 : lhs ! Implicit contributions to xm (tridiagonal matrix)
2344 :
2345 : ! ----------------------- Local Variables -----------------------
2346 : real( kind = core_rknd ) :: &
2347 : invrs_dt ! Inverse of dt, 1/dt, used for computational efficiency
2348 :
2349 : integer :: k, i ! Loop variable
2350 :
2351 : ! ----------------------- Begin Code -----------------------
2352 :
2353 : ! Calculate coefs of eddy diffusivity and inverse of dt
2354 8935056 : invrs_dt = 1.0_core_rknd / dt
2355 :
2356 : ! Set lower boundary, see notes
2357 : !$acc parallel loop gang vector default(present)
2358 149194656 : do i = 1, ngrdcol
2359 140259600 : lhs(1,i,1) = 0.0_core_rknd
2360 140259600 : lhs(2,i,1) = 1.0_core_rknd
2361 149194656 : lhs(3,i,1) = 0.0_core_rknd
2362 : end do
2363 : !$acc end parallel loop
2364 :
2365 : ! Add terms to lhs
2366 : !$acc parallel loop gang vector collapse(2) default(present)
2367 759479760 : do k = 2, nz
2368 12541286160 : do i = 1, ngrdcol
2369 11781806400 : lhs(1,i,k) = 0.5_core_rknd * lhs_diff(1,i,k)
2370 :
2371 11781806400 : lhs(2,i,k) = 0.5_core_rknd * lhs_diff(2,i,k)
2372 :
2373 : ! LHS time tendency.
2374 11781806400 : lhs(2,i,k) = lhs(2,i,k) + invrs_dt
2375 :
2376 12532351104 : lhs(3,i,k) = 0.5_core_rknd * lhs_diff(3,i,k)
2377 : end do
2378 : end do
2379 : !$acc end parallel loop
2380 :
2381 : ! LHS mean advection term.
2382 8935056 : if ( .not. l_implemented ) then
2383 : !$acc parallel loop gang vector collapse(2) default(present)
2384 0 : do k = 2, nz-1
2385 0 : do i = 1, ngrdcol
2386 0 : lhs(1:3,i,k) = lhs(1:3,i,k) + lhs_ma_zt(:,i,k)
2387 : end do
2388 : end do
2389 : !$acc end parallel loop
2390 :
2391 : endif
2392 :
2393 8935056 : if ( l_imp_sfc_momentum_flux ) then
2394 :
2395 : ! LHS momentum surface flux.
2396 : !$acc parallel loop gang vector default(present)
2397 0 : do i = 1, ngrdcol
2398 0 : lhs(2,i,2) = lhs(2,i,2) + invrs_rho_ds_zt(i,2) * gr%invrs_dzt(i,2) &
2399 0 : * rho_ds_zm(i,1) * ( u_star_sqd(i) / wind_speed(i,2) )
2400 : end do
2401 : !$acc end parallel loop
2402 : end if ! l_imp_sfc_momentum_flux
2403 :
2404 8935056 : return
2405 :
2406 : end subroutine windm_edsclrm_lhs
2407 :
2408 : !=============================================================================
2409 116155728 : subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
2410 116155728 : lhs_diff, xm, xm_tndcy, &
2411 116155728 : rho_ds_zm, invrs_rho_ds_zt, &
2412 116155728 : l_imp_sfc_momentum_flux, xpwp_sfc, &
2413 : stats_metadata, &
2414 116155728 : stats_zt, &
2415 116155728 : rhs )
2416 : ! Description:
2417 : ! Calculate the explicit portion of the horizontal wind or eddy-scalar
2418 : ! time-tendency equation. See the description in subroutine
2419 : ! windm_edsclrm_solve for more details.
2420 : !
2421 : ! References:
2422 : ! None
2423 : !
2424 : ! Notes:
2425 : ! The lower boundary condition needs to be applied here at level 2.
2426 : ! The lower boundary condition is a "fixed flux" boundary condition.
2427 : ! The coding is the same as for a zero-flux boundary condition, but with
2428 : ! an extra term added on the right-hand side at the boundary level. For
2429 : ! the rest of the model code, a zero-flux boundary condition is applied
2430 : ! at level 1, and thus subroutine diffusion_zt_lhs is set-up to do that.
2431 : ! In order to apply the same boundary condition code here at level 2, an
2432 : ! adjuster needs to be used to tell diffusion_zt_lhs to use the code at
2433 : ! level 2 that it normally uses at level 1.
2434 : !
2435 : ! --- THIS SUBROUTINE HAS BEEN OPTIMIZED ---
2436 : ! Simple changes to this procedure may adversely affect computational speed
2437 : ! - Gunther Huebler, Aug. 2018, clubb:ticket:834
2438 : !----------------------------------------------------------------------------------
2439 :
2440 : use clubb_precision, only: &
2441 : core_rknd ! Variable(s)
2442 :
2443 : use diffusion, only: &
2444 : diffusion_zt_lhs ! Procedure(s)
2445 :
2446 : use stats_variables, only: &
2447 : stats_metadata_type
2448 :
2449 : use stats_type_utilities, only: &
2450 : stat_begin_update_pt, & ! Procedure(s)
2451 : stat_modify_pt
2452 :
2453 : use grid_class, only: &
2454 : grid ! Type
2455 :
2456 : use stats_type, only: stats ! Type
2457 :
2458 : implicit none
2459 :
2460 : !------------------- Input Variables -------------------
2461 : integer, intent(in) :: &
2462 : nz, &
2463 : ngrdcol
2464 :
2465 : type (grid), target, intent(in) :: &
2466 : gr
2467 :
2468 : integer, intent(in) :: &
2469 : solve_type ! Description of what is being solved for
2470 :
2471 : real( kind = core_rknd ), intent(in) :: &
2472 : dt ! Model timestep [s]
2473 :
2474 : real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(in) :: &
2475 : lhs_diff ! LHS diffustion terms
2476 :
2477 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
2478 : xm, & ! Eddy-scalar variable, xm (thermo. levels) [units vary]
2479 : xm_tndcy, & ! The explicit time-tendency acting on xm [units vary]
2480 : rho_ds_zm, & ! Dry, static density on momentum levels [kg/m^3]
2481 : invrs_rho_ds_zt ! Inv. dry, static density at thermo. levels [m^3/kg]
2482 :
2483 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
2484 : xpwp_sfc ! x'w' at the surface [units vary]
2485 :
2486 : logical, intent(in) :: &
2487 : l_imp_sfc_momentum_flux ! Flag for implicit momentum surface fluxes.
2488 :
2489 : type (stats_metadata_type), intent(in) :: &
2490 : stats_metadata
2491 :
2492 : !------------------- Inout Variable -------------------
2493 : type (stats), target, dimension(ngrdcol), intent(inout) :: &
2494 : stats_zt
2495 :
2496 : !------------------- Output Variable -------------------
2497 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
2498 : rhs ! Right-hand side (explicit) contributions.
2499 :
2500 : !------------------- Local Variables -------------------
2501 : integer :: i, k ! Loop variable
2502 :
2503 : real( kind = core_rknd ) :: invrs_dt
2504 :
2505 : integer :: ixm_ta
2506 :
2507 : !------------------- Begin Code -------------------
2508 :
2509 116155728 : invrs_dt = 1.0_core_rknd / dt ! Precalculate 1.0/dt to avoid redoing the divide
2510 :
2511 116155728 : select case ( solve_type )
2512 : case ( windm_edsclrm_um )
2513 0 : ixm_ta = stats_metadata%ium_ta
2514 : case ( windm_edsclrm_vm )
2515 0 : ixm_ta = stats_metadata%ivm_ta
2516 : case default ! Eddy scalars
2517 116155728 : ixm_ta = 0
2518 : end select
2519 :
2520 : ! For purposes of the matrix equation, rhs(1) is simply set to 0.
2521 : !$acc parallel loop gang vector default(present)
2522 1939530528 : do i = 1, ngrdcol
2523 1939530528 : rhs(i,1) = 0.0_core_rknd
2524 : end do
2525 : !$acc end parallel loop
2526 :
2527 : ! Lower boundary calculation
2528 : !$acc parallel loop gang vector default(present)
2529 1939530528 : do i = 1, ngrdcol
2530 3646749600 : rhs(i,2) = 0.5_core_rknd &
2531 : * ( - lhs_diff(2,i,2) * xm(i,2) &
2532 1823374800 : - lhs_diff(1,i,2) * xm(i,3) ) &
2533 : + xm_tndcy(i,2) & ! RHS forcings
2534 3762905328 : + invrs_dt * xm(i,2) ! RHS time tendency
2535 : end do
2536 : !$acc end parallel loop
2537 :
2538 : ! Non-boundary rhs calculation, this is a highly vectorized loop
2539 : !$acc parallel loop gang vector collapse(2) default(present)
2540 9640925424 : do k = 3, nz-1
2541 >15915*10^7 : do i = 1, ngrdcol
2542 >29903*10^7 : rhs(i,k) = 0.5_core_rknd &
2543 >14951*10^7 : * ( - lhs_diff(3,i,k) * xm(i,k-1) &
2544 : - lhs_diff(2,i,k) * xm(i,k) &
2545 >14951*10^7 : - lhs_diff(1,i,k) * xm(i,k+1) ) &
2546 : + xm_tndcy(i,k) & ! RHS forcings
2547 >75710*10^7 : + invrs_dt * xm(i,k) ! RHS time tendency
2548 : end do
2549 : end do
2550 : !$acc end parallel loop
2551 :
2552 : ! Upper boundary calculation
2553 : !$acc parallel loop gang vector default(present)
2554 1939530528 : do i = 1, ngrdcol
2555 3646749600 : rhs(i,nz) = 0.5_core_rknd &
2556 1823374800 : * ( - lhs_diff(3,i,nz) * xm(i,nz-1) &
2557 : - lhs_diff(2,i,nz) * xm(i,nz) ) &
2558 : + xm_tndcy(i,nz) & ! RHS forcings
2559 5586280128 : + invrs_dt * xm(i,nz) ! RHS time tendency
2560 : end do
2561 : !$acc end parallel loop
2562 :
2563 116155728 : if ( stats_metadata%l_stats_samp .and. ixm_ta > 0 ) then
2564 :
2565 : !$acc update host( lhs_diff, xm )
2566 :
2567 0 : do i = 1, ngrdcol
2568 :
2569 : ! Statistics: explicit contributions for um or vm.
2570 :
2571 : ! xm term ta has both implicit and explicit components; call
2572 : ! stat_begin_update_pt. Since stat_begin_update_pt automatically
2573 : ! subtracts the value sent in, reverse the sign on right-hand side
2574 : ! turbulent advection component.
2575 :
2576 : ! Lower boundary
2577 : call stat_begin_update_pt( ixm_ta, 2, & ! intent(in)
2578 : 0.5_core_rknd &
2579 0 : * ( lhs_diff(2,i,2) * xm(i,2) & ! intent(in)
2580 0 : + lhs_diff(1,i,2) * xm(i,3) ), &
2581 0 : stats_zt(i) ) ! intent(inout)
2582 :
2583 0 : do k = 3, nz-1
2584 :
2585 : call stat_begin_update_pt( ixm_ta, k, & ! intent(in)
2586 : 0.5_core_rknd &
2587 0 : * ( lhs_diff(3,i,k) * xm(i,k-1) &
2588 : + lhs_diff(2,i,k) * xm(i,k) & ! intent(in)
2589 0 : + lhs_diff(1,i,k) * xm(i,k+1) ), &
2590 0 : stats_zt(i) ) ! intent(inout)
2591 : end do
2592 :
2593 : ! Upper boundary
2594 : call stat_begin_update_pt( ixm_ta, nz, &
2595 : 0.5_core_rknd & ! intent(in)
2596 0 : * ( lhs_diff(3,i,nz) * xm(i,nz-1) &
2597 : + lhs_diff(2,i,nz) * xm(i,nz) ), & ! intent(in)
2598 0 : stats_zt(i) ) ! intent(inout)
2599 : end do
2600 : endif
2601 :
2602 116155728 : if ( .not. l_imp_sfc_momentum_flux ) then
2603 :
2604 : ! RHS generalized surface flux.
2605 : !$acc parallel loop gang vector default(present)
2606 1939530528 : do i = 1, ngrdcol
2607 3646749600 : rhs(i,2) = rhs(i,2) + invrs_rho_ds_zt(i,2) &
2608 0 : * gr%invrs_dzt(i,2) &
2609 3762905328 : * rho_ds_zm(i,1) * xpwp_sfc(i)
2610 : end do
2611 : !$acc end parallel loop
2612 :
2613 116155728 : if ( stats_metadata%l_stats_samp .and. ixm_ta > 0 ) then
2614 :
2615 : !$acc update host( invrs_rho_ds_zt, rho_ds_zm, xpwp_sfc )
2616 :
2617 0 : do i = 1, ngrdcol
2618 : call stat_modify_pt( ixm_ta, 2, & ! intent(in)
2619 0 : + invrs_rho_ds_zt(i,2) &
2620 0 : * gr%invrs_dzt(i,2) &
2621 0 : * rho_ds_zm(i,1) * xpwp_sfc(i), & ! intent(in)
2622 0 : stats_zt(i) ) ! intent(inout)
2623 : end do
2624 : end if
2625 :
2626 : endif ! l_imp_sfc_momentum_flux
2627 :
2628 116155728 : return
2629 :
2630 : end subroutine windm_edsclrm_rhs
2631 :
2632 : end module advance_windm_edsclrm_module
|