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

Last change on this file since 1995 was 1874, checked in by maronga, 9 years ago

last commit documented

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