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

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