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

Last change on this file since 1841 was 1832, checked in by hoffmann, 8 years ago

last commit documented

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