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

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

last commit documented

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