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

Last change on this file since 1547 was 1375, checked in by raasch, 11 years ago

last commit documented

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