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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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