Line data Source code
1 : !---------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module turbulent_adv_pdf
5 :
6 : ! Description:
7 : ! Calculates the turbulent advection term in the predictive equation for a
8 : ! variance or covariance where turbulent advection is calculated by
9 : ! integrating over the PDF. This includes the following predictive fields:
10 : ! <w'thl'>, <w'rt'>, <rt'^2>, <thl'^2>, and <rt'thl'>, as well as passive
11 : ! scalar fields <w'sclr'>, <sclr'^2>, <sclr'rt'>, and <sclr'thl'>. CLUBB
12 : ! does not produce a PDF for horizontal wind components u and v. However, the
13 : ! horizontal wind variances <u'^2> and <v'^2> still use this code, as well.
14 :
15 : ! References:
16 : !-------------------------------------------------------------------------
17 :
18 : implicit none
19 :
20 : public :: xpyp_term_ta_pdf_lhs, &
21 : xpyp_term_ta_pdf_lhs_godunov, &
22 : xpyp_term_ta_pdf_rhs, &
23 : xpyp_term_ta_pdf_rhs_godunov
24 :
25 : private ! Set default scope
26 :
27 : integer, parameter :: &
28 : ndiags3 = 3
29 :
30 : contains
31 :
32 : !=============================================================================
33 26805168 : subroutine xpyp_term_ta_pdf_lhs( nz, ngrdcol, gr, coef_wpxpyp_implicit, & ! In
34 26805168 : rho_ds_zt, rho_ds_zm, & ! In
35 26805168 : invrs_rho_ds_zm, & ! In
36 : l_upwind_xpyp_turbulent_adv, & ! In
37 26805168 : sgn_turbulent_vel, & ! In
38 26805168 : coef_wpxpyp_implicit_zm, & ! In
39 26805168 : lhs_ta ) ! Out
40 :
41 : ! Description:
42 : ! Turbulent advection of <w'x'>, <x'^2>, and <x'y'>: implicit portion of
43 : ! the code.
44 : !
45 : ! 1) <w'x'>
46 : !
47 : ! The d<w'x'>/dt equation contains a turbulent advection term:
48 : !
49 : ! - (1/rho_ds) * d( rho_ds * <w'^2 x'> )/dz.
50 : !
51 : ! The value of <w'^2 x'> is found by integrating over the multivariate PDF
52 : ! of w and x, as detailed in function calc_wp2xp_pdf, which is found in
53 : ! module pdf_closure_module in pdf_closure_module.F90.
54 : !
55 : ! The equation obtained for <w'^2 x'> is written in terms of PDF parameters.
56 : ! Substitutions that are specific to the type of PDF used are made for the
57 : ! PDF parameters in order to write the <w'x'> turbulent advection term in
58 : ! the following form:
59 : !
60 : ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit.
61 : !
62 : ! For the ADG1 PDF, coef_wp2xp_implicit is a_1 * < w'^3 > / < w'^2 >, where
63 : ! a_1 is 1 / ( 1 - sigma_sqd_w ). The value of term_wp2xp_explicit is 0, as
64 : ! the <w'x'> turbulent advection term is entirely implicit.
65 : !
66 : ! For the new PDF, the calculations of both coef_wp2xp_implicit and
67 : ! term_wp2xp_explicit are detailed in function calc_coefs_wp2xp_semiimpl,
68 : ! which is found in module new_pdf in new_pdf.F90.
69 : !
70 : ! For the new hybrid PDF, the calculation of coef_wp2xp_implicit is
71 : ! detailed in function calc_coef_wp2xp_implicit, which is found in module
72 : ! new_hybrid_pdf in new_hybrid_pdf.F90. The value of term_wp2xp_explicit
73 : ! is 0, as the <w'x'> turbulent advection term is entirely implicit.
74 : !
75 : ! For explicit turbulent advection, the value of coef_wp2xp_implicit is 0
76 : ! and the value of term_wp2xp_explicit is <w'^2 x'>, as calculated by
77 : ! retaining the equation for <w'^2 x'> that is written in terms of PDF
78 : ! parameters. This is a general form that can be used with any type of PDF.
79 : !
80 : ! The <w'x'> turbulent advection term is rewritten as:
81 : !
82 : ! - (1/rho_ds)
83 : ! * d( rho_ds * ( coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit ) )
84 : ! /dz.
85 : !
86 : ! The variable <w'x'> is evaluated at the (t+1) timestep, which allows the
87 : ! <w'x'> turbulent advection term to be expressed semi-implicitly as:
88 : !
89 : ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz
90 : ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
91 : !
92 : ! The implicit portion of <w'x'> turbulent advection term is:
93 : !
94 : ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz.
95 : !
96 : ! Note: When the implicit term is brought over to the left-hand side, the
97 : ! sign is reversed and the leading "-" in front of the implicit
98 : ! d[ ] / dz term is changed to a "+".
99 : !
100 : ! The timestep index (t+1) means that the value of <w'x'> being used is from
101 : ! the next timestep, which is being advanced to in solving the d<w'x'>/dt
102 : ! equation.
103 : !
104 : ! 2) <x'^2>
105 : !
106 : ! The d<x'^2>/dt equation contains a turbulent advection term:
107 : !
108 : ! - (1/rho_ds) * d( rho_ds * <w'x'^2> )/dz;
109 : !
110 : ! The value of <w'x'^2> is found by integrating over the multivariate PDF of
111 : ! w and x, as detailed in function calc_wpxp2_pdf, which is found in module
112 : ! pdf_closure_module in pdf_closure_module.F90.
113 : !
114 : ! The equation obtained for <w'x'^2> is written in terms of PDF parameters.
115 : ! Substitutions that are specific to the type of PDF used are made for the
116 : ! PDF parameters in order to write the <x'^2> turbulent advection term in
117 : ! the following form:
118 : !
119 : ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit.
120 : !
121 : ! For the ADG1 PDF, the value of coef_wpxp2_implicit is
122 : ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >. The value of term_wpxp2_explicit
123 : ! is (1-(1/3)*beta) * a_1^2 * < w'x' >^2 * < w'^3 > / < w'^2 >^2, where
124 : ! a_1 is 1 / ( 1 - sigma_sqd_w ).
125 : !
126 : ! For the new PDF, the calculation of coef_wpxp2_implicit is detailed in
127 : ! function calc_coef_wpxp2_implicit, which is found in module new_pdf in
128 : ! new_pdf.F90. The value of term_wpxp2_explicit is 0, as the <x'^2>
129 : ! turbulent advection term is entirely implicit.
130 : !
131 : ! For the new hybrid PDF, the calculation of both coef_wpxp2_implicit and
132 : ! term_wpxp2_explicit are detailed in subroutine calc_coefs_wpxp2_semiimpl,
133 : ! which is found in module new_hybrid_pdf in new_hybrid_pdf.F90.
134 : !
135 : ! For explicit turbulent advection, the value of coef_wpxp2_implicit is 0
136 : ! and the value of term_wpxp2_explicit is <w'x'^2>, as calculated by
137 : ! retaining the equation for <w'x'^2> that is written in terms of PDF
138 : ! parameters. This is a general form that can be used with any type of PDF.
139 : !
140 : ! The <x'^2> turbulent advection term is rewritten as:
141 : !
142 : ! - (1/rho_ds)
143 : ! * d( rho_ds * ( coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit ) )
144 : ! /dz;
145 : !
146 : ! The variable <x'^2> is evaluated at the (t+1) timestep, which allows the
147 : ! <x'^2> turbulent advection term to be expressed semi-implicitly as:
148 : !
149 : ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz;
150 : ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
151 : !
152 : ! The implicit portion of <x'^2> turbulent advection term is:
153 : !
154 : ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz.
155 : !
156 : ! Note: When the implicit term is brought over to the left-hand side, the
157 : ! sign is reversed and the leading "-" in front of all implicit
158 : ! d[ ] / dz terms is changed to a "+".
159 : !
160 : ! The timestep index (t+1) means that the value of <x'^2> being used is from
161 : ! the next timestep, which is being advanced to in solving the d<x'^2>/dt
162 : ! equation.
163 : !
164 : ! 3) <x'y'>
165 : !
166 : ! The d<x'y'>/dt equation contains a turbulent advection term:
167 : !
168 : ! - (1/rho_ds) * d( rho_ds * <w'x'y'> )/dz.
169 : !
170 : ! The value of <w'x'y'> is found by integrating over the multivariate PDF of
171 : ! w, x, and y, as detailed in function calc_wpxpyp_pdf, which is found in
172 : ! module pdf_closure_module in pdf_closure_module.F90.
173 : !
174 : ! The equation obtained for <w'x'y'> is written in terms of PDF parameters.
175 : ! Substitutions that are specific to the type of PDF used are made for the
176 : ! PDF parameters in order to write the <x'y'> turbulent advection term in
177 : ! the following form:
178 : !
179 : ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit.
180 : !
181 : ! For the ADG1 PDF, the value of coef_wpxpyp_implicit is
182 : ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >. The value of term_wpxpyp_explicit
183 : ! is (1-(1/3)*beta) * a_1^2 * < w'x' > * < w'y' > * < w'^3 > / < w'^2 >^2,
184 : ! where a_1 is 1 / ( 1 - sigma_sqd_w ).
185 : !
186 : ! For the new PDF, the calculation of both coef_wpxpyp_implicit and
187 : ! term_wpxpyp_explicit are detailed in function calc_coefs_wpxpyp_semiimpl,
188 : ! which is found in module new_pdf in new_pdf.F90.
189 : !
190 : ! For the new hybrid PDF, the calculation of both coef_wpxpyp_implicit
191 : ! and term_wpxpyp_explicit are detailed in subroutine
192 : ! calc_coefs_wpxpyp_semiimpl, which is found in module new_hybrid_pdf in
193 : ! new_hybrid_pdf.F90.
194 : !
195 : ! For explicit turbulent advection, the value of coef_wpxpyp_implicit is 0
196 : ! and the value of term_wpxpyp_explicit is <w'x'y'>, as calculated by
197 : ! retaining the equation for <w'x'y'> that is written in terms of PDF
198 : ! parameters. This is a general form that can be used with any type of PDF.
199 : !
200 : ! The <x'y'> turbulent advection term is rewritten as:
201 : !
202 : ! - (1/rho_ds)
203 : ! * d( rho_ds * ( coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit ) )
204 : ! /dz.
205 : !
206 : ! The variable <x'y'> is evaluated at the (t+1) timestep, which allows the
207 : ! <x'y'> turbulent advection term to be expressed semi-implicitly as:
208 : !
209 : ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz
210 : ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
211 : !
212 : ! The implicit portion of <x'y'> turbulent advection term is:
213 : !
214 : ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz.
215 : !
216 : ! Note: When the implicit term is brought over to the left-hand side, the
217 : ! sign is reversed and the leading "-" in front of all implicit
218 : ! d[ ] / dz terms is changed to a "+".
219 : !
220 : ! The timestep index (t+1) means that the value of <x'y'> being used is from
221 : ! the next timestep, which is being advanced to in solving the d<x'y'>/dt
222 : ! equation.
223 : !
224 : ! When x and y are the same variable, <x'y'> reduces to <x'^2> and <w'x'y'>
225 : ! reduces to <w'x'^2>. Likewise, when y is set equal to w, <x'y'> becomes
226 : ! <w'x'> and <w'x'y'> reduces to <w'^2 x'>. The discretization and the code
227 : ! used in this function will be written generally in terms of <x'y'> and
228 : ! coef_wpxpyp_implicit, but also applies to <x'^2> and coef_wpxp2_implicit,
229 : ! as well as to <w'x'> and coef_wp2xp_implicit.
230 : !
231 : ! The implicit discretization of this term is as follows:
232 : !
233 : ! 1) Centered Discretization
234 : !
235 : ! The values of <x'y'> are found on the momentum levels, while the values of
236 : ! coef_wpxpyp_implicit are found on the thermodynamic levels, which is where
237 : ! they were originally calculated by the PDF. Additionally, the values of
238 : ! rho_ds_zt are found on the thermodynamic levels, and the values of
239 : ! invrs_rho_ds_zm are found on the momentum levels. The values of <x'y'>
240 : ! are interpolated to the intermediate thermodynamic levels as <x'y'>|_zt.
241 : ! At the thermodynamic levels, the values of coef_wpxpyp_implicit are
242 : ! multiplied by <x'y'>|_zt, and their product is multiplied by rho_ds_zt.
243 : ! Then, the derivative (d/dz) of that expression is taken over the central
244 : ! momentum level, where it is multiplied by -invrs_rho_ds_zm. This yields
245 : ! the desired result.
246 : !
247 : ! =xpyp============================================================== m(k+1)
248 : !
249 : ! -xpyp_zt(interp)-------coef_wpxpyp_implicit---------rho_ds_zt------ t(k+1)
250 : !
251 : ! =xpyp=d(rho_ds_zt*coef_wpxpyp_implicit*xpyp_zt)/dz=invrs_rho_ds_zm= m(k)
252 : !
253 : ! -xpyp_zt(interp)-------coef_wpxpyp_implicit---------rho_ds_zt------ t(k)
254 : !
255 : ! =xpyp============================================================== m(k-1)
256 : !
257 : ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
258 : ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
259 : ! The letter "t" is used for thermodynamic levels and the letter "m" is used
260 : ! for momentum levels.
261 : !
262 : ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
263 : !
264 : ! 2) "Upwind" Discretization.
265 : !
266 : ! The values of <x'y'> are found on the momentum levels. The values of
267 : ! coef_wpxpyp_implicit are originally calculated by the PDF on the
268 : ! thermodynamic levels. They are interpolated to the intermediate momentum
269 : ! levels as coef_wpxpyp_implicit_zm. Additionally, the values of rho_ds_zm
270 : ! and the values of invrs_rho_ds_zm are found on the momentum levels. The
271 : ! sign of the turbulent velocity is found on the central momentum level. At
272 : ! the momentum levels, the values of coef_wpxpyp_implicit_zm are multiplied
273 : ! by <x'y'>, and their product is multiplied by rho_ds_zm. Then, the
274 : ! derivative (d/dz) of that expression is taken. When the sign of the
275 : ! turbulent velocity is positive, the "wind" is coming from below, and the
276 : ! derivative involves the central momentum level and the momentum level
277 : ! immediately below it. When the sign of the turbulent velocity is
278 : ! negative, the "wind" is coming from above, and the derivative involves the
279 : ! central momentum level and the momenum level immediately above it. After
280 : ! the derivative is taken, it is multiplied by -invrs_rho_ds_zm at the
281 : ! central momentum level. This yields the desired result.
282 : !
283 : ! The turbulent velocity for <x'y'> is <w'x'y'> / <x'y'>, which has units of
284 : ! m/s. The sign of the turbulent velocity is sgn( <w'x'y'> / <x'y'> ),
285 : ! where:
286 : !
287 : ! sgn(x) = | 1; when x >= 0
288 : ! | -1; when x < 0.
289 : !
290 : ! The sign of the turbulent velocity can also be rewritten as
291 : ! sgn( <w'x'y'> ) / sgn( <x'y'> ). When a variance (<x'^2>) is being solved
292 : ! for, y = x, and sgn( <x'^2> ) is always 1. The sign of the turbulent
293 : ! velocity reduces to simply sgn( <w'x'^2> ).
294 : !
295 : ! ---------coef_wpxpyp_implicit-------------------------------------- t(k+2)
296 : !
297 : ! =xpyp=======coef_wpxpyp_implicit_zm(interp.)=======rho_ds_zm======= m(k+1)
298 : !
299 : ! ---------coef_wpxpyp_implicit-------------------------------------- t(k+1)
300 : !
301 : ! =xpyp===coef_wpxpyp_implicit_zm(interp.)=rho_ds_zm=invrs_rho_ds_zm= m(k)
302 : !
303 : ! ---------coef_wpxpyp_implicit-------------------------------------- t(k)
304 : !
305 : ! =xpyp=======coef_wpxpyp_implicit_zm(interp.)=======rho_ds_zm======= m(k-1)
306 : !
307 : ! ---------coef_wpxpyp_implicit-------------------------------------- t(k-1)
308 : !
309 : ! The vertical indices t(k+2), m(k+1), t(k+1), m(k), t(k), m(k-1), and
310 : ! t(k-1) correspond with altitudes zt(k+2), zm(k+1), zt(k+1), zm(k), zt(k),
311 : ! zm(k-1), and zt(k-1), respectively. The letter "t" is used for
312 : ! thermodynamic levels and the letter "m" is used for momentum levels.
313 : !
314 : ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) ); and
315 : ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) ).
316 :
317 : ! References:
318 : !-----------------------------------------------------------------------
319 :
320 : use grid_class, only: & ! gr%weights_zm2zt
321 : grid ! Type
322 :
323 : use constants_clubb, only: &
324 : zero ! Variable(s)
325 :
326 : use clubb_precision, only: &
327 : core_rknd ! Variable(s)
328 :
329 : implicit none
330 :
331 : ! Constant parameters
332 : integer, parameter :: &
333 : kp1_mdiag = 1, & ! Momentum superdiagonal index.
334 : k_mdiag = 2, & ! Momentum main diagonal index.
335 : km1_mdiag = 3 ! Momentum subdiagonal index.
336 :
337 : integer, parameter :: &
338 : m_above = 1, & ! Index for upper momentum level grid weight.
339 : m_below = 2 ! Index for lower momentum level grid weight.
340 :
341 : ! ------------------------------ Input Variables ------------------------------
342 : integer, intent(in) :: &
343 : nz, &
344 : ngrdcol
345 :
346 : type (grid), target, intent(in) :: gr
347 :
348 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
349 : coef_wpxpyp_implicit, & ! Coef. of <x'y'> in <w'x'y'>; t-levs [m/s]
350 : rho_ds_zt, & ! Dry, static density at t-levels [kg/m^3]
351 : invrs_rho_ds_zm ! Inv dry, static density at m-levels [m^3/kg]
352 :
353 : logical, intent(in) :: &
354 : l_upwind_xpyp_turbulent_adv ! Flag to use "upwind" discretization
355 :
356 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
357 : sgn_turbulent_vel, & ! Sign of the turbulent velocity [-]
358 : coef_wpxpyp_implicit_zm, & ! coef_wpxpyp_implicit interp m-levs [m/s]
359 : rho_ds_zm ! Dry, static density at m-levs [kg/m^3]
360 :
361 : ! ------------------------------ Return Variable ------------------------------
362 : real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
363 : lhs_ta ! LHS coefficient of xpyp turbulent advection [1/s]
364 :
365 : ! ------------------------------ Local Variables ------------------------------
366 : integer :: i, k, b ! Vertical level index
367 :
368 : ! ------------------------------ Begin Code ------------------------------
369 :
370 : !$acc data copyin( gr, gr%weights_zm2zt, gr%invrs_dzm, gr%invrs_dzt, &
371 : !$acc coef_wpxpyp_implicit, rho_ds_zt, invrs_rho_ds_zm, &
372 : !$acc sgn_turbulent_vel, coef_wpxpyp_implicit_zm, rho_ds_zm ) &
373 : !$acc copyout( lhs_ta )
374 :
375 :
376 : ! Set lower boundary array to 0
377 : !$acc parallel loop gang vector collapse(2) default(present)
378 447583968 : do i = 1, ngrdcol
379 1709920368 : do b = 1, ndiags3
380 1683115200 : lhs_ta(b,i,1) = zero
381 : end do
382 : end do
383 : !$acc end parallel loop
384 :
385 26805168 : if ( .not. l_upwind_xpyp_turbulent_adv ) then
386 :
387 : ! Centered discretization.
388 : !$acc parallel loop gang vector collapse(2) default(present)
389 750544704 : do k = 2, nz-1, 1
390 12392091504 : do i = 1, ngrdcol
391 :
392 : ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
393 23283093600 : lhs_ta(kp1_mdiag,i,k) &
394 0 : = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
395 11641546800 : * rho_ds_zt(i,k+1) * coef_wpxpyp_implicit(i,k+1) &
396 46566187200 : * gr%weights_zm2zt(i,k+1,m_above)
397 :
398 : ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
399 : lhs_ta(k_mdiag,i,k) &
400 0 : = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
401 : * ( rho_ds_zt(i,k+1) * coef_wpxpyp_implicit(i,k+1) &
402 0 : * gr%weights_zm2zt(i,k+1,m_below) &
403 : - rho_ds_zt(i,k) * coef_wpxpyp_implicit(i,k) &
404 11641546800 : * gr%weights_zm2zt(i,k,m_above) )
405 :
406 : ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
407 : lhs_ta(km1_mdiag,i,k) &
408 0 : = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
409 : * rho_ds_zt(i,k) * coef_wpxpyp_implicit(i,k) &
410 12383156448 : * gr%weights_zm2zt(i,k,m_below)
411 :
412 : end do
413 : end do
414 : !$acc end parallel loop
415 :
416 : else ! l_upwind_xpyp_turbulent_adv
417 :
418 : ! "Upwind" discretization
419 : !$acc parallel loop gang vector collapse(2) default(present)
420 1501089408 : do k = 2, nz-1, 1
421 24784183008 : do i = 1, ngrdcol
422 :
423 24766312896 : if ( sgn_turbulent_vel(i,k) > zero ) then
424 :
425 : ! The "wind" is blowing upward.
426 :
427 : ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
428 4946413726 : lhs_ta(kp1_mdiag,i,k) = zero
429 :
430 : ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
431 : lhs_ta(k_mdiag,i,k) &
432 0 : = invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k) &
433 4946413726 : * rho_ds_zm(i,k) * coef_wpxpyp_implicit_zm(i,k)
434 :
435 : ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
436 : lhs_ta(km1_mdiag,i,k) &
437 0 : = - invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k) &
438 4946413726 : * rho_ds_zm(i,k-1) * coef_wpxpyp_implicit_zm(i,k-1)
439 :
440 : else ! sgn_turbulent_vel < 0
441 :
442 : ! The "wind" is blowing downward.
443 :
444 : ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
445 : lhs_ta(kp1_mdiag,i,k) &
446 0 : = invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k+1) &
447 18336679874 : * rho_ds_zm(i,k+1) * coef_wpxpyp_implicit_zm(i,k+1)
448 :
449 : ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
450 : lhs_ta(k_mdiag,i,k) &
451 0 : = - invrs_rho_ds_zm(i,k) * gr%invrs_dzt(i,k+1) &
452 18336679874 : * rho_ds_zm(i,k) * coef_wpxpyp_implicit_zm(i,k)
453 :
454 : ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
455 18336679874 : lhs_ta(km1_mdiag,i,k) = zero
456 :
457 : end if ! sgn_turbulent_vel
458 :
459 : end do
460 : end do
461 : !$acc end parallel loop
462 :
463 : endif
464 :
465 : ! Set upper boundary array to
466 : !$acc parallel loop gang vector collapse(2) default(present)
467 447583968 : do i = 1, ngrdcol
468 1709920368 : do b = 1, ndiags3
469 1683115200 : lhs_ta(b,i,nz) = zero
470 : end do
471 : end do
472 : !$acc end parallel loop
473 :
474 : !$acc end data
475 :
476 26805168 : return
477 :
478 : end subroutine xpyp_term_ta_pdf_lhs
479 :
480 : !=============================================================================================
481 0 : subroutine xpyp_term_ta_pdf_lhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
482 0 : coef_wpxpyp_implicit, & ! Intent(in)
483 0 : invrs_rho_ds_zm, rho_ds_zm, & ! Intent(in)
484 0 : lhs_ta )
485 : ! Intent(out)
486 : ! Description:
487 : ! This subroutine is a revised version of xpyp_term_ta_pdf_lhs_all. The
488 : ! revisions are maded to use the Godunov-like upwind scheme for the
489 : ! vertical discretization of the turbulent advection term. This subroutine
490 : ! returns an array of 3 dimensional arrays, one for every grid level not
491 : ! including
492 : ! boundary values.
493 : !
494 : ! Optional Arguements:
495 : ! The optional arguements can be used to override the default indices.
496 : ! from_level - low index, default 2
497 : ! to level - high index, default gr%nz-1
498 : !
499 : ! Notes:
500 : ! This subroutine exists for testing of Godunov-like upwind scheme.
501 : ! THIS SUBROUTINE DOES NOT HANDLE BOUNDARY CONDITIONS AND SETS THEM TO 0
502 : !---------------------------------------------------------------------------------------------
503 :
504 : use grid_class, only: & ! for gr%weights_zm2zt
505 : grid ! Type
506 :
507 : use clubb_precision, only: &
508 : core_rknd ! Variable(s)
509 :
510 : implicit none
511 :
512 : integer, parameter :: &
513 : kp1_mdiag = 1, & ! Momentum superdiagonal index.
514 : k_mdiag = 2, & ! Momentum main diagonal index.
515 : km1_mdiag = 3 ! Momentum subdiagonal index.
516 :
517 : !------------------- Input Variables -------------------
518 : integer, intent(in) :: &
519 : nz, &
520 : ngrdcol
521 :
522 : type (grid), target, intent(in) :: gr
523 :
524 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
525 : coef_wpxpyp_implicit, & ! Coef. of <x'y'> in <w'x'y'>; t-lev [m/s]
526 : invrs_rho_ds_zm, & ! Inv dry, static density @ m-level [m^3/kg]
527 : rho_ds_zm ! Dry, static density at m-lev [kg/m^3]
528 :
529 : !------------------- Output Variables -------------------
530 : real( kind = core_rknd ), dimension(ndiags3,ngrdcol,nz), intent(out) :: &
531 : lhs_ta
532 :
533 : !---------------- Local Variables -------------------
534 : integer :: &
535 : i, k, b ! Loop variable for current grid level
536 :
537 : !---------------- Begin Code -------------------
538 :
539 : !$acc data copyin( gr, gr%invrs_dzm, &
540 : !$acc coef_wpxpyp_implicit, invrs_rho_ds_zm, rho_ds_zm ) &
541 : !$acc copyout( lhs_ta )
542 :
543 : ! Set lower boundary array to 0
544 : !$acc parallel loop gang vector collapse(2) default(present)
545 0 : do i = 1, ngrdcol
546 0 : do b = 1, ndiags3
547 0 : lhs_ta(b,i,1) = 0.0_core_rknd
548 : end do
549 : end do
550 : !$acc end parallel loop
551 :
552 : ! Godunov-like upwind discretization
553 : !$acc parallel loop gang vector collapse(2) default(present)
554 0 : do k = 2, nz-1
555 0 : do i = 1, ngrdcol
556 :
557 : ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
558 0 : lhs_ta(kp1_mdiag,i,k) = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
559 0 : * rho_ds_zm(i,k+1) &
560 0 : * min(0.0_core_rknd,coef_wpxpyp_implicit(i,k+1))
561 :
562 : ! Momentum main diagonal: [ x xpyp(k,<t+1>) ]
563 0 : lhs_ta(k_mdiag,i,k) = invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
564 : * rho_ds_zm(i,k) &
565 : * ( max(0.0_core_rknd,coef_wpxpyp_implicit(i,k+1)) - &
566 0 : min(0.0_core_rknd,coef_wpxpyp_implicit(i,k)) )
567 :
568 : ! Momentum subdiagonal: [ x xpyp(k-1,<t+1>) ]
569 0 : lhs_ta(km1_mdiag,i,k) = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
570 0 : * rho_ds_zm(i,k-1) &
571 0 : * max(0.0_core_rknd,coef_wpxpyp_implicit(i,k) )
572 : end do
573 : end do
574 : !$acc end parallel loop
575 :
576 : ! Set upper boundary array to 0
577 : !$acc parallel loop gang vector collapse(2) default(present)
578 0 : do i = 1, ngrdcol
579 0 : do b = 1, ndiags3
580 0 : lhs_ta(b,i,nz) = 0.0_core_rknd
581 : end do
582 : end do
583 : !$acc end parallel loop
584 :
585 : !$acc end data
586 :
587 0 : return
588 :
589 : end subroutine xpyp_term_ta_pdf_lhs_godunov
590 :
591 : !=============================================================================
592 44675280 : subroutine xpyp_term_ta_pdf_rhs( nz, ngrdcol, gr, term_wpxpyp_explicit, & ! In
593 44675280 : rho_ds_zt, rho_ds_zm, & ! In
594 44675280 : invrs_rho_ds_zm, & ! In
595 : l_upwind_xpyp_turbulent_adv, & ! In
596 44675280 : sgn_turbulent_vel, & ! In
597 44675280 : term_wpxpyp_explicit_zm, & ! In
598 44675280 : rhs_ta ) ! Out
599 :
600 : ! Description:
601 : ! Turbulent advection of <w'x'>, <x'^2>, and <x'y'>: explicit portion of
602 : ! the code.
603 : !
604 : ! 1) <w'x'>
605 : !
606 : ! The d<w'x'>/dt equation contains a turbulent advection term:
607 : !
608 : ! - (1/rho_ds) * d( rho_ds * <w'^2 x'> )/dz.
609 : !
610 : ! The value of <w'^2 x'> is found by integrating over the multivariate PDF
611 : ! of w and x, as detailed in function calc_wp2xp_pdf, which is found in
612 : ! module pdf_closure_module in pdf_closure_module.F90.
613 : !
614 : ! The equation obtained for <w'^2 x'> is written in terms of PDF parameters.
615 : ! Substitutions that are specific to the type of PDF used are made for the
616 : ! PDF parameters in order to write the <w'x'> turbulent advection term in
617 : ! the following form:
618 : !
619 : ! <w'^2 x'> = coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit.
620 : !
621 : ! For the ADG1 PDF, coef_wp2xp_implicit is a_1 * < w'^3 > / < w'^2 >, where
622 : ! a_1 is 1 / ( 1 - sigma_sqd_w ). The value of term_wp2xp_explicit is 0, as
623 : ! the <w'x'> turbulent advection term is entirely implicit.
624 : !
625 : ! For the new PDF, the calculations of both coef_wp2xp_implicit and
626 : ! term_wp2xp_explicit are detailed in function calc_coefs_wp2xp_semiimpl,
627 : ! which is found in module new_pdf in new_pdf.F90.
628 : !
629 : ! For the new hybrid PDF, the calculation of coef_wp2xp_implicit is
630 : ! detailed in function calc_coef_wp2xp_implicit, which is found in module
631 : ! new_hybrid_pdf in new_hybrid_pdf.F90. The value of term_wp2xp_explicit
632 : ! is 0, as the <w'x'> turbulent advection term is entirely implicit.
633 : !
634 : ! For explicit turbulent advection, the value of coef_wp2xp_implicit is 0
635 : ! and the value of term_wp2xp_explicit is <w'^2 x'>, as calculated by
636 : ! retaining the equation for <w'^2 x'> that is written in terms of PDF
637 : ! parameters. This is a general form that can be used with any type of PDF.
638 : !
639 : ! The <w'x'> turbulent advection term is rewritten as:
640 : !
641 : ! - (1/rho_ds)
642 : ! * d( rho_ds * ( coef_wp2xp_implicit * <w'x'> + term_wp2xp_explicit ) )
643 : ! /dz.
644 : !
645 : ! The variable <w'x'> is evaluated at the (t+1) timestep, which allows the
646 : ! <w'x'> turbulent advection term to be expressed semi-implicitly as:
647 : !
648 : ! - (1/rho_ds) * d( rho_ds * coef_wp2xp_implicit * <w'x'>(t+1) )/dz
649 : ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
650 : !
651 : ! The explicit portion of <w'x'> turbulent advection term is:
652 : !
653 : ! - (1/rho_ds) * d( rho_ds * term_wp2xp_explicit )/dz.
654 : !
655 : ! 2) <x'^2>
656 : !
657 : ! The d<x'^2>/dt equation contains a turbulent advection term:
658 : !
659 : ! - (1/rho_ds) * d( rho_ds * <w'x'^2> )/dz;
660 : !
661 : ! The value of <w'x'^2> is found by integrating over the multivariate PDF of
662 : ! w and x, as detailed in function calc_wpxp2_pdf, which is found in module
663 : ! pdf_closure_module in pdf_closure_module.F90.
664 : !
665 : ! The equation obtained for <w'x'^2> is written in terms of PDF parameters.
666 : ! Substitutions that are specific to the type of PDF used are made for the
667 : ! PDF parameters in order to write the <x'^2> turbulent advection term in
668 : ! the following form:
669 : !
670 : ! <w'x'^2> = coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit.
671 : !
672 : ! For the ADG1 PDF, the value of coef_wpxp2_implicit is
673 : ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >. The value of term_wpxp2_explicit
674 : ! is (1-(1/3)*beta) * a_1^2 * < w'x' >^2 * < w'^3 > / < w'^2 >^2, where
675 : ! a_1 is 1 / ( 1 - sigma_sqd_w ).
676 : !
677 : ! For the new PDF, the calculation of coef_wpxp2_implicit is detailed in
678 : ! function calc_coef_wpxp2_implicit, which is found in module new_pdf in
679 : ! new_pdf.F90. The value of term_wpxp2_explicit is 0, as the <x'^2>
680 : ! turbulent advection term is entirely implicit.
681 : !
682 : ! For the new hybrid PDF, the calculation of both coef_wpxp2_implicit and
683 : ! term_wpxp2_explicit are detailed in subroutine calc_coefs_wpxp2_semiimpl,
684 : ! which is found in module new_hybrid_pdf in new_hybrid_pdf.F90.
685 : !
686 : ! For explicit turbulent advection, the value of coef_wpxp2_implicit is 0
687 : ! and the value of term_wpxp2_explicit is <w'x'^2>, as calculated by
688 : ! retaining the equation for <w'x'^2> that is written in terms of PDF
689 : ! parameters. This is a general form that can be used with any type of PDF.
690 : !
691 : ! The <x'^2> turbulent advection term is rewritten as:
692 : !
693 : ! - (1/rho_ds)
694 : ! * d( rho_ds * ( coef_wpxp2_implicit * <x'^2> + term_wpxp2_explicit ) )
695 : ! /dz;
696 : !
697 : ! The variable <x'^2> is evaluated at the (t+1) timestep, which allows the
698 : ! <x'^2> turbulent advection term to be expressed semi-implicitly as:
699 : !
700 : ! - (1/rho_ds) * d( rho_ds * coef_wpxp2_implicit * <x'^2>(t+1) )/dz;
701 : ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
702 : !
703 : ! The explicit portion of <x'^2> turbulent advection term is:
704 : !
705 : ! - (1/rho_ds) * d( rho_ds * term_wpxp2_explicit )/dz.
706 : !
707 : ! 3) <x'y'>
708 : !
709 : ! The d<x'y'>/dt equation contains a turbulent advection term:
710 : !
711 : ! - (1/rho_ds) * d( rho_ds * <w'x'y'> )/dz.
712 : !
713 : ! The value of <w'x'y'> is found by integrating over the multivariate PDF of
714 : ! w, x, and y, as detailed in function calc_wpxpyp_pdf, which is found in
715 : ! module pdf_closure_module in pdf_closure_module.F90.
716 : !
717 : ! The equation obtained for <w'x'y'> is written in terms of PDF parameters.
718 : ! Substitutions that are specific to the type of PDF used are made for the
719 : ! PDF parameters in order to write the <x'y'> turbulent advection term in
720 : ! the following form:
721 : !
722 : ! <w'x'y'> = coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit.
723 : !
724 : ! For the ADG1 PDF, the value of coef_wpxpyp_implicit is
725 : ! (1/3)*beta * a_1 * < w'^3 > / < w'^2 >. The value of term_wpxpyp_explicit
726 : ! is (1-(1/3)*beta) * a_1^2 * < w'x' > * < w'y' > * < w'^3 > / < w'^2 >^2,
727 : ! where a_1 is 1 / ( 1 - sigma_sqd_w ).
728 : !
729 : ! For the new PDF, the calculation of both coef_wpxpyp_implicit and
730 : ! term_wpxpyp_explicit are detailed in function calc_coefs_wpxpyp_semiimpl,
731 : ! which is found in module new_pdf in new_pdf.F90.
732 : !
733 : ! For the new hybrid PDF, the calculation of both coef_wpxpyp_implicit
734 : ! and term_wpxpyp_explicit are detailed in subroutine
735 : ! calc_coefs_wpxpyp_semiimpl, which is found in module new_hybrid_pdf in
736 : ! new_hybrid_pdf.F90.
737 : !
738 : ! For explicit turbulent advection, the value of coef_wpxpyp_implicit is 0
739 : ! and the value of term_wpxpyp_explicit is <w'x'y'>, as calculated by
740 : ! retaining the equation for <w'x'y'> that is written in terms of PDF
741 : ! parameters. This is a general form that can be used with any type of PDF.
742 : !
743 : ! The <x'y'> turbulent advection term is rewritten as:
744 : !
745 : ! - (1/rho_ds)
746 : ! * d( rho_ds * ( coef_wpxpyp_implicit * <x'y'> + term_wpxpyp_explicit ) )
747 : ! /dz.
748 : !
749 : ! The variable <x'y'> is evaluated at the (t+1) timestep, which allows the
750 : ! <x'y'> turbulent advection term to be expressed semi-implicitly as:
751 : !
752 : ! - (1/rho_ds) * d( rho_ds * coef_wpxpyp_implicit * <x'y'>(t+1) )/dz
753 : ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
754 : !
755 : ! The explicit portion of <x'y'> turbulent advection term is:
756 : !
757 : ! - (1/rho_ds) * d( rho_ds * term_wpxpyp_explicit )/dz.
758 : !
759 : ! When x and y are the same variable, <x'y'> reduces to <x'^2> and <w'x'y'>
760 : ! reduces to <w'x'^2>. Likewise, when y is set equal to w, <x'y'> becomes
761 : ! <w'x'> and <w'x'y'> reduces to <w'^2 x'>. The discretization and the code
762 : ! used in this function will be written generally in terms of <x'y'> and
763 : ! term_wpxpyp_explicit, but also applies to <x'^2> and term_wpxp2_explicit,
764 : ! as well as to <w'x'> and term_wp2xp_explicit.
765 : !
766 : ! The explicit discretization of this term is as follows:
767 : !
768 : ! 1) Centered Discretization
769 : !
770 : ! The values of <x'y'> are found on the momentum levels, while the values of
771 : ! term_wpxpyp_explicit are found on the thermodynamic levels, which is where
772 : ! they were originally calculated by the PDF. Additionally, the values of
773 : ! rho_ds_zt are found on the thermodynamic levels, and the values of
774 : ! invrs_rho_ds_zm are found on the momentum levels. At the thermodynamic
775 : ! levels, the values of term_wpxpyp_explicit are multiplied by rho_ds_zt.
776 : ! Then, the derivative (d/dz) of that expression is taken over the central
777 : ! momentum level, where it is multiplied by -invrs_rho_ds_zm. This yields
778 : ! the desired result.
779 : !
780 : ! ---------term_wpxpyp_explicitp1-------rho_ds_ztp1------------------ t(k+1)
781 : !
782 : ! ==xpyp==d( rho_ds_zt * term_wpxpyp_explicit )/dz==invrs_rho_ds_zm== m(k)
783 : !
784 : ! ---------term_wpxpyp_explicit---------rho_ds_zt-------------------- t(k)
785 : !
786 : ! The vertical indices t(k+1), m(k), and t(k) correspond with altitudes
787 : ! zt(k+1), zm(k), and zt(k), respectively. The letter "t" is used for
788 : ! thermodynamic levels and the letter "m" is used for momentum levels.
789 : !
790 : ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
791 : !
792 : ! 2) "Upwind" Discretization.
793 : !
794 : ! The values of <x'y'> are found on the momentum levels. The values of
795 : ! term_wpxpyp_explicit are originally calculated by the PDF on the
796 : ! thermodynamic levels. They are interpolated to the intermediate momentum
797 : ! levels as term_wpxpyp_explicit_zm. Additionally, the values of rho_ds_zm
798 : ! and the values of invrs_rho_ds_zm are found on the momentum levels. The
799 : ! sign of the turbulent velocity is found on the central momentum level. At
800 : ! the momentum levels, the values of term_wpxpyp_explicit_zm are multiplied
801 : ! by rho_ds_zm. Then, the derivative (d/dz) of that expression is taken.
802 : ! When the sign of the turbulent velocity is positive, the "wind" is coming
803 : ! from below, and the derivative involves the central momentum level and the
804 : ! momentum level immediately below it. When the sign of the turbulent
805 : ! velocity is negative, the "wind" is coming from above, and the derivative
806 : ! involves the central momentum level and the momenum level immediately
807 : ! above it. After the derivative is taken, it is multiplied by
808 : ! -invrs_rho_ds_zm at the central momentum level. This yields the desired
809 : ! result.
810 : !
811 : ! The turbulent velocity for <x'y'> is <w'x'y'> / <x'y'>, which has units of
812 : ! m/s. The sign of the turbulent velocity is sgn( <w'x'y'> / <x'y'> ),
813 : ! where:
814 : !
815 : ! sgn(x) = | 1; when x >= 0
816 : ! | -1; when x < 0.
817 : !
818 : ! The sign of the turbulent velocity can also be rewritten as
819 : ! sgn( <w'x'y'> ) / sgn( <x'y'> ). When a variance (<x'^2>) is being solved
820 : ! for, y = x, and sgn( <x'^2> ) is always 1. The sign of the turbulent
821 : ! velocity reduces to simply sgn( <w'x'^2> ).
822 : !
823 : ! ---------term_wpxpyp_explicit-------------------------------------- t(k+2)
824 : !
825 : ! ========term_wpxpyp_explicit_zm(interp.)=======rho_ds_zm=========== m(k+1)
826 : !
827 : ! ---------term_wpxpyp_explicit-------------------------------------- t(k+1)
828 : !
829 : ! =xpyp===term_wpxpyp_explicit_zm(interp.)=rho_ds_zm=invrs_rho_ds_zm= m(k)
830 : !
831 : ! ---------term_wpxpyp_explicit-------------------------------------- t(k)
832 : !
833 : ! ========term_wpxpyp_explicit_zm(interp.)=======rho_ds_zm=========== m(k-1)
834 : !
835 : ! ---------term_wpxpyp_explicit-------------------------------------- t(k-1)
836 : !
837 : ! The vertical indices t(k+2), m(k+1), t(k+1), m(k), t(k), m(k-1), and
838 : ! t(k-1) correspond with altitudes zt(k+2), zm(k+1), zt(k+1), zm(k), zt(k),
839 : ! zm(k-1), and zt(k-1), respectively. The letter "t" is used for
840 : ! thermodynamic levels and the letter "m" is used for momentum levels.
841 : !
842 : ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) ); and
843 : ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) ).
844 :
845 : ! References:
846 : !-----------------------------------------------------------------------
847 :
848 : use grid_class, only: &
849 : grid ! Type
850 :
851 : use constants_clubb, only: &
852 : zero ! Variable(s)
853 :
854 : use clubb_precision, only: &
855 : core_rknd ! Variable(s)
856 :
857 : implicit none
858 :
859 : ! ----------------------------- Input Variables -----------------------------
860 : integer, intent(in) :: &
861 : nz, &
862 : ngrdcol
863 :
864 : type (grid), target, intent(in) :: gr
865 :
866 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
867 : term_wpxpyp_explicit, & ! RHS: <w'x'y'> eq; t-levs [m/s(x un)(y un)]
868 : rho_ds_zt, & ! Dry, static density at t-levs [kg/m^3]
869 : invrs_rho_ds_zm ! Inv dry, static density at m-levs [m^3/kg]
870 :
871 : logical, intent(in) :: &
872 : l_upwind_xpyp_turbulent_adv ! Flag to use "upwind" discretization
873 :
874 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
875 : sgn_turbulent_vel, & ! Sign of the turbulent velocity [-]
876 : term_wpxpyp_explicit_zm, & ! term_wpxpyp_expl. zm [m/s(x un)(y un)]
877 : rho_ds_zm ! Dry, static density at m-levs [kg/m^3]
878 :
879 : ! ----------------------------- Return Variable -----------------------------
880 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
881 : rhs_ta ! RHS portion of xpyp turbulent advection [(x un)(y un)/s]
882 :
883 : ! ----------------------------- Local Variables -----------------------------
884 : integer :: i, k ! Vertical level index
885 :
886 : ! ----------------------------- Begin Code -----------------------------
887 :
888 : !$acc data copyin( gr, gr%invrs_dzm, gr%invrs_dzt, &
889 : !$acc term_wpxpyp_explicit, rho_ds_zt, invrs_rho_ds_zm, &
890 : !$acc sgn_turbulent_vel, term_wpxpyp_explicit_zm, rho_ds_zm ) &
891 : !$acc copyout( rhs_ta )
892 :
893 :
894 : ! Set lower boundary value to 0
895 : !$acc parallel loop gang vector default(present)
896 745973280 : do i = 1, ngrdcol
897 745973280 : rhs_ta(i,1) = zero
898 : end do
899 : !$acc end parallel loop
900 :
901 44675280 : if ( .not. l_upwind_xpyp_turbulent_adv ) then
902 :
903 : ! Centered discretization.
904 : !$acc parallel loop gang vector collapse(2) default(present)
905 0 : do k = 2, nz-1, 1
906 0 : do i = 1, ngrdcol
907 :
908 0 : rhs_ta(i,k) &
909 : = - invrs_rho_ds_zm(i,k) &
910 0 : * gr%invrs_dzm(i,k) * ( rho_ds_zt(i,k+1) * term_wpxpyp_explicit(i,k+1) &
911 0 : - rho_ds_zt(i,k) * term_wpxpyp_explicit(i,k) )
912 : end do
913 : end do ! k = 2, nz-1, 1
914 : !$acc end parallel loop
915 :
916 : else ! l_upwind_xpyp_turbulent_adv
917 :
918 : ! "Upwind" discretization
919 : !$acc parallel loop gang vector collapse(2) default(present)
920 3752723520 : do k = 2, nz-1, 1
921 61960457520 : do i = 1, ngrdcol
922 :
923 61915782240 : if ( sgn_turbulent_vel(i,k) > zero ) then
924 :
925 : ! The "wind" is blowing upward.
926 : rhs_ta(i,k) &
927 : = - invrs_rho_ds_zm(i,k) &
928 0 : * gr%invrs_dzt(i,k) &
929 : * ( rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) &
930 12366034315 : - rho_ds_zm(i,k-1) * term_wpxpyp_explicit_zm(i,k-1) )
931 :
932 : else ! sgn_turbulent_vel < 0
933 :
934 : ! The "wind" is blowing downward.
935 : rhs_ta(i,k) &
936 : = - invrs_rho_ds_zm(i,k) &
937 0 : * gr%invrs_dzt(i,k+1) &
938 45841699685 : * ( rho_ds_zm(i,k+1) * term_wpxpyp_explicit_zm(i,k+1) &
939 91683399370 : - rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) )
940 :
941 : endif ! sgn_turbulent_vel
942 :
943 : end do
944 : end do ! k = 2, nz-1, 1
945 : !$acc end parallel loop
946 :
947 : end if
948 :
949 : ! Set upper boundary value to 0
950 : !$acc parallel loop gang vector default(present)
951 745973280 : do i = 1, ngrdcol
952 745973280 : rhs_ta(i,nz) = zero
953 : end do
954 : !$acc end parallel loop
955 :
956 : !$acc end data
957 :
958 44675280 : return
959 :
960 : end subroutine xpyp_term_ta_pdf_rhs
961 :
962 : !=============================================================================
963 0 : subroutine xpyp_term_ta_pdf_rhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
964 0 : term_wpxpyp_explicit_zm, & ! Intent(in)
965 0 : invrs_rho_ds_zm, & ! Intent(in)
966 0 : sgn_turbulent_vel, & ! Intent(in)
967 0 : rho_ds_zm, & ! Intent(in)
968 0 : rhs_ta ) ! Intent(out)
969 : ! Description:
970 : ! This subroutine intends to add godunov upwind difference scheme based
971 : ! on xpyp_term_ta_pdf_rhs. The revisions are maded to use the Godunov-like
972 : ! upwind scheme for the vertical discretization.
973 : ! This subroutine returns an array of values for every grid level.
974 : !
975 : ! Optional Arguements:
976 : ! The optional arguements can be used to override the default indices.
977 : ! from_level - low index, default 2
978 : ! to level - high index, default gr%nz-1
979 : !
980 : ! Notes:
981 : ! This subroutine exists for testing of Godunov-like upwind scheme.
982 : ! THIS SUBROUTINE DOES NOT HANDLE BOUNDARY CONDITIONS AND SETS THEM TO 0
983 : !--------------------------------------------------------------------------------------------
984 : use clubb_precision, only: &
985 : core_rknd ! Variable(s)
986 :
987 : use grid_class, only: &
988 : grid ! Type
989 :
990 : implicit none
991 :
992 : !------------------- Input Variables -------------------
993 : integer, intent(in) :: &
994 : nz, &
995 : ngrdcol
996 :
997 : type (grid), target, intent(in) :: gr
998 :
999 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1000 : term_wpxpyp_explicit_zm, & ! RHS: <w'x'y'> eq; m-lev(k) [m/s(x un)(y un)]
1001 : invrs_rho_ds_zm, & ! Inv dry, static density at m-lev (k) [m^3/kg]
1002 : sgn_turbulent_vel, & ! Sign of the turbulent velocity [-]
1003 : rho_ds_zm ! Dry, static density at m-lev (k) [kg/m^3]
1004 :
1005 : !------------------- Output Variables -------------------
1006 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1007 : rhs_ta
1008 :
1009 : !---------------- Local Variables -------------------
1010 : integer :: &
1011 : i, k ! Loop variable for current grid level
1012 :
1013 : !---------------- Begin Code -------------------
1014 :
1015 : !$acc data copyin( gr, gr%invrs_dzm, &
1016 : !$acc term_wpxpyp_explicit_zm, invrs_rho_ds_zm, &
1017 : !$acc sgn_turbulent_vel, rho_ds_zm ) &
1018 : !$acc copyout( rhs_ta )
1019 :
1020 : ! Set lower boundary value to 0
1021 : !$acc parallel loop gang vector default(present)
1022 0 : do i = 1, ngrdcol
1023 0 : rhs_ta(i,1) = 0.0_core_rknd
1024 : end do
1025 : !$acc end parallel loop
1026 :
1027 : !$acc parallel loop gang vector collapse(2) default(present)
1028 0 : do k = 2, nz-1
1029 0 : do i = 1, ngrdcol
1030 0 : rhs_ta(i,k) = - invrs_rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
1031 0 : * ( min(0.0_core_rknd, sgn_turbulent_vel(i,k+1)) &
1032 : * rho_ds_zm(i,k+1) * term_wpxpyp_explicit_zm(i,k+1) &
1033 : + max(0.0_core_rknd, sgn_turbulent_vel(i,k+1)) &
1034 : * rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) &
1035 : - min(0.0_core_rknd, sgn_turbulent_vel(i,k))&
1036 : * rho_ds_zm(i,k) * term_wpxpyp_explicit_zm(i,k) &
1037 : - max(0.0_core_rknd, sgn_turbulent_vel(i,k)) &
1038 0 : * rho_ds_zm(i,k-1) * term_wpxpyp_explicit_zm(i,k-1) &
1039 0 : )
1040 : end do
1041 : end do
1042 : !$acc end parallel loop
1043 :
1044 : ! Set upper boundary value to 0
1045 : !$acc parallel loop gang vector default(present)
1046 0 : do i = 1, ngrdcol
1047 0 : rhs_ta(i,nz) = 0.0_core_rknd
1048 : end do
1049 : !$acc end parallel loop
1050 :
1051 : !$acc end data
1052 :
1053 0 : return
1054 :
1055 : end subroutine xpyp_term_ta_pdf_rhs_godunov
1056 :
1057 : !===============================================================================
1058 :
1059 : end module turbulent_adv_pdf
|