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

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

r1028 documented

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