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

Last change on this file since 1832 was 1832, checked in by hoffmann, 8 years ago

last commit documented

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