source: palm/trunk/SOURCE/diffusion_e_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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