source: palm/trunk/SOURCE/diffusion_w.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 8 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:


openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:


bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

  • Property svn:keywords set to Id
File size: 19.2 KB
Line 
1 MODULE diffusion_w_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! openacc loop and loop vector clauses removed, declare create moved after
23! the FORTRAN declaration statement
24!
25! Former revisions:
26! -----------------
27! $Id: diffusion_w.f90 1257 2013-11-08 15:18:40Z raasch $
28!
29! 1128 2013-04-12 06:19:32Z raasch
30! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
31! j_north
32!
33! 1036 2012-10-22 13:43:42Z raasch
34! code put under GPL (PALM 3.9)
35!
36! 1015 2012-09-27 09:23:24Z raasch
37! accelerator version (*_acc) added
38!
39! 1001 2012-09-13 14:08:46Z raasch
40! arrays comunicated by module instead of parameter list
41!
42! 978 2012-08-09 08:28:32Z fricke
43! outflow damping layer removed
44! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
45! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
46!
47! 667 2010-12-23 12:06:00Z suehring/gryschka
48! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
49!
50! 366 2009-08-25 08:06:27Z raasch
51! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
52!
53! 75 2007-03-22 09:54:05Z raasch
54! Wall functions now include diabatic conditions, call of routine wall_fluxes,
55! z0 removed from argument list
56!
57! 20 2007-02-26 00:12:32Z raasch
58! Bugfix: ddzw dimensioned 1:nzt"+1"
59!
60! RCS Log replace by Id keyword, revision history cleaned up
61!
62! Revision 1.12  2006/02/23 10:38:03  raasch
63! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
64! +z0 in argument list
65! WARNING: loops containing the MAX function are still not properly vectorized!
66!
67! Revision 1.1  1997/09/12 06:24:11  raasch
68! Initial revision
69!
70!
71! Description:
72! ------------
73! Diffusion term of the w-component
74!------------------------------------------------------------------------------!
75
76    USE wall_fluxes_mod
77
78    PRIVATE
79    PUBLIC diffusion_w, diffusion_w_acc
80
81    INTERFACE diffusion_w
82       MODULE PROCEDURE diffusion_w
83       MODULE PROCEDURE diffusion_w_ij
84    END INTERFACE diffusion_w
85
86    INTERFACE diffusion_w_acc
87       MODULE PROCEDURE diffusion_w_acc
88    END INTERFACE diffusion_w_acc
89
90 CONTAINS
91
92
93!------------------------------------------------------------------------------!
94! Call for all grid points
95!------------------------------------------------------------------------------!
96    SUBROUTINE diffusion_w
97
98       USE arrays_3d
99       USE control_parameters
100       USE grid_variables
101       USE indices
102
103       IMPLICIT NONE
104
105       INTEGER ::  i, j, k
106       REAL    ::  kmxm, kmxp, kmym, kmyp
107
108       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
109
110
111!
112!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
113!--    walls, if neccessary
114       IF ( topography /= 'flat' )  THEN
115          CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
116                            nzb_w_outer, wall_w_x )
117          CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
118                            nzb_w_outer, wall_w_y )
119       ENDIF
120
121       DO  i = nxl, nxr
122          DO  j = nys, nyn
123             DO  k = nzb_w_outer(j,i)+1, nzt-1
124!
125!--             Interpolate eddy diffusivities on staggered gridpoints
126                kmxp = 0.25 * &
127                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
128                kmxm = 0.25 * &
129                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
130                kmyp = 0.25 * &
131                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
132                kmym = 0.25 * &
133                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
134
135                tend(k,j,i) = tend(k,j,i)                                      &
136                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
137                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
138                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
139                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
140                      &   ) * ddx                                              &
141                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
142                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
143                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
144                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
145                      &   ) * ddy                                              &
146                      & + 2.0 * (                                              &
147                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
148                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
149                      &         ) * ddzu(k+1)
150             ENDDO
151
152!
153!--          Wall functions at all vertical walls, where necessary
154             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
155
156                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
157!
158!--                Interpolate eddy diffusivities on staggered gridpoints
159                   kmxp = 0.25 * &
160                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
161                   kmxm = 0.25 * &
162                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
163                   kmyp = 0.25 * &
164                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
165                   kmym = 0.25 * &
166                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
167
168                   tend(k,j,i) = tend(k,j,i)                                   &
169                                 + (   fwxp(j,i) * (                           &
170                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
171                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
172                                                   )                           &
173                                     - fwxm(j,i) * (                           &
174                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
175                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
176                                                   )                           &
177                                     + wall_w_x(j,i) * wsus(k,j,i)             &
178                                   ) * ddx                                     &
179                                 + (   fwyp(j,i) * (                           &
180                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
181                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
182                                                   )                           &
183                                     - fwym(j,i) * (                           &
184                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
185                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
186                                                   )                           &
187                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
188                                   ) * ddy                                     &
189                                 + 2.0 * (                                     &
190                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
191                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
192                                         ) * ddzu(k+1)
193                ENDDO
194             ENDIF
195
196          ENDDO
197       ENDDO
198
199    END SUBROUTINE diffusion_w
200
201
202!------------------------------------------------------------------------------!
203! Call for all grid points - accelerator version
204!------------------------------------------------------------------------------!
205    SUBROUTINE diffusion_w_acc
206
207       USE arrays_3d
208       USE control_parameters
209       USE grid_variables
210       USE indices
211
212       IMPLICIT NONE
213
214       INTEGER ::  i, j, k
215       REAL    ::  kmxm, kmxp, kmym, kmyp
216
217       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
218       !$acc declare create ( wsus, wsvs )
219
220!
221!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
222!--    walls, if neccessary
223       IF ( topography /= 'flat' )  THEN
224          CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
225                                nzb_w_outer, wall_w_x )
226          CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
227                                nzb_w_outer, wall_w_y )
228       ENDIF
229
230       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )    &
231       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
232       !$acc         present ( nzb_w_inner, nzb_w_outer )
233       DO  i = i_left, i_right
234          DO  j = j_south, j_north
235             DO  k = 1, nzt
236                IF ( k > nzb_w_outer(j,i) )  THEN
237!
238!--                Interpolate eddy diffusivities on staggered gridpoints
239                   kmxp = 0.25 * &
240                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
241                   kmxm = 0.25 * &
242                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
243                   kmyp = 0.25 * &
244                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
245                   kmym = 0.25 * &
246                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
247
248                   tend(k,j,i) = tend(k,j,i)                                     &
249                         & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx        &
250                         &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
251                         &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx          &
252                         &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)    &
253                         &   ) * ddx                                             &
254                         & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy        &
255                         &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
256                         &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy          &
257                         &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
258                         &   ) * ddy                                             &
259                         & + 2.0 * (                                             &
260                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
261                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
262                         &         ) * ddzu(k+1)
263                ENDIF
264             ENDDO
265
266!
267!--          Wall functions at all vertical walls, where necessary
268             DO  k = 1,nzt
269
270                IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
271                     wall_w_x(j,i) /= 0.0  .AND.  wall_w_y(j,i) /= 0.0 )  THEN
272!
273!--                Interpolate eddy diffusivities on staggered gridpoints
274                   kmxp = 0.25 * &
275                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
276                   kmxm = 0.25 * &
277                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
278                   kmyp = 0.25 * &
279                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
280                   kmym = 0.25 * &
281                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
282
283                   tend(k,j,i) = tend(k,j,i)                                   &
284                                 + (   fwxp(j,i) * (                           &
285                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
286                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
287                                                   )                           &
288                                     - fwxm(j,i) * (                           &
289                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
290                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
291                                                   )                           &
292                                     + wall_w_x(j,i) * wsus(k,j,i)             &
293                                   ) * ddx                                     &
294                                 + (   fwyp(j,i) * (                           &
295                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
296                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
297                                                   )                           &
298                                     - fwym(j,i) * (                           &
299                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
300                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
301                                                   )                           &
302                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
303                                   ) * ddy                                     &
304                                 + 2.0 * (                                     &
305                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
306                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
307                                         ) * ddzu(k+1)
308                ENDIF
309             ENDDO
310
311          ENDDO
312       ENDDO
313       !$acc end kernels
314
315    END SUBROUTINE diffusion_w_acc
316
317
318!------------------------------------------------------------------------------!
319! Call for grid point i,j
320!------------------------------------------------------------------------------!
321    SUBROUTINE diffusion_w_ij( i, j )
322
323       USE arrays_3d
324       USE control_parameters
325       USE grid_variables
326       USE indices
327
328       IMPLICIT NONE
329
330       INTEGER ::  i, j, k
331       REAL    ::  kmxm, kmxp, kmym, kmyp
332
333       REAL, DIMENSION(nzb:nzt+1) ::  wsus, wsvs
334
335
336       DO  k = nzb_w_outer(j,i)+1, nzt-1
337!
338!--       Interpolate eddy diffusivities on staggered gridpoints
339          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
340          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
341          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
342          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
343
344          tend(k,j,i) = tend(k,j,i)                                            &
345                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
346                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
347                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
348                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
349                      &   ) * ddx                                              &
350                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
351                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
352                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
353                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
354                      &   ) * ddy                                              &
355                      & + 2.0 * (                                              &
356                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
357                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
358                      &         ) * ddzu(k+1)
359       ENDDO
360
361!
362!--    Wall functions at all vertical walls, where necessary
363       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
364
365!
366!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
367          IF ( wall_w_x(j,i) /= 0.0 )  THEN
368             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
369                               wsus, 0.0, 0.0, 0.0, 1.0 )
370          ELSE
371             wsus = 0.0
372          ENDIF
373
374          IF ( wall_w_y(j,i) /= 0.0 )  THEN
375             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),  &
376                               wsvs, 0.0, 0.0, 1.0, 0.0 )
377          ELSE
378             wsvs = 0.0
379          ENDIF
380
381          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
382!
383!--          Interpolate eddy diffusivities on staggered gridpoints
384             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
385             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
386             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
387             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
388
389             tend(k,j,i) = tend(k,j,i)                                         &
390                                 + (   fwxp(j,i) * (                           &
391                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
392                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
393                                                   )                           &
394                                     - fwxm(j,i) * (                           &
395                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
396                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
397                                                   )                           &
398                                     + wall_w_x(j,i) * wsus(k)                 &
399                                   ) * ddx                                     &
400                                 + (   fwyp(j,i) * (                           &
401                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
402                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
403                                                   )                           &
404                                     - fwym(j,i) * (                           &
405                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
406                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
407                                                   )                           &
408                                     + wall_w_y(j,i) * wsvs(k)                 &
409                                   ) * ddy                                     &
410                                 + 2.0 * (                                     &
411                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
412                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
413                                         ) * ddzu(k+1)
414          ENDDO
415       ENDIF
416
417    END SUBROUTINE diffusion_w_ij
418
419 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.