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

Last change on this file since 1941 was 1874, checked in by maronga, 9 years ago

last commit documented

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