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

Last change on this file since 1344 was 1341, checked in by kanani, 11 years ago

last commit documented

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