source: palm/trunk/SOURCE/diffusion_e_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago

added _mod string to several filenames to meet the naming convection for modules

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