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

Last change on this file since 1851 was 1851, checked in by maronga, 8 years ago

last commit documented

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