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

Last change on this file since 1322 was 1321, checked in by raasch, 10 years ago

last commit documented

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