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

Last change on this file since 1015 was 1015, checked in by raasch, 12 years ago

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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