Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module diffusion
5 :
6 : ! Description:
7 : ! Module diffusion computes the eddy diffusion terms for all of the
8 : ! time-tendency (prognostic) equations in the CLUBB parameterization. Most of
9 : ! the eddy diffusion terms are solved for completely implicitly, and therefore
10 : ! become part of the left-hand side of their respective equations. However,
11 : ! wp2 and wp3 have an option to use a Crank-Nicholson eddy diffusion scheme,
12 : ! which has both implicit and explicit components.
13 : !
14 : ! Function diffusion_zt_lhs handles the eddy diffusion terms for the variables
15 : ! located at thermodynamic grid levels. These variables are: wp3 and all
16 : ! hydrometeor species. The variables um and vm also use the Crank-Nicholson
17 : ! eddy-diffusion scheme for their turbulent advection term.
18 : !
19 : ! Function diffusion_zm_lhs handles the eddy diffusion terms for the variables
20 : ! located at momentum grid levels. The variables are: wprtp, wpthlp, wp2,
21 : ! rtp2, thlp2, rtpthlp, up2, vp2, wpsclrp, sclrprtp, sclrpthlp, and sclrp2.
22 :
23 : implicit none
24 :
25 : private ! Default Scope
26 :
27 : public :: diffusion_zt_lhs, &
28 : diffusion_cloud_frac_zt_lhs, &
29 : diffusion_zm_lhs
30 :
31 : contains
32 :
33 : !=============================================================================
34 17870112 : subroutine diffusion_zt_lhs( nz, ngrdcol, gr, K_zm, K_zt, nu, & ! In
35 17870112 : invrs_rho_ds_zt, rho_ds_zm, & ! In
36 17870112 : lhs ) ! Out
37 :
38 : ! Description:
39 : ! Vertical eddy diffusion of var_zt: implicit portion of the code.
40 : !
41 : ! The variable "var_zt" stands for a variable that is located at
42 : ! thermodynamic grid levels.
43 : !
44 : ! The d(var_zt)/dt equation contains an eddy diffusion term:
45 : !
46 : ! + d [ ( K_zm + nu ) * d(var_zt)/dz ] / dz.
47 : !
48 : ! This term is usually solved for completely implicitly, such that:
49 : !
50 : ! + d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
51 : !
52 : ! However, when a Crank-Nicholson scheme is used, the eddy diffusion term
53 : ! has both implicit and explicit components, such that:
54 : !
55 : ! + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz
56 : ! + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t) )/dz ] / dz;
57 : !
58 : ! for which the implicit component is:
59 : !
60 : ! + (1/2) * d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
61 : !
62 : ! Note: When the implicit term is brought over to the left-hand side,
63 : ! the sign is reversed and the leading "+" in front of the term
64 : ! is changed to a "-".
65 : !
66 : ! Timestep index (t) stands for the index of the current timestep, while
67 : ! timestep index (t+1) stands for the index of the next timestep, which is
68 : ! being advanced to in solving the d(var_zt)/dt equation.
69 : !
70 : ! The implicit portion of this term is discretized as follows:
71 : !
72 : ! The values of var_zt are found on the thermodynamic levels, while the
73 : ! values of K_zm are found on the momentum levels. The derivatives (d/dz)
74 : ! of var_zt are taken over the intermediate momentum levels. At the
75 : ! intermediate momentum levels, d(var_zt)/dz is multiplied by ( K_zm + nu ).
76 : ! Then, the derivative of the whole mathematical expression is taken over
77 : ! the central thermodynamic level, which yields the desired result.
78 : !
79 : ! --var_zt------------------------------------------------- t(k+1)
80 : !
81 : ! ==========d(var_zt)/dz==(K_zm+nu)======================== m(k)
82 : !
83 : ! --var_zt-------------------d[(K_zm+nu)*d(var_zt)/dz]/dz-- t(k)
84 : !
85 : ! ==========d(var_zt)/dz==(K_zm+nu)======================== m(k-1)
86 : !
87 : ! --var_zt------------------------------------------------- t(k-1)
88 : !
89 : ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
90 : ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
91 : ! The letter "t" is used for thermodynamic levels and the letter "m" is used
92 : ! for momentum levels.
93 : !
94 : ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) )
95 : ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
96 : ! invrs_dzm(k-1) = 1 / ( zt(k) - zt(k-1) )
97 : !
98 : ! Note: This function only computes the general implicit form:
99 : ! + d [ ( K_zm + nu ) * d( var_zt(t+1) )/dz ] / dz.
100 : ! For a Crank-Nicholson scheme, the left-hand side result of this
101 : ! function will have to be multiplied by (1/2). For a
102 : ! Crank-Nicholson scheme, the right-hand side (explicit) component
103 : ! needs to be computed by multiplying the left-hand side results by
104 : ! (1/2), reversing the sign on each left-hand side element, and then
105 : ! multiplying each element by the appropriate var_zt(t) value from
106 : ! the appropriate vertical level.
107 : !
108 : !
109 : ! Boundary Conditions:
110 : !
111 : ! 1) Zero-flux boundary conditions.
112 : ! This function is set up to use zero-flux boundary conditions at both
113 : ! the lower boundary level and the upper boundary level. The flux, F,
114 : ! is the amount of var_zt flowing normal through the boundary per unit
115 : ! time per unit surface area. The derivative of the flux effects the
116 : ! time-tendency of var_zt, such that:
117 : !
118 : ! d(var_zt)/dt = -dF/dz.
119 : !
120 : ! For the 2nd-order eddy-diffusion term, +d[(K_zm+nu)*d(var_zt)/dz]/dz,
121 : ! the flux is:
122 : !
123 : ! F = -(K_zm+nu)*d(var_zt)/dz.
124 : !
125 : ! In order to have zero-flux boundary conditions, the derivative of
126 : ! var_zt, d(var_zt)/dz, needs to equal 0 at both the lower boundary and
127 : ! the upper boundary.
128 : !
129 : ! In order to discretize the lower boundary condition, consider a new
130 : ! level outside the model (thermodynamic level 0) just below the lower
131 : ! boundary level (thermodynamic level 1). The value of var_zt at the
132 : ! level just outside the model is defined to be the same as the value of
133 : ! var_zt at the lower boundary level. Therefore, the value of
134 : ! d(var_zt)/dz between the level just outside the model and the lower
135 : ! boundary level is 0, satisfying the zero-flux boundary condition. The
136 : ! other value for d(var_zt)/dz (between thermodynamic level 2 and
137 : ! thermodynamic level 1) is taken over the intermediate momentum level
138 : ! (momentum level 1), where it is multiplied by the factor
139 : ! ( K_zm(1) + nu ). Then, the derivative of the whole expression is
140 : ! taken over the central thermodynamic level.
141 : !
142 : ! -var_zt(2)-------------------------------------------- t(2)
143 : !
144 : ! ==========d(var_zt)/dz==(K_zm(1)+nu)================== m(1)
145 : !
146 : ! -var_zt(1)---------------d[(K_zm+nu)*d(var_zt)/dz]/dz- t(1) Boundary
147 : !
148 : ! [d(var_zt)/dz = 0]
149 : !
150 : ! -[var_zt(0) = var_zt(1)]-----(level outside model)---- t(0)
151 : !
152 : ! The result is dependent only on values of K_zm found at momentum
153 : ! level 1 and values of var_zt found at thermodynamic levels 1 and 2.
154 : ! Thus, it only affects 2 diagonals on the left-hand side matrix.
155 : !
156 : ! The same method can be used to discretize the upper boundary by
157 : ! considering a new level outside the model just above the upper boundary
158 : ! level.
159 : !
160 : ! 2) Fixed-point boundary conditions.
161 : ! Many equations in the model use fixed-point boundary conditions rather
162 : ! than zero-flux boundary conditions. This means that the value of
163 : ! var_zt stays the same over the course of the timestep at the lower
164 : ! boundary, as well as at the upper boundary.
165 : !
166 : ! In order to discretize the boundary conditions for equations requiring
167 : ! fixed-point boundary conditions, either:
168 : ! a) in the parent subroutine or function (that calls this function),
169 : ! loop over all vertical levels from the second-lowest to the
170 : ! second-highest, ignoring the boundary levels. Then set the values
171 : ! at the boundary levels in the parent subroutine; or
172 : ! b) in the parent subroutine or function, loop over all vertical levels
173 : ! and then overwrite the results at the boundary levels.
174 : !
175 : ! Either way, at the boundary levels, an array with a value of 1 at the
176 : ! main diagonal on the left-hand side and with values of 0 at all other
177 : ! diagonals on the left-hand side will preserve the right-hand side value
178 : ! at that level, thus satisfying the fixed-point boundary conditions.
179 : !
180 : !
181 : ! Conservation Properties:
182 : !
183 : ! When zero-flux boundary conditions are used, this technique of
184 : ! discretizing the eddy diffusion term leads to conservative differencing.
185 : ! When conservative differencing is in place, the column totals for each
186 : ! column in the left-hand side matrix (for the eddy diffusion term) should
187 : ! be equal to 0. This ensures that the total amount of the quantity var_zt
188 : ! over the entire vertical domain is being conserved, meaning that nothing
189 : ! is lost due to diffusional effects.
190 : !
191 : ! To see that this conservation law is satisfied, compute the eddy diffusion
192 : ! of var_zt and integrate vertically. In discretized matrix notation (where
193 : ! "i" stands for the matrix column and "j" stands for the matrix row):
194 : !
195 : ! 0 = Sum_j Sum_i ( 1/invrs_dzt )_i
196 : ! ( invrs_dzt * ((K_zm+nu)*invrs_dzm) )_ij (var_zt)_j.
197 : !
198 : ! The left-hand side matrix, ( invrs_dzt * ((K_zm+nu)*invrs_dzm) )_ij, is
199 : ! partially written below. The sum over i in the above equation removes
200 : ! invrs_dzt everywhere from the matrix below. The sum over j leaves the
201 : ! column totals that are desired.
202 : !
203 : ! Left-hand side matrix contributions from eddy diffusion term; first four
204 : ! vertical levels:
205 : !
206 : ! -------------------------------------------------------------------------->
207 : !k=1 | +invrs_dzt(k) -invrs_dzt(k) 0
208 : ! | *(K_zm(k)+nu) *(K_zm(k)+nu)
209 : ! | *invrs_dzm(k) *invrs_dzm(k)
210 : ! |
211 : !k=2 | -invrs_dzt(k) +invrs_dzt(k) -invrs_dzt(k)
212 : ! | *(K_zm(k-1)+nu) *[ (K_zm(k)+nu) *(K_zm(k)+nu)
213 : ! | *invrs_dzm(k-1) *invrs_dzm(k) *invrs_dzm(k)
214 : ! | +(K_zm(k-1)+nu)
215 : ! | *invrs_dzm(k-1) ]
216 : ! |
217 : !k=3 | 0 -invrs_dzt(k) +invrs_dzt(k)
218 : ! | *(K_zm(k-1)+nu) *[ (K_zm(k)+nu)
219 : ! | *invrs_dzm(k-1) *invrs_dzm(k)
220 : ! | +(K_zm(k-1)+nu)
221 : ! | *invrs_dzm(k-1) ]
222 : ! |
223 : !k=4 | 0 0 -invrs_dzt(k)
224 : ! | *(K_zm(k-1)+nu)
225 : ! | *invrs_dzm(k-1)
226 : ! \ /
227 : !
228 : ! Note: The superdiagonal term from level 3 and both the main diagonal and
229 : ! superdiagonal terms from level 4 are not shown on this diagram.
230 : !
231 : ! Note: The matrix shown is a tridiagonal matrix. For a band diagonal
232 : ! matrix (with 5 diagonals), there would be an extra row between each
233 : ! of the rows shown and an extra column between each of the columns
234 : ! shown. However, for the purposes of the var_zt eddy diffusion
235 : ! term, those extra row and column values are all 0, and the
236 : ! conservation properties of the matrix aren't effected.
237 : !
238 : ! If fixed-point boundary conditions are used, the matrix entries at
239 : ! level 1 (k=1) read: 1 0 0; which means that conservative differencing
240 : ! is not in play. The total amount of var_zt over the entire vertical
241 : ! domain is not being conserved, as amounts of var_zt may be fluxed out
242 : ! through the upper boundary or lower boundary through the effects of
243 : ! diffusion.
244 : !
245 : ! Brian Griffin. April 26, 2008.
246 :
247 : ! References:
248 : ! None
249 : !-----------------------------------------------------------------------
250 :
251 : use grid_class, only: &
252 : grid, & ! Type
253 : ddzm ! Procedure
254 :
255 : use constants_clubb, only: &
256 : zero ! Constant(s)
257 :
258 : use clubb_precision, only: &
259 : core_rknd ! Variable(s)
260 :
261 : use model_flags, only: &
262 : l_upwind_Kh_dp_term
263 :
264 : implicit none
265 :
266 : ! Constant parameters
267 : integer, parameter :: &
268 : kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
269 : k_tdiag = 2, & ! Thermodynamic main diagonal index.
270 : km1_tdiag = 3 ! Thermodynamic subdiagonal index.
271 :
272 : !------------------------ Input Variables ------------------------
273 : integer, intent(in) :: &
274 : nz, &
275 : ngrdcol
276 :
277 : type (grid), target, intent(in) :: gr
278 :
279 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
280 : K_zm, & ! Coef. of eddy diffusivity at momentum levels [m^2/s]
281 : K_zt, & ! Coef. of eddy diffusivity at thermo. levels [m^2/s]
282 : rho_ds_zm, & ! Dry statis density on momentum levels [kg]
283 : invrs_rho_ds_zt ! Inverse dry statis density on thermo. levels [1/kg]
284 :
285 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
286 : nu ! Background constant coef. of eddy diffusivity [m^2/s]
287 :
288 : !------------------------ Return Variable ------------------------
289 : real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
290 : lhs ! LHS coefficient of diffusion term [1/s]
291 :
292 : !------------------------ Local Variable ------------------------
293 : integer :: i,k ! Vertical level index
294 :
295 : real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
296 17870112 : lhs_upwind ! LHS coefficient diffusion term due to upwind [1/s]
297 :
298 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
299 17870112 : drhoKdz_zt, &
300 35740224 : K_zm_nu, &
301 17870112 : rho_K_zm_nu, &
302 17870112 : ddzm_rho_K_zm_nu
303 :
304 : !------------------------ Begin code ------------------------
305 :
306 : !$acc data copyin( gr, gr%invrs_dzm, gr%invrs_dzt, &
307 : !$acc K_zm, K_zt, rho_ds_zm, invrs_rho_ds_zt, nu ) &
308 : !$acc copyout( lhs ) &
309 : !$acc create( lhs_upwind, drhoKdz_zt, K_zm_nu, rho_K_zm_nu, ddzm_rho_K_zm_nu )
310 :
311 : !$acc parallel loop gang vector collapse(2) default(present)
312 1536829632 : do k = 1, nz
313 25380961632 : do i = 1, ngrdcol
314 25363091520 : K_zm_nu(i,k) = K_zm(i,k) + nu(i)
315 : end do
316 : end do
317 : !$acc end parallel loop
318 :
319 : if ( l_upwind_Kh_dp_term ) then
320 :
321 : !$acc parallel loop gang vector collapse(2) default(present)
322 : do k = 1, nz
323 : do i = 1, ngrdcol
324 : rho_K_zm_nu(i,k) = rho_ds_zm(i,k) * K_zm_nu(i,k)
325 : end do
326 : end do
327 : !$acc end parallel loop
328 :
329 : ! calculate the dKh_zt/dz
330 : ddzm_rho_K_zm_nu = ddzm( nz, ngrdcol, gr, rho_K_zm_nu )
331 :
332 : !$acc parallel loop gang vector collapse(2) default(present)
333 : do k = 1, nz
334 : do i = 1, ngrdcol
335 : drhoKdz_zt(i,k) = - invrs_rho_ds_zt(i,k) * ddzm_rho_K_zm_nu(i,k)
336 : end do
337 : end do
338 : !$acc end parallel loop
339 :
340 :
341 : ! extra terms with upwind scheme
342 : ! k = 1 (bottom level); lower boundary level
343 : !$acc parallel loop gang vector default(present)
344 : do i = 1, ngrdcol
345 : lhs_upwind(kp1_tdiag,i,1) = zero
346 : lhs_upwind(k_tdiag,i,1) = zero
347 : lhs_upwind(km1_tdiag,i,1) = zero
348 : end do
349 : !$acc end parallel loop
350 :
351 : ! extra terms with upwind scheme
352 : ! k = 2 (bottom level); lower boundary level
353 : !$acc parallel loop gang vector default(present)
354 : do i = 1, ngrdcol
355 : lhs_upwind(kp1_tdiag,i,2) = + min( drhoKdz_zt(i,2) , zero ) * gr%invrs_dzm(i,2)
356 : lhs_upwind(k_tdiag,i,2) = - min( drhoKdz_zt(i,2) , zero ) * gr%invrs_dzm(i,2)
357 : lhs_upwind(km1_tdiag,i,2) = zero
358 : end do
359 : !$acc end parallel loop
360 :
361 : ! Most of the interior model; normal conditions.
362 : !$acc parallel loop gang vector collapse(2) default(present)
363 : do k = 3, nz-1, 1
364 : do i = 1, ngrdcol
365 :
366 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
367 : lhs_upwind(kp1_tdiag,i,k) = + min( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k)
368 :
369 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
370 : lhs_upwind(k_tdiag,i,k) = - min( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k) &
371 : + max( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k-1)
372 :
373 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
374 : lhs_upwind(km1_tdiag,i,k) = - max( drhoKdz_zt(i,k) , zero ) * gr%invrs_dzm(i,k-1)
375 :
376 : end do
377 : end do
378 : !$acc end parallel loop
379 :
380 : ! k = nz (top level); upper boundary level.
381 : ! Only relevant if zero-flux boundary conditions are used.
382 : !$acc parallel loop gang vector default(present)
383 : do i = 1, ngrdcol
384 : lhs_upwind(kp1_tdiag,i,nz) = zero
385 : lhs_upwind(k_tdiag,i,nz) = + max( drhoKdz_zt(i,nz) , zero ) * gr%invrs_dzm(i,nz-1)
386 : lhs_upwind(km1_tdiag,i,nz) = - max( drhoKdz_zt(i,nz) , zero ) * gr%invrs_dzm(i,nz-1)
387 : end do
388 : !$acc end parallel loop
389 :
390 : end if
391 :
392 : !$acc parallel loop gang vector default(present)
393 298389312 : do i = 1, ngrdcol
394 280519200 : lhs(kp1_tdiag,i,1) = zero
395 280519200 : lhs(k_tdiag,i,1) = zero
396 298389312 : lhs(km1_tdiag,i,1) = zero
397 : end do
398 : !$acc end parallel loop
399 :
400 : ! k = 2 (bottom level); lower boundary level.
401 : ! Only relevant if zero-flux boundary conditions are used.
402 : ! These k=2 lines currently do not have any effect on model results.
403 : ! These k=2 level of this "lhs" array is not fed into
404 : ! the final LHS matrix that will be used to solve for the next timestep.
405 :
406 : if ( .not. l_upwind_Kh_dp_term ) then
407 :
408 : !$acc parallel loop gang vector default(present)
409 298389312 : do i = 1, ngrdcol
410 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
411 561038400 : lhs(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
412 561038400 : * ( K_zm(i,2) + nu(i) ) * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
413 :
414 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
415 0 : lhs(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * invrs_rho_ds_zt(i,2) &
416 280519200 : * ( K_zm(i,2) + nu(i) ) * rho_ds_zm(i,2) * gr%invrs_dzm(i,2)
417 :
418 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
419 298389312 : lhs(km1_tdiag,i,2) = zero
420 : end do
421 : !$acc end parallel loop
422 :
423 : else
424 :
425 : !$acc parallel loop gang vector default(present)
426 : do i = 1, ngrdcol
427 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
428 : lhs(kp1_tdiag,i,2) = - gr%invrs_dzt(i,2) * ( K_zt(i,2) + nu(i) ) * gr%invrs_dzm(i,2) &
429 : + lhs_upwind(kp1_tdiag,i,2)
430 :
431 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
432 : lhs(k_tdiag,i,2) = + gr%invrs_dzt(i,2) * ( K_zt(i,2) + nu(i) ) * gr%invrs_dzm(i,2) &
433 : + lhs_upwind(k_tdiag,i,1)
434 :
435 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
436 : lhs(km1_tdiag,i,2) = zero + lhs_upwind(km1_tdiag,i,2)
437 : end do
438 : !$acc end parallel loop
439 :
440 : end if
441 :
442 : if ( .not. l_upwind_Kh_dp_term ) then
443 :
444 : ! Most of the interior model; normal conditions.
445 : !$acc parallel loop gang vector collapse(2) default(present)
446 1483219296 : do k = 3, nz-1, 1
447 24485793696 : do i = 1, ngrdcol
448 :
449 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
450 46005148800 : lhs(kp1_tdiag,i,k) &
451 0 : = - gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
452 69007723200 : * K_zm_nu(i,k) * rho_ds_zm(i,k) * gr%invrs_dzm(i,k)
453 :
454 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
455 : lhs(k_tdiag,i,k) &
456 0 : = + gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
457 0 : * ( K_zm_nu(i,k) * rho_ds_zm(i,k) * gr%invrs_dzm(i,k) &
458 23002574400 : + K_zm_nu(i,k-1) * rho_ds_zm(i,k-1) * gr%invrs_dzm(i,k-1) )
459 :
460 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
461 : lhs(km1_tdiag,i,k) &
462 0 : = - gr%invrs_dzt(i,k) * invrs_rho_ds_zt(i,k) &
463 24467923584 : * K_zm_nu(i,k-1) * rho_ds_zm(i,k-1) * gr%invrs_dzm(i,k-1)
464 : end do
465 : end do ! k = 2, nz-1, 1
466 : !$acc end parallel loop
467 :
468 : else
469 :
470 : !$acc parallel loop gang vector collapse(2) default(present)
471 : do k = 2, nz-1, 1
472 : do i = 1, ngrdcol
473 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
474 : lhs(kp1_tdiag,i,k) &
475 : = - gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * gr%invrs_dzm(i,k) &
476 : + lhs_upwind(kp1_tdiag,i,k)
477 :
478 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
479 : lhs(k_tdiag,i,k) &
480 : = + gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * ( gr%invrs_dzm(i,k) + gr%invrs_dzm(i,k-1) ) &
481 : + lhs_upwind(k_tdiag,i,k)
482 :
483 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
484 : lhs(km1_tdiag,i,k) &
485 : = - gr%invrs_dzt(i,k) * ( K_zt(i,k) + nu(i) ) * gr%invrs_dzm(i,k-1) &
486 : + lhs_upwind(km1_tdiag,i,k)
487 : end do
488 : end do ! k = 2, nz-1, 1
489 : !$acc end parallel loop
490 :
491 : end if
492 :
493 : ! k = nz (top level); upper boundary level.
494 : ! Only relevant if zero-flux boundary conditions are used.
495 :
496 : if ( .not. l_upwind_Kh_dp_term ) then
497 :
498 : !$acc parallel loop gang vector default(present)
499 298389312 : do i = 1, ngrdcol
500 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
501 280519200 : lhs(kp1_tdiag,i,nz) = zero
502 :
503 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
504 : lhs(k_tdiag,i,nz) &
505 0 : = + gr%invrs_dzt(i,nz) * invrs_rho_ds_zt(i,nz) &
506 280519200 : * ( K_zm(i,nz-1) + nu(i) ) * rho_ds_zm(i,nz-1) * gr%invrs_dzm(i,nz-1)
507 :
508 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
509 : lhs(km1_tdiag,i,nz) &
510 0 : = - gr%invrs_dzt(i,nz) * invrs_rho_ds_zt(i,nz) &
511 298389312 : * ( K_zm(i,nz-1) + nu(i) ) * rho_ds_zm(i,nz-1) * gr%invrs_dzm(i,nz-1)
512 : end do
513 : !$acc end parallel loop
514 :
515 : else
516 :
517 : !$acc parallel loop gang vector default(present)
518 : do i = 1, ngrdcol
519 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
520 : lhs(kp1_tdiag,i,nz) = zero + lhs_upwind(kp1_tdiag,i,nz)
521 :
522 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
523 : lhs(k_tdiag,i,nz) &
524 : = + gr%invrs_dzt(i,nz) * ( K_zt(i,nz) + nu(i) ) * gr%invrs_dzm(i,nz-1) &
525 : + lhs_upwind(k_tdiag,i,nz)
526 :
527 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
528 : lhs(km1_tdiag,i,nz) &
529 : = - gr%invrs_dzt(i,nz) * ( K_zt(i,nz) + nu(i) ) * gr%invrs_dzm(i,nz-1) &
530 : + lhs_upwind(km1_tdiag,i,nz)
531 : end do
532 : !$acc end parallel loop
533 :
534 : end if
535 :
536 : !$acc end data
537 :
538 17870112 : return
539 :
540 : end subroutine diffusion_zt_lhs
541 :
542 : !=============================================================================
543 0 : function diffusion_cloud_frac_zt_lhs &
544 : ( gr, K_zm, K_zmm1, cloud_frac_zt, cloud_frac_ztm1, &
545 : cloud_frac_ztp1, cloud_frac_zm, &
546 : cloud_frac_zmm1, &
547 : nu, invrs_dzmm1, invrs_dzm, invrs_dzt, level ) &
548 0 : result( lhs )
549 :
550 : ! Description:
551 : ! This function adds a weight of cloud fraction to the existing diffusion
552 : ! function for number concentration variables (e.g. cloud droplet number
553 : ! concentration). This code should be considered experimental and may
554 : ! contain bugs.
555 : ! References:
556 : ! This algorithm uses equations derived from Guo, et al. 2009.
557 : !-----------------------------------------------------------------------------
558 :
559 : use grid_class, only: &
560 : grid ! Type
561 :
562 : use clubb_precision, only: &
563 : core_rknd ! Variable(s)
564 :
565 : implicit none
566 :
567 : type (grid), target, intent(in) :: gr
568 :
569 : ! External
570 : intrinsic :: min
571 :
572 : ! Constant parameters
573 : real( kind = core_rknd ), parameter :: &
574 : cf_ratio = 10._core_rknd ! Maximum cloud-fraction coefficient applied to Kh_zm
575 :
576 : integer, parameter :: &
577 : kp1_tdiag = 1, & ! Thermodynamic superdiagonal index.
578 : k_tdiag = 2, & ! Thermodynamic main diagonal index.
579 : km1_tdiag = 3 ! Thermodynamic subdiagonal index.
580 :
581 : ! Input Variables
582 : real( kind = core_rknd ), intent(in) :: &
583 : K_zm, & ! Coef. of eddy diffusivity at mom. level (k) [m^2/s]
584 : K_zmm1, & ! Coef. of eddy diffusivity at mom. level (k-1) [m^2/s]
585 : cloud_frac_zt, & ! Cloud fraction at the thermo. level (k) [-]
586 : cloud_frac_ztm1, & ! Cloud fraction at the thermo. level (k-1) [-]
587 : cloud_frac_ztp1, & ! Cloud fraction at the thermo. level (k+1) [-]
588 : cloud_frac_zm, & ! Cloud fraction at the momentum level (k) [-]
589 : cloud_frac_zmm1, & ! Cloud fraction at the momentum level (k-1) [-]
590 : invrs_dzt, & ! Inverse of grid spacing over thermo. lev. (k) [1/m]
591 : invrs_dzm, & ! Inverse of grid spacing over mom. level (k) [1/m]
592 : invrs_dzmm1 ! Inverse of grid spacing over mom. level (k-1) [1/m]
593 :
594 : real( kind = core_rknd ), intent(in) :: &
595 : nu ! Background constant coef. of eddy diffusivity [m^2/s]
596 :
597 : integer, intent(in) :: &
598 : level ! Thermodynamic level where calculation occurs. [-]
599 :
600 : ! Return Variable
601 : real( kind = core_rknd ), dimension(3) :: lhs
602 :
603 : ! ---- Begin Code ----
604 :
605 0 : if ( level == 1 ) then
606 :
607 : ! k = 1 (bottom level); lower boundary level.
608 : ! Only relevant if zero-flux boundary conditions are used.
609 :
610 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
611 : ! lhs(kp1_tdiag) = - invrs_dzt &
612 : ! * (K_zm+nu) &
613 : ! * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
614 : lhs(kp1_tdiag) = - invrs_dzt &
615 : * (K_zm &
616 : * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
617 0 : + nu) * invrs_dzm
618 :
619 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
620 : ! lhs(k_tdiag) = + invrs_dzt &
621 : ! * (K_zm+nu) &
622 : ! * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
623 : lhs(k_tdiag) = + invrs_dzt &
624 : * (K_zm &
625 : * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
626 0 : + nu) * invrs_dzm
627 :
628 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
629 0 : lhs(km1_tdiag) = 0.0_core_rknd
630 :
631 :
632 0 : else if ( level > 1 .and. level < gr%nz ) then
633 :
634 : ! Most of the interior model; normal conditions.
635 :
636 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
637 : ! lhs(kp1_tdiag) = - invrs_dzt &
638 : ! * (K_zm+nu) &
639 : ! * ( cloud_frac_zm / cloud_frac_ztp1 ) * invrs_dzm
640 : ! lhs(kp1_tdiag) = - invrs_dzt &
641 : ! * (K_zm &
642 : ! * ( cloud_frac_zm / cloud_frac_ztp1 ) &
643 : ! + nu ) * invrs_dzm
644 : lhs(kp1_tdiag) = - invrs_dzt &
645 : * (K_zm &
646 : * min( cloud_frac_zm / cloud_frac_ztp1, cf_ratio ) &
647 0 : + nu ) * invrs_dzm
648 :
649 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
650 : ! lhs(k_tdiag) = + invrs_dzt &
651 : ! * ( ((K_zm+nu)*cloud_frac_zm)*invrs_dzm &
652 : ! + ((K_zmm1+nu)*cloud_frac_zmm1)*invrs_dzmm1 ) &
653 : ! / cloud_frac_zt
654 : ! lhs(k_tdiag) = + invrs_dzt &
655 : ! * ( nu*(invrs_dzm+invrs_dzmm1) + &
656 : ! ( ((K_zm*cloud_frac_zm)*invrs_dzm +
657 : ! (K_zmm1*cloud_frac_zmm1)*invrs_dzmm1)&
658 : ! / cloud_frac_zt &
659 : ! ) &
660 : ! )
661 : lhs(k_tdiag) = + invrs_dzt &
662 : * ( nu*(invrs_dzm+invrs_dzmm1) + &
663 : ( K_zm*invrs_dzm* &
664 : min( cloud_frac_zm / cloud_frac_zt, &
665 : cf_ratio ) &
666 : + K_zmm1*invrs_dzmm1* &
667 : min( cloud_frac_zmm1 / cloud_frac_zt, &
668 : cf_ratio ) &
669 : ) &
670 0 : )
671 :
672 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
673 : ! lhs(km1_tdiag) = - invrs_dzt * (K_zmm1+nu) * &
674 : ! ( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
675 : lhs(km1_tdiag) = - invrs_dzt &
676 : * (K_zmm1 &
677 : * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
678 : cf_ratio ) &
679 0 : + nu ) * invrs_dzmm1
680 :
681 0 : else if ( level == gr%nz ) then
682 :
683 : ! k = gr%nz (top level); upper boundary level.
684 : ! Only relevant if zero-flux boundary conditions are used.
685 :
686 : ! Thermodynamic superdiagonal: [ x var_zt(k+1,<t+1>) ]
687 0 : lhs(kp1_tdiag) = 0.0_core_rknd
688 :
689 : ! Thermodynamic main diagonal: [ x var_zt(k,<t+1>) ]
690 : ! lhs(k_tdiag) = + invrs_dzt &
691 : ! *(K_zmm1+nu) &
692 : ! *( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
693 : lhs(k_tdiag) = + invrs_dzt &
694 : * (K_zmm1 &
695 : * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
696 : cf_ratio ) &
697 0 : + nu) * invrs_dzmm1
698 :
699 : ! Thermodynamic subdiagonal: [ x var_zt(k-1,<t+1>) ]
700 : ! lhs(km1_tdiag) = - invrs_dzt * (K_zmm1+nu) * &
701 : ! ( cloud_frac_zmm1 / cloud_frac_ztm1 ) * invrs_dzmm1
702 : lhs(km1_tdiag) = - invrs_dzt &
703 : * (K_zmm1 &
704 : * min( cloud_frac_zmm1 / cloud_frac_ztm1, &
705 : cf_ratio ) &
706 0 : + nu) * invrs_dzmm1
707 :
708 : end if
709 :
710 0 : return
711 0 : end function diffusion_cloud_frac_zt_lhs
712 :
713 : !=============================================================================
714 35740224 : subroutine diffusion_zm_lhs( nz, ngrdcol, gr, K_zt, K_zm, nu, & ! In
715 35740224 : invrs_rho_ds_zm, rho_ds_zt, & ! In
716 35740224 : lhs ) ! Out
717 :
718 : ! Description:
719 : ! Vertical eddy diffusion of var_zm: implicit portion of the code.
720 : !
721 : ! The variable "var_zm" stands for a variable that is located at momentum
722 : ! grid levels.
723 : !
724 : ! The d(var_zm)/dt equation contains an eddy diffusion term:
725 : !
726 : ! + d [ ( K_zt + nu ) * d(var_zm)/dz ] / dz.
727 : !
728 : ! This term is usually solved for completely implicitly, such that:
729 : !
730 : ! + d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
731 : !
732 : ! However, when a Crank-Nicholson scheme is used, the eddy diffusion term
733 : ! has both implicit and explicit components, such that:
734 : !
735 : ! + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz
736 : ! + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t) )/dz ] / dz;
737 : !
738 : ! for which the implicit component is:
739 : !
740 : ! + (1/2) * d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
741 : !
742 : ! Note: When the implicit term is brought over to the left-hand side,
743 : ! the sign is reversed and the leading "+" in front of the term
744 : ! is changed to a "-".
745 : !
746 : ! Timestep index (t) stands for the index of the current timestep, while
747 : ! timestep index (t+1) stands for the index of the next timestep, which is
748 : ! being advanced to in solving the d(var_zm)/dt equation.
749 : !
750 : ! The implicit portion of this term is discretized as follows:
751 : !
752 : ! The values of var_zm are found on the momentum levels, while the values of
753 : ! K_zt are found on the thermodynamic levels. The derivatives (d/dz) of
754 : ! var_zm are taken over the intermediate thermodynamic levels. At the
755 : ! intermediate thermodynamic levels, d(var_zm)/dz is multiplied by
756 : ! ( K_zt + nu ). Then, the derivative of the whole mathematical expression
757 : ! is taken over the central momentum level, which yields the desired result.
758 : !
759 : ! ==var_zm================================================= m(k+1)
760 : !
761 : ! ----------d(var_zm)/dz--(K_zt+nu)------------------------ t(k+1)
762 : !
763 : ! ==var_zm===================d[(K_zt+nu)*d(var_zm)/dz]/dz== m(k)
764 : !
765 : ! ----------d(var_zm)/dz--(K_zt+nu)------------------------ t(k)
766 : !
767 : ! ==var_zm================================================= m(k-1)
768 : !
769 : ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
770 : ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
771 : ! The letter "t" is used for thermodynamic levels and the letter "m" is used
772 : ! for momentum levels.
773 : !
774 : ! invrs_dzm(k) = 1 / ( zt(k+1) - zt(k) )
775 : ! invrs_dzt(k+1) = 1 / ( zm(k+1) - zm(k) )
776 : ! invrs_dzt(k) = 1 / ( zm(k) - zm(k-1) )
777 : !
778 : ! Note: This function only computes the general implicit form:
779 : ! + d [ ( K_zt + nu ) * d( var_zm(t+1) )/dz ] / dz.
780 : ! For a Crank-Nicholson scheme, the left-hand side result of this
781 : ! function will have to be multiplied by (1/2). For a
782 : ! Crank-Nicholson scheme, the right-hand side (explicit) component
783 : ! needs to be computed by multiplying the left-hand side results by
784 : ! (1/2), reversing the sign on each left-hand side element, and then
785 : ! multiplying each element by the appropriate var_zm(t) value from
786 : ! the appropriate vertical level.
787 : !
788 : !
789 : ! Boundary Conditions:
790 : !
791 : ! 1) Zero-flux boundary conditions.
792 : ! This function is set up to use zero-flux boundary conditions at both
793 : ! the lower boundary level and the upper boundary level. The flux, F,
794 : ! is the amount of var_zm flowing normal through the boundary per unit
795 : ! time per unit surface area. The derivative of the flux effects the
796 : ! time-tendency of var_zm, such that:
797 : !
798 : ! d(var_zm)/dt = -dF/dz.
799 : !
800 : ! For the 2nd-order eddy-diffusion term, +d[(K_zt+nu)*d(var_zm)/dz]/dz,
801 : ! the flux is:
802 : !
803 : ! F = -(K_zt+nu)*d(var_zm)/dz.
804 : !
805 : ! In order to have zero-flux boundary conditions, the derivative of
806 : ! var_zm, d(var_zm)/dz, needs to equal 0 at both the lower boundary and
807 : ! the upper boundary.
808 : !
809 : ! In order to discretize the lower boundary condition, consider a new
810 : ! level outside the model (momentum level 0) just below the lower
811 : ! boundary level (momentum level 1). The value of var_zm at the level
812 : ! just outside the model is defined to be the same as the value of var_zm
813 : ! at the lower boundary level. Therefore, the value of d(var_zm)/dz
814 : ! between the level just outside the model and the lower boundary level
815 : ! is 0, satisfying the zero-flux boundary condition. The other value for
816 : ! d(var_zm)/dz (between momentum level 2 and momentum level 1) is taken
817 : ! over the intermediate thermodynamic level (thermodynamic level 2),
818 : ! where it is multiplied by the factor ( K_zt(2) + nu ). Then, the
819 : ! derivative of the whole expression is taken over the central momentum
820 : ! level.
821 : !
822 : ! =var_zm(2)============================================ m(2)
823 : !
824 : ! ----------d(var_zm)/dz==(K_zt(2)+nu)------------------ t(2)
825 : !
826 : ! =var_zm(1)===============d[(K_zt+nu)*d(var_zm)/dz]/dz= m(1) Boundary
827 : !
828 : ! ----------[d(var_zm)/dz = 0]-------------------------- t(1)
829 : !
830 : ! =[var_zm(0) = var_zm(1)]=====(level outside model)==== m(0)
831 : !
832 : ! The result is dependent only on values of K_zt found at thermodynamic
833 : ! level 2 and values of var_zm found at momentum levels 1 and 2. Thus,
834 : ! it only affects 2 diagonals on the left-hand side matrix.
835 : !
836 : ! The same method can be used to discretize the upper boundary by
837 : ! considering a new level outside the model just above the upper boundary
838 : ! level.
839 : !
840 : ! 2) Fixed-point boundary conditions.
841 : ! Many equations in the model use fixed-point boundary conditions rather
842 : ! than zero-flux boundary conditions. This means that the value of
843 : ! var_zm stays the same over the course of the timestep at the lower
844 : ! boundary, as well as at the upper boundary.
845 : !
846 : ! In order to discretize the boundary conditions for equations requiring
847 : ! fixed-point boundary conditions, either:
848 : ! a) in the parent subroutine or function (that calls this function),
849 : ! loop over all vertical levels from the second-lowest to the
850 : ! second-highest, ignoring the boundary levels. Then set the values
851 : ! at the boundary levels in the parent subroutine; or
852 : ! b) in the parent subroutine or function, loop over all vertical levels
853 : ! and then overwrite the results at the boundary levels.
854 : !
855 : ! Either way, at the boundary levels, an array with a value of 1 at the
856 : ! main diagonal on the left-hand side and with values of 0 at all other
857 : ! diagonals on the left-hand side will preserve the right-hand side value
858 : ! at that level, thus satisfying the fixed-point boundary conditions.
859 : !
860 : !
861 : ! Conservation Properties:
862 : !
863 : ! When zero-flux boundary conditions are used, this technique of
864 : ! discretizing the eddy diffusion term leads to conservative differencing.
865 : ! When conservative differencing is in place, the column totals for each
866 : ! column in the left-hand side matrix (for the eddy diffusion term) should
867 : ! be equal to 0. This ensures that the total amount of the quantity var_zm
868 : ! over the entire vertical domain is being conserved, meaning that nothing
869 : ! is lost due to diffusional effects.
870 : !
871 : ! To see that this conservation law is satisfied, compute the eddy diffusion
872 : ! of var_zm and integrate vertically. In discretized matrix notation (where
873 : ! "i" stands for the matrix column and "j" stands for the matrix row):
874 : !
875 : ! 0 = Sum_j Sum_i ( 1/invrs_dzm )_i
876 : ! ( invrs_dzm * ((K_zt+nu)*invrs_dzt) )_ij (var_zm)_j.
877 : !
878 : ! The left-hand side matrix, ( invrs_dzm * ((K_zt+nu)*invrs_dzt) )_ij, is
879 : ! partially written below. The sum over i in the above equation removes
880 : ! invrs_dzm everywhere from the matrix below. The sum over j leaves the
881 : ! column totals that are desired.
882 : !
883 : ! Left-hand side matrix contributions from eddy diffusion term; first four
884 : ! vertical levels:
885 : !
886 : ! ---------------------------------------------------------------------->
887 : !k=1 | +invrs_dzm(k) -invrs_dzm(k) 0
888 : ! | *(K_zt(k+1)+nu) *(K_zt(k+1)+nu)
889 : ! | *invrs_dzt(k+1) *invrs_dzt(k+1)
890 : ! |
891 : !k=2 | -invrs_dzm(k) +invrs_dzm(k) -invrs_dzm(k)
892 : ! | *(K_zt(k)+nu) *[ (K_zt(k+1)+nu) *(K_zt(k+1)+nu)
893 : ! | *invrs_dzt(k) *invrs_dzt(k+1) *invrs_dzt(k+1)
894 : ! | +(K_zt(k)+nu)
895 : ! | *invrs_dzt(k) ]
896 : ! |
897 : !k=3 | 0 -invrs_dzm(k) +invrs_dzm(k)
898 : ! | *(K_zt(k)+nu) *[ (K_zt(k+1)+nu)
899 : ! | *invrs_dzt(k) *invrs_dzt(k+1)
900 : ! | +(K_zt(k)+nu)
901 : ! | *invrs_dzt(k) ]
902 : ! |
903 : !k=4 | 0 0 -invrs_dzm(k)
904 : ! | *(K_zt(k)+nu)
905 : ! | *invrs_dzt(k)
906 : ! \ /
907 : !
908 : ! Note: The superdiagonal term from level 3 and both the main diagonal and
909 : ! superdiagonal terms from level 4 are not shown on this diagram.
910 : !
911 : ! Note: The matrix shown is a tridiagonal matrix. For a band diagonal
912 : ! matrix (with 5 diagonals), there would be an extra row between each
913 : ! of the rows shown and an extra column between each of the columns
914 : ! shown. However, for the purposes of the var_zm eddy diffusion
915 : ! term, those extra row and column values are all 0, and the
916 : ! conservation properties of the matrix aren't effected.
917 : !
918 : ! If fixed-point boundary conditions are used, the matrix entries at
919 : ! level 1 (k=1) read: 1 0 0; which means that conservative differencing
920 : ! is not in play. The total amount of var_zm over the entire vertical
921 : ! domain is not being conserved, as amounts of var_zm may be fluxed out
922 : ! through the upper boundary or lower boundary through the effects of
923 : ! diffusion.
924 : !
925 : ! Brian Griffin. April 26, 2008.
926 :
927 : ! References:
928 : ! None
929 : !-----------------------------------------------------------------------
930 :
931 : use grid_class, only: &
932 : grid, & ! Type
933 : ddzt ! Procedure
934 :
935 : use constants_clubb, only: &
936 : zero ! Constant(s)
937 :
938 : use clubb_precision, only: &
939 : core_rknd ! Variable(s)
940 :
941 : use model_flags, only: &
942 : l_upwind_Kh_dp_term
943 :
944 : implicit none
945 :
946 : ! Constant parameters
947 : integer, parameter :: &
948 : kp1_mdiag = 1, & ! Momentum superdiagonal index.
949 : k_mdiag = 2, & ! Momentum main diagonal index.
950 : km1_mdiag = 3 ! Momentum subdiagonal index.
951 :
952 : ! Input Variables
953 : integer, intent(in) :: &
954 : nz, &
955 : ngrdcol
956 :
957 : type (grid), target, intent(in) :: gr
958 :
959 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
960 : K_zm, &
961 : K_zt, & ! Coef. of eddy diffusivity at thermo. levels [m^2/s]
962 : rho_ds_zt, &
963 : invrs_rho_ds_zm
964 :
965 : real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
966 : nu ! Background constant coef. of eddy diffusivity [m^2/s]
967 :
968 : ! Return Variable
969 : real( kind = core_rknd ), dimension(3,ngrdcol,nz), intent(out) :: &
970 : lhs ! LHS coefficient of diffusion term [1/s]
971 :
972 : ! Local Variable
973 : integer :: i, k ! Vertical level index
974 :
975 : real( kind = core_rknd ), dimension(3,ngrdcol,nz) :: &
976 71480448 : lhs_upwind ! LHS coefficient diffusion term due to upwind [1/s]
977 :
978 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
979 71480448 : drhoKdz_zm, &
980 71480448 : rho_K_zt_nu, &
981 71480448 : ddzt_rho_K_zt_nu
982 :
983 : !------------Begin code------------------------------
984 :
985 : !$acc data copyin( gr, gr%invrs_dzt, gr%invrs_dzm, &
986 : !$acc K_zm, K_zt, rho_ds_zt, invrs_rho_ds_zm, nu ) &
987 : !$acc copyout( lhs ) &
988 : !$acc create( lhs_upwind, drhoKdz_zm, rho_K_zt_nu, ddzt_rho_K_zt_nu )
989 :
990 : !$acc parallel loop gang vector collapse(2) default(present)
991 3073659264 : do k = 1, nz
992 50761923264 : do i = 1, ngrdcol
993 50726183040 : rho_K_zt_nu(i,k) = rho_ds_zt(i,k) * ( K_zt(i,k) + nu(i) )
994 : end do
995 : end do
996 : !$acc end parallel loop
997 :
998 : ! calculate the dKh_zm/dz
999 35740224 : ddzt_rho_K_zt_nu = ddzt( nz, ngrdcol, gr, rho_K_zt_nu )
1000 :
1001 : !$acc parallel loop gang vector collapse(2) default(present)
1002 3073659264 : do k = 1, nz
1003 50761923264 : do i = 1, ngrdcol
1004 50726183040 : drhoKdz_zm(i,k) = - invrs_rho_ds_zm(i,k) * ddzt_rho_K_zt_nu(i,k)
1005 : end do
1006 : end do
1007 : !$acc end parallel loop
1008 :
1009 : ! extra terms with upwind scheme
1010 : ! k = 1 (bottom level); lowere boundary level
1011 : !$acc parallel loop gang vector default(present)
1012 596778624 : do i = 1, ngrdcol
1013 561038400 : lhs_upwind(kp1_mdiag,i,1) = + min( drhoKdz_zm(i,1) , zero ) * gr%invrs_dzt(i,2)
1014 561038400 : lhs_upwind(k_mdiag,i,1) = - min( drhoKdz_zm(i,1) , zero ) * gr%invrs_dzt(i,2)
1015 596778624 : lhs_upwind(km1_mdiag,i,1) = zero
1016 : end do
1017 : !$acc end parallel loop
1018 :
1019 : ! Most of the interior model; normal conditions.
1020 : !$acc parallel loop gang vector collapse(2) default(present)
1021 3002178816 : do k = 2, nz-1, 1
1022 49568366016 : do i = 1, ngrdcol
1023 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1024 46566187200 : lhs_upwind(kp1_mdiag,i,k) = + min( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k+1)
1025 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1026 0 : lhs_upwind(k_mdiag,i,k) = - min( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k+1) &
1027 46566187200 : + max( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k)
1028 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1029 49532625792 : lhs_upwind(km1_mdiag,i,k) = - max( drhoKdz_zm(i,k) , zero ) * gr%invrs_dzt(i,k)
1030 : end do
1031 : end do
1032 : !$acc end parallel loop
1033 :
1034 : ! k = nz (top level); upper boundary level.
1035 : ! Only relevant if zero-flux boundary conditions are used.
1036 : !$acc parallel loop gang vector default(present)
1037 596778624 : do i = 1, ngrdcol
1038 561038400 : lhs_upwind(kp1_mdiag,i,nz) = zero
1039 561038400 : lhs_upwind(k_mdiag,i,nz) = + max( drhoKdz_zm(i,nz) , zero ) * gr%invrs_dzt(i,nz)
1040 596778624 : lhs_upwind(km1_mdiag,i,nz) = - max( drhoKdz_zm(i,nz) , zero ) * gr%invrs_dzt(i,nz)
1041 : end do
1042 : !$acc end parallel loop
1043 :
1044 : ! k = 1; lower boundary level at surface.
1045 : ! Only relevant if zero-flux boundary conditions are used.
1046 : ! These k=1 lines currently do not have any effect on model results.
1047 : ! These k=1 level of this "lhs" array is not fed into
1048 : ! the final LHS matrix that will be used to solve for the next timestep.
1049 :
1050 : if ( .not. l_upwind_Kh_dp_term ) then
1051 :
1052 : !$acc parallel loop gang vector default(present)
1053 596778624 : do i = 1, ngrdcol
1054 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1055 1122076800 : lhs(kp1_mdiag,i,1) = - gr%invrs_dzm(i,1) * invrs_rho_ds_zm(i,1) &
1056 1122076800 : * ( K_zt(i,2) + nu(i) ) * rho_ds_zt(i,2) * gr%invrs_dzt(i,2)
1057 :
1058 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1059 0 : lhs(k_mdiag,i,1) = + gr%invrs_dzm(i,1) * invrs_rho_ds_zm(i,1) &
1060 561038400 : * ( K_zt(i,2) + nu(i) ) * rho_ds_zt(i,2) * gr%invrs_dzt(i,2)
1061 :
1062 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1063 596778624 : lhs(km1_mdiag,i,1) = zero
1064 : end do
1065 : !$acc end parallel loop
1066 :
1067 : else
1068 :
1069 : !$acc parallel loop gang vector default(present)
1070 : do i = 1, ngrdcol
1071 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1072 : lhs(kp1_mdiag,i,1) = - gr%invrs_dzm(i,1) * ( K_zm(i,1) + nu(i) ) * gr%invrs_dzt(i,2) &
1073 : + lhs_upwind(kp1_mdiag,i,1)
1074 :
1075 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1076 : lhs(k_mdiag,i,1) = + gr%invrs_dzm(i,1) * ( K_zm(i,1) + nu(i) ) * gr%invrs_dzt(i,2) &
1077 : + lhs_upwind(k_mdiag,i,1)
1078 :
1079 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1080 : lhs(km1_mdiag,i,1) = zero + lhs_upwind(km1_mdiag,i,1)
1081 : end do
1082 : !$acc end parallel loop
1083 :
1084 : end if
1085 :
1086 : ! Most of the interior model; normal conditions.
1087 : if ( .not. l_upwind_Kh_dp_term ) then
1088 :
1089 : !$acc parallel loop gang vector collapse(2) default(present)
1090 3002178816 : do k = 2, nz-1, 1
1091 49568366016 : do i = 1, ngrdcol
1092 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1093 93132374400 : lhs(kp1_mdiag,i,k) &
1094 0 : = - gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
1095 >13969*10^7 : * ( K_zt(i,k+1) + nu(i) ) * rho_ds_zt(i,k+1) * gr%invrs_dzt(i,k+1)
1096 :
1097 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1098 : lhs(k_mdiag,i,k) &
1099 0 : = + gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
1100 0 : * ( ( K_zt(i,k+1) + nu(i) ) * rho_ds_zt(i,k+1) * gr%invrs_dzt(i,k+1) &
1101 46566187200 : + ( K_zt(i,k) + nu(i) ) * rho_ds_zt(i,k) * gr%invrs_dzt(i,k) )
1102 :
1103 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1104 : lhs(km1_mdiag,i,k) &
1105 0 : = - gr%invrs_dzm(i,k) * invrs_rho_ds_zm(i,k) &
1106 49532625792 : * ( K_zt(i,k) + nu(i) ) * rho_ds_zt(i,k) * gr%invrs_dzt(i,k)
1107 : end do
1108 : end do
1109 : !$acc end parallel loop
1110 :
1111 : else
1112 :
1113 : !$acc parallel loop gang vector collapse(2) default(present)
1114 : do k = 2, nz-1, 1
1115 : do i = 1, ngrdcol
1116 :
1117 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1118 : lhs(kp1_mdiag,i,k) &
1119 : = - gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * gr%invrs_dzt(i,k+1) &
1120 : + lhs_upwind(kp1_mdiag,i,k)
1121 :
1122 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1123 : lhs(k_mdiag,i,k) &
1124 : = + gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * ( gr%invrs_dzt(i,k+1) + gr%invrs_dzt(i,k) ) &
1125 : + lhs_upwind(k_mdiag,i,k)
1126 :
1127 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1128 : lhs(km1_mdiag,i,k) &
1129 : = - gr%invrs_dzm(i,k) * ( K_zm(i,k) + nu(i) ) * gr%invrs_dzt(i,k) &
1130 : + lhs_upwind(km1_mdiag,i,k)
1131 :
1132 : end do
1133 : end do
1134 : !$acc end parallel loop
1135 :
1136 : end if
1137 :
1138 : ! k = nz (top level); upper boundary level.
1139 : ! Only relevant if zero-flux boundary conditions are used.
1140 :
1141 : if ( .not. l_upwind_Kh_dp_term ) then
1142 :
1143 : !$acc parallel loop gang vector default(present)
1144 596778624 : do i = 1, ngrdcol
1145 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1146 561038400 : lhs(kp1_mdiag,i,nz) = zero
1147 :
1148 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1149 : lhs(k_mdiag,i,nz) &
1150 0 : = + gr%invrs_dzm(i,nz) * invrs_rho_ds_zm(i,nz) &
1151 561038400 : * ( K_zt(i,nz) + nu(i) ) * rho_ds_zt(i,nz) * gr%invrs_dzt(i,nz)
1152 :
1153 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1154 : lhs(km1_mdiag,i,nz) &
1155 0 : = - gr%invrs_dzm(i,nz) * invrs_rho_ds_zm(i,nz) &
1156 596778624 : * ( K_zt(i,nz) + nu(i) ) * rho_ds_zt(i,nz) * gr%invrs_dzt(i,nz)
1157 : end do
1158 : !$acc end parallel loop
1159 :
1160 : else
1161 :
1162 : !$acc parallel loop gang vector default(present)
1163 : do i = 1, ngrdcol
1164 : ! Momentum superdiagonal: [ x var_zm(k+1,<t+1>) ]
1165 : lhs(kp1_mdiag,i,nz) = zero &
1166 : + lhs_upwind(kp1_mdiag,i,nz)
1167 :
1168 : ! Momentum main diagonal: [ x var_zm(k,<t+1>) ]
1169 : lhs(k_mdiag,i,nz) &
1170 : = + gr%invrs_dzm(i,nz) * ( K_zm(i,nz) + nu(i) ) * gr%invrs_dzt(i,nz) &
1171 : + lhs_upwind(k_mdiag,i,nz)
1172 :
1173 : ! Momentum subdiagonal: [ x var_zm(k-1,<t+1>) ]
1174 : lhs(km1_mdiag,i,nz) &
1175 : = - gr%invrs_dzm(i,nz) * ( K_zm(i,nz) + nu(i) ) * gr%invrs_dzt(i,nz) &
1176 : + lhs_upwind(km1_mdiag,i,nz)
1177 : end do
1178 : !$acc end parallel loop
1179 :
1180 : end if
1181 :
1182 : !$acc end data
1183 :
1184 35740224 : return
1185 :
1186 : end subroutine diffusion_zm_lhs
1187 :
1188 : !=============================================================================
1189 :
1190 : end module diffusion
|
