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

Last change on this file since 2371 was 2233, checked in by suehring, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 19.7 KB
Line 
1!> @file diffusion_e.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_e.f90 2233 2017-05-30 18:08:54Z kanani $
27!
28! 2232 2017-05-30 17:47:52Z suehring
29! Adjustments to new topography and surface concept
30!
31! 2118 2017-01-17 16:38:49Z raasch
32! OpenACC version of subroutine removed
33!
34! 2037 2016-10-26 11:15:40Z knoop
35! Anelastic approximation implemented
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1873 2016-04-18 14:50:06Z maronga
41! Module renamed (removed _mod)
42!
43! 1850 2016-04-08 13:29:27Z maronga
44! Module renamed
45! Adapted for modularization of microphysics
46!
47! 1831 2016-04-07 13:15:51Z hoffmann
48! turbulence renamed collision_turbulence
49!
50! 1682 2015-10-07 23:56:08Z knoop
51! Code annotations made doxygen readable
52!
53! 1374 2014-04-25 12:55:07Z raasch
54! missing variables added to ONLY list
55! rif removed from acc-present-list
56!
57! 1340 2014-03-25 19:45:13Z kanani
58! REAL constants defined as wp-kind
59!
60! 1320 2014-03-20 08:40:49Z raasch
61! ONLY-attribute added to USE-statements,
62! kind-parameters added to all INTEGER and REAL declaration statements,
63! kinds are defined in new module kinds,
64! revision history before 2012 removed,
65! comment fields (!:) to be used for variable explanations added to
66! all variable declaration statements
67!
68! 1257 2013-11-08 15:18:40Z raasch
69! openacc loop and loop vector clauses removed
70!
71! 1179 2013-06-14 05:57:58Z raasch
72! use_reference renamed use_single_reference_value
73!
74! 1171 2013-05-30 11:27:45Z raasch
75! use_reference-case activated in accelerator version
76!
77! 1128 2013-04-12 06:19:32Z raasch
78! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
79! j_north
80!
81! 1065 2012-11-22 17:42:36Z hoffmann
82! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
83! effects of turbulence on autoconversion and accretion in two-moments cloud
84! physics scheme).
85!
86! 1036 2012-10-22 13:43:42Z raasch
87! code put under GPL (PALM 3.9)
88!
89! 1015 2012-09-27 09:23:24Z raasch
90! accelerator version (*_acc) added,
91! adjustment of mixing length to the Prandtl mixing length at first grid point
92! above ground removed
93!
94! 1010 2012-09-20 07:59:54Z raasch
95! cpp switch __nopointer added for pointer free version
96!
97! 1001 2012-09-13 14:08:46Z raasch
98! most arrays comunicated by module instead of parameter list
99!
100! 825 2012-02-19 03:03:44Z raasch
101! wang_collision_kernel renamed wang_kernel
102!
103! Revision 1.1  1997/09/19 07:40:24  raasch
104! Initial revision
105!
106!
107! Description:
108! ------------
109!> Diffusion- and dissipation terms for the TKE
110!------------------------------------------------------------------------------!
111 MODULE diffusion_e_mod
112 
113
114    PRIVATE
115    PUBLIC diffusion_e
116   
117
118    INTERFACE diffusion_e
119       MODULE PROCEDURE diffusion_e
120       MODULE PROCEDURE diffusion_e_ij
121    END INTERFACE diffusion_e
122 
123 CONTAINS
124
125
126!------------------------------------------------------------------------------!
127! Description:
128! ------------
129!> Call for all grid points
130!------------------------------------------------------------------------------!
131    SUBROUTINE diffusion_e( var, var_reference )
132
133       USE arrays_3d,                                                          &
134           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
135                  drho_air, rho_air_zw
136           
137       USE control_parameters,                                                 &
138           ONLY:  atmos_ocean_sign, g, use_single_reference_value, &
139                  wall_adjustment, wall_adjustment_factor
140
141       USE grid_variables,                                                     &
142           ONLY:  ddx2, ddy2
143           
144       USE indices,                                                            &
145           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,    &
146                  nzt, wall_flags_0
147           
148       USE kinds
149
150       USE microphysics_mod,                                                   &
151           ONLY:  collision_turbulence
152
153       USE particle_attributes,                                                &
154           ONLY:  use_sgs_for_particles, wang_kernel
155
156       USE surface_mod,                                                        &
157          ONLY :  bc_h
158
159       IMPLICIT NONE
160
161       INTEGER(iwp) ::  i              !< running index x direction
162       INTEGER(iwp) ::  j              !< running index y direction
163       INTEGER(iwp) ::  k              !< running index z direction
164       INTEGER(iwp) ::  m              !< running index surface elements
165 
166       REAL(wp)     ::  dvar_dz        !<
167       REAL(wp)     ::  flag           !< flag to mask topography
168       REAL(wp)     ::  l_stable       !<
169       REAL(wp)     ::  var_reference  !<
170
171#if defined( __nopointer )
172       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
173#else
174       REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
175#endif
176       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !<
177       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  l            !<
178       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  ll           !<
179 
180
181!
182!--    This if clause must be outside the k-loop because otherwise
183!--    runtime errors occur with -C hopt on NEC
184       IF ( use_single_reference_value )  THEN
185
186          DO  i = nxl, nxr
187             DO  j = nys, nyn
188                DO  k = nzb+1, nzt
189!
190!--                Calculate the mixing length (for dissipation)
191                   dvar_dz = atmos_ocean_sign * &
192                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
193                   IF ( dvar_dz > 0.0_wp ) THEN
194                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
195                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
196                   ELSE
197                      l_stable = l_grid(k)
198                   ENDIF
199!
200!--                Adjustment of the mixing length
201                   IF ( wall_adjustment )  THEN
202                      l(k,j)  = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
203                                     l_grid(k), l_stable )
204                      ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
205                                     l_grid(k) )
206                   ELSE
207                      l(k,j)  = MIN( l_grid(k), l_stable )
208                      ll(k,j) = l_grid(k)
209                   ENDIF
210
211                ENDDO
212             ENDDO
213
214!
215!--          Calculate the tendency terms
216             DO  j = nys, nyn
217                DO  k = nzb+1, nzt
218!
219!--                Predetermine flag to mask topography
220                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
221
222                   dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
223                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
224
225                   tend(k,j,i) = tend(k,j,i)                                   &
226                                        + (                                    &
227                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
228                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
229                                          ) * ddx2  * flag                     &
230                                        + (                                    &
231                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
232                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
233                                          ) * ddy2  * flag                     &
234                                        + (                                    &
235               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
236                                                             * rho_air_zw(k)   &
237             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
238                                                             * rho_air_zw(k-1) &
239                                          ) * ddzw(k) * drho_air(k) * flag     &
240                             - dissipation(k,j) * flag
241
242                ENDDO
243             ENDDO
244
245!
246!--          Store dissipation if needed for calculating the sgs particle
247!--          velocities
248             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
249                  collision_turbulence )  THEN
250                DO  j = nys, nyn
251                   DO  k = nzb+1, nzt
252                      diss(k,j,i) = dissipation(k,j) *                         &
253                                        MERGE( 1.0_wp, 0.0_wp,                 &
254                                               BTEST( wall_flags_0(k,j,i), 0 ) )
255                   ENDDO
256                ENDDO
257             ENDIF
258
259          ENDDO
260
261       ELSE
262
263          DO  i = nxl, nxr
264             DO  j = nys, nyn
265                DO  k = nzb+1, nzt
266!
267!--                Calculate the mixing length (for dissipation)
268                   dvar_dz = atmos_ocean_sign * &
269                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
270                   IF ( dvar_dz > 0.0_wp ) THEN
271                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
272                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
273                   ELSE
274                      l_stable = l_grid(k)
275                   ENDIF
276!
277!--                Adjustment of the mixing length
278                   IF ( wall_adjustment )  THEN
279                      l(k,j)  = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
280                                     l_grid(k), l_stable )
281                      ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
282                                     l_grid(k) )
283                   ELSE
284                      l(k,j)  = MIN( l_grid(k), l_stable )
285                      ll(k,j) = l_grid(k)
286                   ENDIF
287
288                ENDDO
289             ENDDO
290
291!
292!--          Calculate the tendency terms
293             DO  j = nys, nyn
294                DO  k = nzb+1, nzt
295!
296!--                Predetermine flag to mask topography
297                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
298
299                   dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) )   *&
300                                             e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
301
302                   tend(k,j,i) = tend(k,j,i)                                   &
303                                        + (                                    &
304                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
305                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
306                                          ) * ddx2  * flag                     &
307                                        + (                                    &
308                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
309                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
310                                          ) * ddy2  * flag                     &
311                                        + (                                    &
312               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
313                                                             * rho_air_zw(k)   &
314             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
315                                                             * rho_air_zw(k-1) &
316                                          ) * ddzw(k) * drho_air(k) * flag     &
317                             - dissipation(k,j) * flag
318
319                ENDDO
320             ENDDO
321
322!
323!--          Store dissipation if needed for calculating the sgs particle
324!--          velocities
325             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
326                  collision_turbulence )  THEN
327                DO  j = nys, nyn
328                   DO  k = nzb+1, nzt
329                      diss(k,j,i) = dissipation(k,j) *                         &
330                                        MERGE( 1.0_wp, 0.0_wp,                 &
331                                               BTEST( wall_flags_0(k,j,i), 0 ) )
332                   ENDDO
333                ENDDO
334             ENDIF
335
336          ENDDO
337
338       ENDIF
339
340!
341!--    Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
342       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  &
343            collision_turbulence )  THEN
344!
345!--       Upward facing surfaces
346          DO  m = 1, bc_h(0)%ns
347             i = bc_h(0)%i(m)           
348             j = bc_h(0)%j(m)
349             k = bc_h(0)%k(m)
350             diss(k-1,j,i) = diss(k,j,i)
351          ENDDO
352!
353!--       Downward facing surfaces
354          DO  m = 1, bc_h(1)%ns
355             i = bc_h(1)%i(m)           
356             j = bc_h(1)%j(m)
357             k = bc_h(1)%k(m)
358             diss(k+1,j,i) = diss(k,j,i)
359          ENDDO
360
361       ENDIF
362
363    END SUBROUTINE diffusion_e
364
365
366!------------------------------------------------------------------------------!
367! Description:
368! ------------
369!> Call for grid point i,j
370!------------------------------------------------------------------------------!
371    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
372
373       USE arrays_3d,                                                          &
374           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
375                  drho_air, rho_air_zw
376         
377       USE control_parameters,                                                 &
378           ONLY:  atmos_ocean_sign, g, use_single_reference_value,             &
379                  wall_adjustment, wall_adjustment_factor
380
381       USE grid_variables,                                                     &
382           ONLY:  ddx2, ddy2
383           
384       USE indices,                                                            &
385           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
386           
387       USE kinds
388
389       USE microphysics_mod,                                                   &
390           ONLY:  collision_turbulence
391
392       USE particle_attributes,                                                &
393           ONLY:  use_sgs_for_particles, wang_kernel
394
395       USE surface_mod,                                                        &
396          ONLY :  bc_h
397
398       IMPLICIT NONE
399
400       INTEGER(iwp) ::  i              !< running index x direction
401       INTEGER(iwp) ::  j              !< running index y direction
402       INTEGER(iwp) ::  k              !< running index z direction
403       INTEGER(iwp) ::  m              !< running index surface elements
404       INTEGER(iwp) ::  surf_e         !< End index of surface elements at (j,i)-gridpoint
405       INTEGER(iwp) ::  surf_s         !< Start index of surface elements at (j,i)-gridpoint
406
407       REAL(wp)     ::  dvar_dz        !<
408       REAL(wp)     ::  flag           !< flag to mask topography
409       REAL(wp)     ::  l_stable       !<
410       REAL(wp)     ::  var_reference  !<
411
412#if defined( __nopointer )
413       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
414#else
415       REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !<
416#endif
417       REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !<
418       REAL(wp), DIMENSION(nzb+1:nzt) ::  l            !<
419       REAL(wp), DIMENSION(nzb+1:nzt) ::  ll           !<
420
421!
422!--    Calculate the mixing length (for dissipation)
423       DO  k = nzb+1, nzt
424!
425!--       Predetermine flag to mask topography
426          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
427
428          dvar_dz = atmos_ocean_sign * &
429                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
430          IF ( dvar_dz > 0.0_wp ) THEN
431             IF ( use_single_reference_value )  THEN
432                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
433                                     SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
434             ELSE
435                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
436                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
437             ENDIF
438          ELSE
439             l_stable = l_grid(k)
440          ENDIF
441!
442!--       Adjustment of the mixing length
443          IF ( wall_adjustment )  THEN
444             l(k)  = MIN( wall_adjustment_factor * l_wall(k,j,i),              &
445                          l_grid(k), l_stable )
446             ll(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
447          ELSE
448             l(k)  = MIN( l_grid(k), l_stable )
449             ll(k) = l_grid(k)
450          ENDIF
451!
452!--       Calculate the tendency term
453          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) *   &
454                                 SQRT( e(k,j,i) ) / l(k)
455
456          tend(k,j,i) = tend(k,j,i)                                            &
457                                       + (                                     &
458                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )   &
459                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )   &
460                                         ) * ddx2  * flag                      &
461                                       + (                                     &
462                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )   &
463                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )   &
464                                         ) * ddy2  * flag                      &
465                                       + (                                     &
466              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)  &
467                                                            * rho_air_zw(k)    &
468            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)    &
469                                                            * rho_air_zw(k-1)  &
470                                         ) * ddzw(k) * drho_air(k) * flag      &
471                                       - dissipation(k) * flag
472
473       ENDDO
474
475!
476!--    Store dissipation if needed for calculating the sgs particle velocities
477       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
478            collision_turbulence )  THEN
479          DO  k = nzb+1, nzt
480             diss(k,j,i) = dissipation(k) *                                    &
481                                        MERGE( 1.0_wp, 0.0_wp,                 &
482                                               BTEST( wall_flags_0(k,j,i), 0 ) )
483          ENDDO
484!
485!--       Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
486!--       For each surface type determine start and end index (in case of elevated
487!--       toopography several up/downward facing surfaces may exist.
488          surf_s = bc_h(0)%start_index(j,i)   
489          surf_e = bc_h(0)%end_index(j,i)   
490          DO  m = surf_s, surf_e
491             k             = bc_h(0)%k(m)
492             diss(k-1,j,i) = diss(k,j,i)
493          ENDDO
494!
495!--       Downward facing surfaces
496          surf_s = bc_h(1)%start_index(j,i)   
497          surf_e = bc_h(1)%end_index(j,i)   
498          DO  m = surf_s, surf_e
499             k             = bc_h(1)%k(m)
500             diss(k+1,j,i) = diss(k,j,i)
501          ENDDO
502       ENDIF
503
504    END SUBROUTINE diffusion_e_ij
505
506 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.