source: palm/trunk/SOURCE/diffusion_e.f90 @ 2021

Last change on this file since 2021 was 2001, checked in by knoop, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 25.5 KB
Line 
1!> @file diffusion_e.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_e.f90 2001 2016-08-20 18:41:22Z suehring $
27!
28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
31! 1873 2016-04-18 14:50:06Z maronga
32! Module renamed (removed _mod)
33!
34!
35! 1850 2016-04-08 13:29:27Z maronga
36! Module renamed
37! Adapted for modularization of microphysics
38!
39! 1831 2016-04-07 13:15:51Z hoffmann
40! turbulence renamed collision_turbulence
41!
42! 1682 2015-10-07 23:56:08Z knoop
43! Code annotations made doxygen readable
44!
45! 1374 2014-04-25 12:55:07Z raasch
46! missing variables added to ONLY list
47! rif removed from acc-present-list
48!
49! 1340 2014-03-25 19:45:13Z kanani
50! REAL constants defined as wp-kind
51!
52! 1320 2014-03-20 08:40:49Z raasch
53! ONLY-attribute added to USE-statements,
54! kind-parameters added to all INTEGER and REAL declaration statements,
55! kinds are defined in new module kinds,
56! revision history before 2012 removed,
57! comment fields (!:) to be used for variable explanations added to
58! all variable declaration statements
59!
60! 1257 2013-11-08 15:18:40Z raasch
61! openacc loop and loop vector clauses removed
62!
63! 1179 2013-06-14 05:57:58Z raasch
64! use_reference renamed use_single_reference_value
65!
66! 1171 2013-05-30 11:27:45Z raasch
67! use_reference-case activated in accelerator version
68!
69! 1128 2013-04-12 06:19:32Z raasch
70! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
71! j_north
72!
73! 1065 2012-11-22 17:42:36Z hoffmann
74! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
75! effects of turbulence on autoconversion and accretion in two-moments cloud
76! physics scheme).
77!
78! 1036 2012-10-22 13:43:42Z raasch
79! code put under GPL (PALM 3.9)
80!
81! 1015 2012-09-27 09:23:24Z raasch
82! accelerator version (*_acc) added,
83! adjustment of mixing length to the Prandtl mixing length at first grid point
84! above ground removed
85!
86! 1010 2012-09-20 07:59:54Z raasch
87! cpp switch __nopointer added for pointer free version
88!
89! 1001 2012-09-13 14:08:46Z raasch
90! most arrays comunicated by module instead of parameter list
91!
92! 825 2012-02-19 03:03:44Z raasch
93! wang_collision_kernel renamed wang_kernel
94!
95! Revision 1.1  1997/09/19 07:40:24  raasch
96! Initial revision
97!
98!
99! Description:
100! ------------
101!> Diffusion- and dissipation terms for the TKE
102!------------------------------------------------------------------------------!
103 MODULE diffusion_e_mod
104 
105
106    PRIVATE
107    PUBLIC diffusion_e, diffusion_e_acc
108   
109
110    INTERFACE diffusion_e
111       MODULE PROCEDURE diffusion_e
112       MODULE PROCEDURE diffusion_e_ij
113    END INTERFACE diffusion_e
114 
115    INTERFACE diffusion_e_acc
116       MODULE PROCEDURE diffusion_e_acc
117    END INTERFACE diffusion_e_acc
118
119 CONTAINS
120
121
122!------------------------------------------------------------------------------!
123! Description:
124! ------------
125!> Call for all grid points
126!------------------------------------------------------------------------------!
127    SUBROUTINE diffusion_e( var, var_reference )
128
129       USE arrays_3d,                                                          &
130           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
131           
132       USE control_parameters,                                                 &
133           ONLY:  atmos_ocean_sign, g, use_single_reference_value, &
134                  wall_adjustment, wall_adjustment_factor
135
136       USE grid_variables,                                                     &
137           ONLY:  ddx2, ddy2
138           
139       USE indices,                                                            &
140           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
141                  nzt
142           
143       USE kinds
144
145       USE microphysics_mod,                                                   &
146           ONLY:  collision_turbulence
147
148       USE particle_attributes,                                                &
149           ONLY:  use_sgs_for_particles, wang_kernel
150
151       IMPLICIT NONE
152
153       INTEGER(iwp) ::  i              !<
154       INTEGER(iwp) ::  j              !<
155       INTEGER(iwp) ::  k              !<
156       REAL(wp)     ::  dvar_dz        !<
157       REAL(wp)     ::  l_stable       !<
158       REAL(wp)     ::  var_reference  !<
159
160#if defined( __nopointer )
161       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
162#else
163       REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
164#endif
165       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !<
166       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  l            !<
167       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  ll           !<
168 
169
170!
171!--    This if clause must be outside the k-loop because otherwise
172!--    runtime errors occur with -C hopt on NEC
173       IF ( use_single_reference_value )  THEN
174
175          DO  i = nxl, nxr
176             DO  j = nys, nyn
177                DO  k = nzb_s_inner(j,i)+1, nzt
178!
179!--                Calculate the mixing length (for dissipation)
180                   dvar_dz = atmos_ocean_sign * &
181                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
182                   IF ( dvar_dz > 0.0_wp ) THEN
183                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
184                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
185                   ELSE
186                      l_stable = l_grid(k)
187                   ENDIF
188!
189!--                Adjustment of the mixing length
190                   IF ( wall_adjustment )  THEN
191                      l(k,j)  = MIN( wall_adjustment_factor *          &
192                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
193                                     l_grid(k), l_stable )
194                      ll(k,j) = MIN( wall_adjustment_factor *          &
195                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
196                                     l_grid(k) )
197                   ELSE
198                      l(k,j)  = MIN( l_grid(k), l_stable )
199                      ll(k,j) = l_grid(k)
200                   ENDIF
201
202                ENDDO
203             ENDDO
204
205!
206!--          Calculate the tendency terms
207             DO  j = nys, nyn
208                DO  k = nzb_s_inner(j,i)+1, nzt
209
210                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
211                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
212
213                    tend(k,j,i) = tend(k,j,i)                                  &
214                                        + (                                    &
215                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
216                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
217                                          ) * ddx2                             &
218                                        + (                                    &
219                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
220                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
221                                          ) * ddy2                             &
222                                        + (                                    &
223               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
224             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
225                                          ) * ddzw(k)                          &
226                             - dissipation(k,j)
227
228                ENDDO
229             ENDDO
230
231!
232!--          Store dissipation if needed for calculating the sgs particle
233!--          velocities
234             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
235                  collision_turbulence )  THEN
236                DO  j = nys, nyn
237                   DO  k = nzb_s_inner(j,i)+1, nzt
238                      diss(k,j,i) = dissipation(k,j)
239                   ENDDO
240                ENDDO
241             ENDIF
242
243          ENDDO
244
245       ELSE
246
247          DO  i = nxl, nxr
248             DO  j = nys, nyn
249                DO  k = nzb_s_inner(j,i)+1, nzt
250!
251!--                Calculate the mixing length (for dissipation)
252                   dvar_dz = atmos_ocean_sign * &
253                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
254                   IF ( dvar_dz > 0.0_wp ) THEN
255                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
256                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
257                   ELSE
258                      l_stable = l_grid(k)
259                   ENDIF
260!
261!--                Adjustment of the mixing length
262                   IF ( wall_adjustment )  THEN
263                      l(k,j)  = MIN( wall_adjustment_factor *          &
264                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
265                                     l_grid(k), l_stable )
266                      ll(k,j) = MIN( wall_adjustment_factor *          &
267                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
268                                     l_grid(k) )
269                   ELSE
270                      l(k,j)  = MIN( l_grid(k), l_stable )
271                      ll(k,j) = l_grid(k)
272                   ENDIF
273
274                ENDDO
275             ENDDO
276
277!
278!--          Calculate the tendency terms
279             DO  j = nys, nyn
280                DO  k = nzb_s_inner(j,i)+1, nzt
281
282                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
283                                             e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
284
285                    tend(k,j,i) = tend(k,j,i)                                  &
286                                        + (                                    &
287                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
288                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
289                                          ) * ddx2                             &
290                                        + (                                    &
291                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
292                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
293                                          ) * ddy2                             &
294                                        + (                                    &
295               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
296             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
297                                          ) * ddzw(k)                          &
298                             - dissipation(k,j)
299
300                ENDDO
301             ENDDO
302
303!
304!--          Store dissipation if needed for calculating the sgs particle
305!--          velocities
306             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
307                  collision_turbulence )  THEN
308                DO  j = nys, nyn
309                   DO  k = nzb_s_inner(j,i)+1, nzt
310                      diss(k,j,i) = dissipation(k,j)
311                   ENDDO
312                ENDDO
313             ENDIF
314
315          ENDDO
316
317       ENDIF
318
319!
320!--    Boundary condition for dissipation
321       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  &
322            collision_turbulence )  THEN
323          DO  i = nxl, nxr
324             DO  j = nys, nyn
325                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
326             ENDDO
327          ENDDO
328       ENDIF
329
330    END SUBROUTINE diffusion_e
331
332
333!------------------------------------------------------------------------------!
334! Description:
335! ------------
336!> Call for all grid points - accelerator version
337!------------------------------------------------------------------------------!
338    SUBROUTINE diffusion_e_acc( var, var_reference )
339
340       USE arrays_3d,                                                          &
341           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
342         
343       USE control_parameters,                                                 &
344           ONLY:  atmos_ocean_sign, g, use_single_reference_value,             &
345                  wall_adjustment, wall_adjustment_factor
346
347       USE grid_variables,                                                     &
348           ONLY:  ddx2, ddy2
349           
350       USE indices,                                                            &
351           ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,   &
352                  nzb, nzb_s_inner, nzt
353           
354       USE kinds
355
356       USE microphysics_mod,                                                   &
357           ONLY:  collision_turbulence
358
359       USE particle_attributes,                                                &
360           ONLY:  use_sgs_for_particles, wang_kernel
361
362       IMPLICIT NONE
363
364       INTEGER(iwp) ::  i              !<
365       INTEGER(iwp) ::  j              !<
366       INTEGER(iwp) ::  k              !<
367       REAL(wp)     ::  dissipation    !<
368       REAL(wp)     ::  dvar_dz        !<
369       REAL(wp)     ::  l              !<
370       REAL(wp)     ::  ll             !<
371       REAL(wp)     ::  l_stable       !<
372       REAL(wp)     ::  var_reference  !<
373
374#if defined( __nopointer )
375       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
376#else
377       REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
378#endif
379
380
381!
382!--    This if clause must be outside the k-loop because otherwise
383!--    runtime errors occur with -C hopt on NEC
384       IF ( use_single_reference_value )  THEN
385
386          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
387          !$acc         present( nzb_s_inner, tend, var, zu, zw )
388          DO  i = i_left, i_right
389             DO  j = j_south, j_north
390                DO  k = 1, nzt
391
392                   IF ( k > nzb_s_inner(j,i) )  THEN
393!
394!--                   Calculate the mixing length (for dissipation)
395                      dvar_dz = atmos_ocean_sign * &
396                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
397                      IF ( dvar_dz > 0.0_wp ) THEN
398                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
399                                       SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
400                      ELSE
401                         l_stable = l_grid(k)
402                      ENDIF
403!
404!--                   Adjustment of the mixing length
405                      IF ( wall_adjustment )  THEN
406                         l  = MIN( wall_adjustment_factor *          &
407                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
408                                   l_grid(k), l_stable )
409                         ll = MIN( wall_adjustment_factor *          &
410                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
411                                   l_grid(k) )
412                      ELSE
413                         l  = MIN( l_grid(k), l_stable )
414                         ll = l_grid(k)
415                      ENDIF
416!
417!--                   Calculate the tendency terms
418                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
419                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
420
421                      tend(k,j,i) = tend(k,j,i)                                  &
422                                        + (                                    &
423                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
424                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
425                                          ) * ddx2                             &
426                                        + (                                    &
427                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
428                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
429                                          ) * ddy2                             &
430                                        + (                                    &
431               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
432             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
433                                          ) * ddzw(k)                          &
434                                  - dissipation
435
436!
437!--                   Store dissipation if needed for calculating the sgs particle
438!--                   velocities
439                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
440                           collision_turbulence )  THEN
441                         diss(k,j,i) = dissipation
442                      ENDIF
443
444                   ENDIF
445
446                ENDDO
447             ENDDO
448          ENDDO
449          !$acc end kernels
450
451       ELSE
452
453          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
454          !$acc         present( nzb_s_inner, tend, var, zu, zw )
455          DO  i = i_left, i_right
456             DO  j = j_south, j_north
457                DO  k = 1, nzt
458
459                   IF ( k > nzb_s_inner(j,i) )  THEN
460!
461!--                   Calculate the mixing length (for dissipation)
462                      dvar_dz = atmos_ocean_sign * &
463                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
464                      IF ( dvar_dz > 0.0_wp ) THEN
465                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
466                                              SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
467                      ELSE
468                         l_stable = l_grid(k)
469                      ENDIF
470!
471!--                   Adjustment of the mixing length
472                      IF ( wall_adjustment )  THEN
473                         l  = MIN( wall_adjustment_factor *          &
474                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
475                                     l_grid(k), l_stable )
476                         ll = MIN( wall_adjustment_factor *          &
477                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
478                                   l_grid(k) )
479                      ELSE
480                         l  = MIN( l_grid(k), l_stable )
481                         ll = l_grid(k)
482                      ENDIF
483!
484!--                   Calculate the tendency terms
485                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
486                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
487
488                      tend(k,j,i) = tend(k,j,i)                                &
489                                        + (                                    &
490                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
491                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
492                                          ) * ddx2                             &
493                                        + (                                    &
494                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
495                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
496                                          ) * ddy2                             &
497                                        + (                                    &
498               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
499             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
500                                          ) * ddzw(k)                          &
501                                  - dissipation
502
503!
504!--                   Store dissipation if needed for calculating the sgs
505!--                   particle  velocities
506                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
507                           collision_turbulence )  THEN
508                         diss(k,j,i) = dissipation
509                      ENDIF
510
511                   ENDIF
512
513                ENDDO
514             ENDDO
515          ENDDO
516          !$acc end kernels
517
518       ENDIF
519
520!
521!--    Boundary condition for dissipation
522       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
523            collision_turbulence )  THEN
524          !$acc kernels present( diss, nzb_s_inner )
525          DO  i = i_left, i_right
526             DO  j = j_south, j_north
527                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
528             ENDDO
529          ENDDO
530          !$acc end kernels
531       ENDIF
532
533    END SUBROUTINE diffusion_e_acc
534
535
536!------------------------------------------------------------------------------!
537! Description:
538! ------------
539!> Call for grid point i,j
540!------------------------------------------------------------------------------!
541    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
542
543       USE arrays_3d,                                                          &
544           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
545         
546       USE control_parameters,                                                 &
547           ONLY:  atmos_ocean_sign, g, use_single_reference_value,             &
548                  wall_adjustment, wall_adjustment_factor
549
550       USE grid_variables,                                                     &
551           ONLY:  ddx2, ddy2
552           
553       USE indices,                                                            &
554           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
555           
556       USE kinds
557
558       USE microphysics_mod,                                                   &
559           ONLY:  collision_turbulence
560
561       USE particle_attributes,                                                &
562           ONLY:  use_sgs_for_particles, wang_kernel
563
564       IMPLICIT NONE
565
566       INTEGER(iwp) ::  i              !<
567       INTEGER(iwp) ::  j              !<
568       INTEGER(iwp) ::  k              !<
569       REAL(wp)     ::  dvar_dz        !<
570       REAL(wp)     ::  l_stable       !<
571       REAL(wp)     ::  var_reference  !<
572
573#if defined( __nopointer )
574       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
575#else
576       REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !<
577#endif
578       REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !<
579       REAL(wp), DIMENSION(nzb+1:nzt) ::  l            !<
580       REAL(wp), DIMENSION(nzb+1:nzt) ::  ll           !<
581
582
583!
584!--    Calculate the mixing length (for dissipation)
585       DO  k = nzb_s_inner(j,i)+1, nzt
586          dvar_dz = atmos_ocean_sign * &
587                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
588          IF ( dvar_dz > 0.0_wp ) THEN
589             IF ( use_single_reference_value )  THEN
590                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
591                                     SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
592             ELSE
593                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
594                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
595             ENDIF
596          ELSE
597             l_stable = l_grid(k)
598          ENDIF
599!
600!--       Adjustment of the mixing length
601          IF ( wall_adjustment )  THEN
602             l(k)  = MIN( wall_adjustment_factor *                     &
603                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
604                          l_stable )
605             ll(k) = MIN( wall_adjustment_factor *                     &
606                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
607          ELSE
608             l(k)  = MIN( l_grid(k), l_stable )
609             ll(k) = l_grid(k)
610          ENDIF
611!
612!--       Calculate the tendency term
613          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
614                                 SQRT( e(k,j,i) ) / l(k)
615
616          tend(k,j,i) = tend(k,j,i)                                           &
617                                       + (                                    &
618                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
619                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
620                                         ) * ddx2                             &
621                                       + (                                    &
622                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
623                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
624                                         ) * ddy2                             &
625                                       + (                                    &
626              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
627            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
628                                         ) * ddzw(k)                          &
629                                       - dissipation(k)
630
631       ENDDO
632
633!
634!--    Store dissipation if needed for calculating the sgs particle velocities
635       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
636            collision_turbulence )  THEN
637          DO  k = nzb_s_inner(j,i)+1, nzt
638             diss(k,j,i) = dissipation(k)
639          ENDDO
640!
641!--       Boundary condition for dissipation
642          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
643       ENDIF
644
645    END SUBROUTINE diffusion_e_ij
646
647 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.