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

Last change on this file since 1179 was 1179, checked in by raasch, 8 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:


small bugfix for background communication (time_integration)

  • Property svn:keywords set to Id
File size: 21.9 KB
Line 
1 MODULE diffusion_e_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! use_reference renamed use_single_reference_value
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_e.f90 1179 2013-06-14 05:57:58Z raasch $
27!
28! 1171 2013-05-30 11:27:45Z raasch
29! use_reference-case activated in accelerator version
30!
31! 1128 2013-04-12 06:19:32Z raasch
32! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
33! j_north
34!
35! 1065 2012-11-22 17:42:36Z hoffmann
36! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
37! effects of turbulence on autoconversion and accretion in two-moments cloud
38! physics scheme).
39!
40! 1036 2012-10-22 13:43:42Z raasch
41! code put under GPL (PALM 3.9)
42!
43! 1015 2012-09-27 09:23:24Z raasch
44! accelerator version (*_acc) added,
45! adjustment of mixing length to the Prandtl mixing length at first grid point
46! above ground removed
47!
48! 1010 2012-09-20 07:59:54Z raasch
49! cpp switch __nopointer added for pointer free version
50!
51! 1001 2012-09-13 14:08:46Z raasch
52! most arrays comunicated by module instead of parameter list
53!
54! 825 2012-02-19 03:03:44Z raasch
55! wang_collision_kernel renamed wang_kernel
56!
57! 790 2011-11-29 03:11:20Z raasch
58! diss is also calculated in case that the Wang kernel is used
59!
60! 667 2010-12-23 12:06:00Z suehring/gryschka
61! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
62!
63! 97 2007-06-21 08:23:15Z raasch
64! Adjustment of mixing length calculation for the ocean version. zw added to
65! argument list.
66! This is also a bugfix, because the height above the topography is now
67! used instead of the height above level k=0.
68! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
69! use_pt_reference renamed use_reference
70!
71! 65 2007-03-13 12:11:43Z raasch
72! Reference temperature pt_reference can be used in buoyancy term
73!
74! 20 2007-02-26 00:12:32Z raasch
75! Bugfix: ddzw dimensioned 1:nzt"+1"
76! Calculation extended for gridpoint nzt
77!
78! RCS Log replace by Id keyword, revision history cleaned up
79!
80! Revision 1.18  2006/08/04 14:29:43  raasch
81! dissipation is stored in extra array diss if needed later on for calculating
82! the sgs particle velocities
83!
84! Revision 1.1  1997/09/19 07:40:24  raasch
85! Initial revision
86!
87!
88! Description:
89! ------------
90! Diffusion- and dissipation terms for the TKE
91!------------------------------------------------------------------------------!
92
93    PRIVATE
94    PUBLIC diffusion_e, diffusion_e_acc
95   
96
97    INTERFACE diffusion_e
98       MODULE PROCEDURE diffusion_e
99       MODULE PROCEDURE diffusion_e_ij
100    END INTERFACE diffusion_e
101 
102    INTERFACE diffusion_e_acc
103       MODULE PROCEDURE diffusion_e_acc
104    END INTERFACE diffusion_e_acc
105
106 CONTAINS
107
108
109!------------------------------------------------------------------------------!
110! Call for all grid points
111!------------------------------------------------------------------------------!
112    SUBROUTINE diffusion_e( var, var_reference )
113
114       USE arrays_3d
115       USE control_parameters
116       USE grid_variables
117       USE indices
118       USE particle_attributes
119
120       IMPLICIT NONE
121
122       INTEGER ::  i, j, k
123       REAL    ::  dvar_dz, l_stable, var_reference
124
125#if defined( __nopointer )
126       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
127#else
128       REAL, DIMENSION(:,:,:), POINTER ::  var
129#endif
130       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
131 
132
133!
134!--    This if clause must be outside the k-loop because otherwise
135!--    runtime errors occur with -C hopt on NEC
136       IF ( use_single_reference_value )  THEN
137
138          DO  i = nxl, nxr
139             DO  j = nys, nyn
140                DO  k = nzb_s_inner(j,i)+1, nzt
141!
142!--                Calculate the mixing length (for dissipation)
143                   dvar_dz = atmos_ocean_sign * &
144                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
145                   IF ( dvar_dz > 0.0 ) THEN
146                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
147                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
148                   ELSE
149                      l_stable = l_grid(k)
150                   ENDIF
151!
152!--                Adjustment of the mixing length
153                   IF ( wall_adjustment )  THEN
154                      l(k,j)  = MIN( wall_adjustment_factor *          &
155                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
156                                     l_grid(k), l_stable )
157                      ll(k,j) = MIN( wall_adjustment_factor *          &
158                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
159                                     l_grid(k) )
160                   ELSE
161                      l(k,j)  = MIN( l_grid(k), l_stable )
162                      ll(k,j) = l_grid(k)
163                   ENDIF
164
165                ENDDO
166             ENDDO
167
168!
169!--          Calculate the tendency terms
170             DO  j = nys, nyn
171                DO  k = nzb_s_inner(j,i)+1, nzt
172
173                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
174                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
175
176                    tend(k,j,i) = tend(k,j,i)                                  &
177                                        + (                                    &
178                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
179                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
180                                          ) * ddx2                             &
181                                        + (                                    &
182                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
183                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
184                                          ) * ddy2                             &
185                                        + (                                    &
186               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
187             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
188                                          ) * ddzw(k)                          &
189                             - dissipation(k,j)
190
191                ENDDO
192             ENDDO
193
194!
195!--          Store dissipation if needed for calculating the sgs particle
196!--          velocities
197             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
198                  turbulence )  THEN
199                DO  j = nys, nyn
200                   DO  k = nzb_s_inner(j,i)+1, nzt
201                      diss(k,j,i) = dissipation(k,j)
202                   ENDDO
203                ENDDO
204             ENDIF
205
206          ENDDO
207
208       ELSE
209
210          DO  i = nxl, nxr
211             DO  j = nys, nyn
212                DO  k = nzb_s_inner(j,i)+1, nzt
213!
214!--                Calculate the mixing length (for dissipation)
215                   dvar_dz = atmos_ocean_sign * &
216                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
217                   IF ( dvar_dz > 0.0 ) THEN
218                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
219                                        SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
220                   ELSE
221                      l_stable = l_grid(k)
222                   ENDIF
223!
224!--                Adjustment of the mixing length
225                   IF ( wall_adjustment )  THEN
226                      l(k,j)  = MIN( wall_adjustment_factor *          &
227                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
228                                     l_grid(k), l_stable )
229                      ll(k,j) = MIN( wall_adjustment_factor *          &
230                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
231                                     l_grid(k) )
232                   ELSE
233                      l(k,j)  = MIN( l_grid(k), l_stable )
234                      ll(k,j) = l_grid(k)
235                   ENDIF
236
237                ENDDO
238             ENDDO
239
240!
241!--          Calculate the tendency terms
242             DO  j = nys, nyn
243                DO  k = nzb_s_inner(j,i)+1, nzt
244
245                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
246                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
247
248                    tend(k,j,i) = tend(k,j,i)                                  &
249                                        + (                                    &
250                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
251                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
252                                          ) * ddx2                             &
253                                        + (                                    &
254                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
255                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
256                                          ) * ddy2                             &
257                                        + (                                    &
258               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
259             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
260                                          ) * ddzw(k)                          &
261                             - dissipation(k,j)
262
263                ENDDO
264             ENDDO
265
266!
267!--          Store dissipation if needed for calculating the sgs particle
268!--          velocities
269             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
270                  turbulence )  THEN
271                DO  j = nys, nyn
272                   DO  k = nzb_s_inner(j,i)+1, nzt
273                      diss(k,j,i) = dissipation(k,j)
274                   ENDDO
275                ENDDO
276             ENDIF
277
278          ENDDO
279
280       ENDIF
281
282!
283!--    Boundary condition for dissipation
284       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
285          DO  i = nxl, nxr
286             DO  j = nys, nyn
287                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
288             ENDDO
289          ENDDO
290       ENDIF
291
292    END SUBROUTINE diffusion_e
293
294
295!------------------------------------------------------------------------------!
296! Call for all grid points - accelerator version
297!------------------------------------------------------------------------------!
298    SUBROUTINE diffusion_e_acc( var, var_reference )
299
300       USE arrays_3d
301       USE control_parameters
302       USE grid_variables
303       USE indices
304       USE particle_attributes
305
306       IMPLICIT NONE
307
308       INTEGER ::  i, j, k
309       REAL    ::  dissipation, dvar_dz, l, ll, l_stable, var_reference
310
311#if defined( __nopointer )
312       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
313#else
314       REAL, DIMENSION(:,:,:), POINTER ::  var
315#endif
316
317
318!
319!--    This if clause must be outside the k-loop because otherwise
320!--    runtime errors occur with -C hopt on NEC
321       IF ( use_single_reference_value )  THEN
322
323          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
324          !$acc         present( nzb_s_inner, rif, tend, var, zu, zw )
325          !$acc loop
326          DO  i = i_left, i_right
327             DO  j = j_south, j_north
328                !$acc loop vector( 32 )
329                DO  k = 1, nzt
330
331                   IF ( k > nzb_s_inner(j,i) )  THEN
332!
333!--                   Calculate the mixing length (for dissipation)
334                      dvar_dz = atmos_ocean_sign * &
335                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
336                      IF ( dvar_dz > 0.0 ) THEN
337                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
338                                    SQRT( g / var_reference * dvar_dz ) + 1E-5
339                      ELSE
340                         l_stable = l_grid(k)
341                      ENDIF
342!
343!--                   Adjustment of the mixing length
344                      IF ( wall_adjustment )  THEN
345                         = MIN( wall_adjustment_factor *          &
346                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
347                                   l_grid(k), l_stable )
348                         ll = MIN( wall_adjustment_factor *          &
349                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
350                                   l_grid(k) )
351                      ELSE
352                         = MIN( l_grid(k), l_stable )
353                         ll = l_grid(k)
354                      ENDIF
355!
356!--                   Calculate the tendency terms
357                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
358                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
359
360                      tend(k,j,i) = tend(k,j,i)                                  &
361                                        + (                                    &
362                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
363                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
364                                          ) * ddx2                             &
365                                        + (                                    &
366                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
367                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
368                                          ) * ddy2                             &
369                                        + (                                    &
370               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
371             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
372                                          ) * ddzw(k)                          &
373                                  - dissipation
374
375!
376!--                   Store dissipation if needed for calculating the sgs particle
377!--                   velocities
378                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
379                           turbulence )  THEN
380                         diss(k,j,i) = dissipation
381                      ENDIF
382
383                   ENDIF
384
385                ENDDO
386             ENDDO
387          ENDDO
388          !$acc end kernels
389
390       ELSE
391
392          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
393          !$acc         present( nzb_s_inner, rif, tend, var, zu, zw )
394          !$acc loop
395          DO  i = i_left, i_right
396             DO  j = j_south, j_north
397                !$acc loop vector( 32 )
398                DO  k = 1, nzt
399
400                   IF ( k > nzb_s_inner(j,i) )  THEN
401!
402!--                   Calculate the mixing length (for dissipation)
403                      dvar_dz = atmos_ocean_sign * &
404                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
405                      IF ( dvar_dz > 0.0 ) THEN
406                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
407                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
408                      ELSE
409                         l_stable = l_grid(k)
410                      ENDIF
411!
412!--                   Adjustment of the mixing length
413                      IF ( wall_adjustment )  THEN
414                         = MIN( wall_adjustment_factor *          &
415                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
416                                     l_grid(k), l_stable )
417                         ll = MIN( wall_adjustment_factor *          &
418                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
419                                   l_grid(k) )
420                      ELSE
421                         = MIN( l_grid(k), l_stable )
422                         ll = l_grid(k)
423                      ENDIF
424!
425!--                   Calculate the tendency terms
426                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
427                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
428
429                      tend(k,j,i) = tend(k,j,i)                                &
430                                        + (                                    &
431                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
432                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
433                                          ) * ddx2                             &
434                                        + (                                    &
435                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
436                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
437                                          ) * ddy2                             &
438                                        + (                                    &
439               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
440             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
441                                          ) * ddzw(k)                          &
442                                  - dissipation
443
444!
445!--                   Store dissipation if needed for calculating the sgs
446!--                   particle  velocities
447                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
448                           turbulence )  THEN
449                         diss(k,j,i) = dissipation
450                      ENDIF
451
452                   ENDIF
453
454                ENDDO
455             ENDDO
456          ENDDO
457          !$acc end kernels
458
459       ENDIF
460
461!
462!--    Boundary condition for dissipation
463       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
464          !$acc kernels present( diss, nzb_s_inner )
465          !$acc loop
466          DO  i = i_left, i_right
467             !$acc loop vector( 32 )
468             DO  j = j_south, j_north
469                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
470             ENDDO
471          ENDDO
472          !$acc end kernels
473       ENDIF
474
475    END SUBROUTINE diffusion_e_acc
476
477
478!------------------------------------------------------------------------------!
479! Call for grid point i,j
480!------------------------------------------------------------------------------!
481    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
482
483       USE arrays_3d
484       USE control_parameters
485       USE grid_variables
486       USE indices
487       USE particle_attributes
488
489       IMPLICIT NONE
490
491       INTEGER ::  i, j, k
492       REAL    ::  dvar_dz, l_stable, var_reference
493
494#if defined( __nopointer )
495       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
496#else
497       REAL, DIMENSION(:,:,:), POINTER ::  var
498#endif
499       REAL, DIMENSION(nzb+1:nzt) ::  dissipation, l, ll
500
501
502!
503!--    Calculate the mixing length (for dissipation)
504       DO  k = nzb_s_inner(j,i)+1, nzt
505          dvar_dz = atmos_ocean_sign * &
506                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
507          IF ( dvar_dz > 0.0 ) THEN
508             IF ( use_single_reference_value )  THEN
509                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
510                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
511             ELSE
512                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
513                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
514             ENDIF
515          ELSE
516             l_stable = l_grid(k)
517          ENDIF
518!
519!--       Adjustment of the mixing length
520          IF ( wall_adjustment )  THEN
521             l(k)  = MIN( wall_adjustment_factor *                     &
522                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
523                          l_stable )
524             ll(k) = MIN( wall_adjustment_factor *                     &
525                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
526          ELSE
527             l(k)  = MIN( l_grid(k), l_stable )
528             ll(k) = l_grid(k)
529          ENDIF
530!
531!--       Calculate the tendency term
532          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
533                           SQRT( e(k,j,i) ) / l(k)
534
535          tend(k,j,i) = tend(k,j,i)                                           &
536                                       + (                                    &
537                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
538                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
539                                         ) * ddx2                             &
540                                       + (                                    &
541                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
542                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
543                                         ) * ddy2                             &
544                                       + (                                    &
545              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
546            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
547                                         ) * ddzw(k)                          &
548                                       - dissipation(k)
549
550       ENDDO
551
552!
553!--    Store dissipation if needed for calculating the sgs particle velocities
554       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
555          DO  k = nzb_s_inner(j,i)+1, nzt
556             diss(k,j,i) = dissipation(k)
557          ENDDO
558!
559!--       Boundary condition for dissipation
560          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
561       ENDIF
562
563    END SUBROUTINE diffusion_e_ij
564
565 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.