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

Last change on this file since 2118 was 2118, checked in by raasch, 5 years ago

all OpenACC directives and related parts removed from the code

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