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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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