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

Last change on this file since 2119 was 2119, checked in by raasch, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 17.4 KB
RevLine 
[1873]1!> @file diffusion_e.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
4!
[2000]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.
[1036]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!
[2101]17! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
[2119]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_e.f90 2119 2017-01-17 16:51:50Z raasch $
27!
[2119]28! 2118 2017-01-17 16:38:49Z raasch
29! OpenACC version of subroutine removed
30!
[2038]31! 2037 2016-10-26 11:15:40Z knoop
32! Anelastic approximation implemented
33!
[2001]34! 2000 2016-08-20 18:09:15Z knoop
35! Forced header and separation lines into 80 columns
36!
[1874]37! 1873 2016-04-18 14:50:06Z maronga
38! Module renamed (removed _mod)
39!
[1851]40! 1850 2016-04-08 13:29:27Z maronga
41! Module renamed
42! Adapted for modularization of microphysics
43!
[1832]44! 1831 2016-04-07 13:15:51Z hoffmann
45! turbulence renamed collision_turbulence
46!
[1683]47! 1682 2015-10-07 23:56:08Z knoop
48! Code annotations made doxygen readable
49!
[1375]50! 1374 2014-04-25 12:55:07Z raasch
51! missing variables added to ONLY list
52! rif removed from acc-present-list
53!
[1341]54! 1340 2014-03-25 19:45:13Z kanani
55! REAL constants defined as wp-kind
56!
[1321]57! 1320 2014-03-20 08:40:49Z raasch
[1320]58! ONLY-attribute added to USE-statements,
59! kind-parameters added to all INTEGER and REAL declaration statements,
60! kinds are defined in new module kinds,
61! revision history before 2012 removed,
62! comment fields (!:) to be used for variable explanations added to
63! all variable declaration statements
[98]64!
[1258]65! 1257 2013-11-08 15:18:40Z raasch
66! openacc loop and loop vector clauses removed
67!
[1182]68! 1179 2013-06-14 05:57:58Z raasch
69! use_reference renamed use_single_reference_value
70!
[1172]71! 1171 2013-05-30 11:27:45Z raasch
72! use_reference-case activated in accelerator version
73!
[1132]74! 1128 2013-04-12 06:19:32Z raasch
75! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
76! j_north
77!
[1066]78! 1065 2012-11-22 17:42:36Z hoffmann
79! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
80! effects of turbulence on autoconversion and accretion in two-moments cloud
81! physics scheme).
82!
[1037]83! 1036 2012-10-22 13:43:42Z raasch
84! code put under GPL (PALM 3.9)
85!
[1017]86! 1015 2012-09-27 09:23:24Z raasch
87! accelerator version (*_acc) added,
88! adjustment of mixing length to the Prandtl mixing length at first grid point
89! above ground removed
90!
[1011]91! 1010 2012-09-20 07:59:54Z raasch
92! cpp switch __nopointer added for pointer free version
93!
[1002]94! 1001 2012-09-13 14:08:46Z raasch
95! most arrays comunicated by module instead of parameter list
96!
[826]97! 825 2012-02-19 03:03:44Z raasch
98! wang_collision_kernel renamed wang_kernel
99!
[1]100! Revision 1.1  1997/09/19 07:40:24  raasch
101! Initial revision
102!
103!
104! Description:
105! ------------
[1682]106!> Diffusion- and dissipation terms for the TKE
[1]107!------------------------------------------------------------------------------!
[1682]108 MODULE diffusion_e_mod
109 
[1]110
111    PRIVATE
[2118]112    PUBLIC diffusion_e
[1]113   
114
115    INTERFACE diffusion_e
116       MODULE PROCEDURE diffusion_e
117       MODULE PROCEDURE diffusion_e_ij
118    END INTERFACE diffusion_e
119 
120 CONTAINS
121
122
123!------------------------------------------------------------------------------!
[1682]124! Description:
125! ------------
126!> Call for all grid points
[1]127!------------------------------------------------------------------------------!
[1001]128    SUBROUTINE diffusion_e( var, var_reference )
[1]129
[1320]130       USE arrays_3d,                                                          &
[2037]131           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
132                  drho_air, rho_air_zw
[1320]133           
134       USE control_parameters,                                                 &
[1831]135           ONLY:  atmos_ocean_sign, g, use_single_reference_value, &
[1320]136                  wall_adjustment, wall_adjustment_factor
[1831]137
[1320]138       USE grid_variables,                                                     &
139           ONLY:  ddx2, ddy2
140           
141       USE indices,                                                            &
[1374]142           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
143                  nzt
[1320]144           
145       USE kinds
[1849]146
147       USE microphysics_mod,                                                   &
148           ONLY:  collision_turbulence
149
[1320]150       USE particle_attributes,                                                &
151           ONLY:  use_sgs_for_particles, wang_kernel
[1]152
153       IMPLICIT NONE
154
[1682]155       INTEGER(iwp) ::  i              !<
156       INTEGER(iwp) ::  j              !<
157       INTEGER(iwp) ::  k              !<
158       REAL(wp)     ::  dvar_dz        !<
159       REAL(wp)     ::  l_stable       !<
160       REAL(wp)     ::  var_reference  !<
[1001]161
[1010]162#if defined( __nopointer )
[1682]163       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
[1010]164#else
[1682]165       REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
[1010]166#endif
[1682]167       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !<
168       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  l            !<
169       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  ll           !<
[1]170 
171
172!
[65]173!--    This if clause must be outside the k-loop because otherwise
174!--    runtime errors occur with -C hopt on NEC
[1179]175       IF ( use_single_reference_value )  THEN
[65]176
177          DO  i = nxl, nxr
178             DO  j = nys, nyn
179                DO  k = nzb_s_inner(j,i)+1, nzt
[1]180!
[65]181!--                Calculate the mixing length (for dissipation)
[97]182                   dvar_dz = atmos_ocean_sign * &
183                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]184                   IF ( dvar_dz > 0.0_wp ) THEN
185                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
186                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
[57]187                   ELSE
[65]188                      l_stable = l_grid(k)
[57]189                   ENDIF
[1]190!
[65]191!--                Adjustment of the mixing length
192                   IF ( wall_adjustment )  THEN
[94]193                      l(k,j)  = MIN( wall_adjustment_factor *          &
194                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
195                                     l_grid(k), l_stable )
196                      ll(k,j) = MIN( wall_adjustment_factor *          &
197                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
198                                     l_grid(k) )
[65]199                   ELSE
200                      l(k,j)  = MIN( l_grid(k), l_stable )
201                      ll(k,j) = l_grid(k)
202                   ENDIF
[1]203
[65]204                ENDDO
[1]205             ENDDO
[65]206
[1]207!
[65]208!--          Calculate the tendency terms
209             DO  j = nys, nyn
210                DO  k = nzb_s_inner(j,i)+1, nzt
[1]211
[1340]212                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
[65]213                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]214
[65]215                    tend(k,j,i) = tend(k,j,i)                                  &
[1]216                                        + (                                    &
217                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
218                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
219                                          ) * ddx2                             &
220                                        + (                                    &
221                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
222                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
223                                          ) * ddy2                             &
224                                        + (                                    &
225               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
[2037]226                                                             * rho_air_zw(k)   &
[1]227             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
[2037]228                                                             * rho_air_zw(k-1) &
229                                          ) * ddzw(k) * drho_air(k)            &
[1]230                             - dissipation(k,j)
231
[65]232                ENDDO
[1]233             ENDDO
[65]234
235!
236!--          Store dissipation if needed for calculating the sgs particle
237!--          velocities
[1065]238             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
[1831]239                  collision_turbulence )  THEN
[65]240                DO  j = nys, nyn
241                   DO  k = nzb_s_inner(j,i)+1, nzt
242                      diss(k,j,i) = dissipation(k,j)
243                   ENDDO
244                ENDDO
245             ENDIF
246
[1]247          ENDDO
248
[65]249       ELSE
250
251          DO  i = nxl, nxr
252             DO  j = nys, nyn
253                DO  k = nzb_s_inner(j,i)+1, nzt
254!
255!--                Calculate the mixing length (for dissipation)
[97]256                   dvar_dz = atmos_ocean_sign * &
257                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]258                   IF ( dvar_dz > 0.0_wp ) THEN
259                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
260                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
[65]261                   ELSE
262                      l_stable = l_grid(k)
263                   ENDIF
264!
265!--                Adjustment of the mixing length
266                   IF ( wall_adjustment )  THEN
[94]267                      l(k,j)  = MIN( wall_adjustment_factor *          &
268                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
269                                     l_grid(k), l_stable )
270                      ll(k,j) = MIN( wall_adjustment_factor *          &
271                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
272                                     l_grid(k) )
[65]273                   ELSE
274                      l(k,j)  = MIN( l_grid(k), l_stable )
275                      ll(k,j) = l_grid(k)
276                   ENDIF
277
278                ENDDO
279             ENDDO
280
281!
282!--          Calculate the tendency terms
[1]283             DO  j = nys, nyn
[19]284                DO  k = nzb_s_inner(j,i)+1, nzt
[65]285
[1340]286                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
287                                             e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[65]288
289                    tend(k,j,i) = tend(k,j,i)                                  &
290                                        + (                                    &
291                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
292                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
293                                          ) * ddx2                             &
294                                        + (                                    &
295                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
296                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
297                                          ) * ddy2                             &
298                                        + (                                    &
299               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
[2037]300                                                             * rho_air_zw(k)   &
[65]301             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
[2037]302                                                             * rho_air_zw(k-1) &
303                                          ) * ddzw(k) * drho_air(k)            &
[65]304                             - dissipation(k,j)
305
[1]306                ENDDO
307             ENDDO
308
[65]309!
310!--          Store dissipation if needed for calculating the sgs particle
311!--          velocities
[1065]312             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
[1831]313                  collision_turbulence )  THEN
[65]314                DO  j = nys, nyn
315                   DO  k = nzb_s_inner(j,i)+1, nzt
316                      diss(k,j,i) = dissipation(k,j)
317                   ENDDO
318                ENDDO
319             ENDIF
[1]320
[65]321          ENDDO
322
323       ENDIF
324
[1]325!
326!--    Boundary condition for dissipation
[1849]327       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  &
[1831]328            collision_turbulence )  THEN
[1]329          DO  i = nxl, nxr
330             DO  j = nys, nyn
331                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
332             ENDDO
333          ENDDO
334       ENDIF
335
336    END SUBROUTINE diffusion_e
337
338
339!------------------------------------------------------------------------------!
[1682]340! Description:
341! ------------
342!> Call for grid point i,j
[1]343!------------------------------------------------------------------------------!
[1001]344    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]345
[1320]346       USE arrays_3d,                                                          &
[2037]347           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
348                  drho_air, rho_air_zw
[1320]349         
350       USE control_parameters,                                                 &
[1831]351           ONLY:  atmos_ocean_sign, g, use_single_reference_value,             &
[1320]352                  wall_adjustment, wall_adjustment_factor
[1831]353
[1320]354       USE grid_variables,                                                     &
355           ONLY:  ddx2, ddy2
356           
357       USE indices,                                                            &
[1374]358           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
[1320]359           
360       USE kinds
[1849]361
362       USE microphysics_mod,                                                   &
363           ONLY:  collision_turbulence
364
[1320]365       USE particle_attributes,                                                &
366           ONLY:  use_sgs_for_particles, wang_kernel
[1]367
368       IMPLICIT NONE
369
[1682]370       INTEGER(iwp) ::  i              !<
371       INTEGER(iwp) ::  j              !<
372       INTEGER(iwp) ::  k              !<
373       REAL(wp)     ::  dvar_dz        !<
374       REAL(wp)     ::  l_stable       !<
375       REAL(wp)     ::  var_reference  !<
[1]376
[1010]377#if defined( __nopointer )
[1682]378       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
[1010]379#else
[1682]380       REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !<
[1010]381#endif
[1682]382       REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !<
383       REAL(wp), DIMENSION(nzb+1:nzt) ::  l            !<
384       REAL(wp), DIMENSION(nzb+1:nzt) ::  ll           !<
[1]385
[1001]386
[1]387!
388!--    Calculate the mixing length (for dissipation)
[19]389       DO  k = nzb_s_inner(j,i)+1, nzt
[97]390          dvar_dz = atmos_ocean_sign * &
391                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]392          IF ( dvar_dz > 0.0_wp ) THEN
[1179]393             IF ( use_single_reference_value )  THEN
[1340]394                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
395                                     SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
[57]396             ELSE
[1340]397                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
398                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
[57]399             ENDIF
[1]400          ELSE
401             l_stable = l_grid(k)
402          ENDIF
403!
404!--       Adjustment of the mixing length
405          IF ( wall_adjustment )  THEN
[94]406             l(k)  = MIN( wall_adjustment_factor *                     &
407                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
408                          l_stable )
409             ll(k) = MIN( wall_adjustment_factor *                     &
410                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]411          ELSE
412             l(k)  = MIN( l_grid(k), l_stable )
413             ll(k) = l_grid(k)
414          ENDIF
415!
416!--       Calculate the tendency term
[1340]417          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
418                                 SQRT( e(k,j,i) ) / l(k)
[1]419
420          tend(k,j,i) = tend(k,j,i)                                           &
421                                       + (                                    &
422                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
423                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
424                                         ) * ddx2                             &
425                                       + (                                    &
426                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
427                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
428                                         ) * ddy2                             &
429                                       + (                                    &
430              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
[2037]431                                                            * rho_air_zw(k)   &
[1]432            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
[2037]433                                                            * rho_air_zw(k-1) &
434                                         ) * ddzw(k) * drho_air(k)            &
[1]435                                       - dissipation(k)
436
437       ENDDO
438
439!
440!--    Store dissipation if needed for calculating the sgs particle velocities
[1849]441       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
442            collision_turbulence )  THEN
[19]443          DO  k = nzb_s_inner(j,i)+1, nzt
[1]444             diss(k,j,i) = dissipation(k)
445          ENDDO
446!
447!--       Boundary condition for dissipation
448          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
449       ENDIF
450
451    END SUBROUTINE diffusion_e_ij
452
453 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.