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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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