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

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