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

Last change on this file since 1745 was 1683, checked in by knoop, 9 years ago

last commit documented

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