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

Last change on this file since 1079 was 1066, checked in by hoffmann, 12 years ago

last commit documented

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