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

Last change on this file since 2037 was 2037, checked in by knoop, 5 years ago

Anelastic approximation implemented

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