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

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

last commit documented

  • Property svn:keywords set to Id
File size: 21.8 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!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_e.f90 1172 2013-05-30 11:46:00Z 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_reference )  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_reference )  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_reference )  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.