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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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