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

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

last commit documented, rc-file for example run updated

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