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

Last change on this file since 1289 was 1258, checked in by raasch, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 19.3 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 1258 2013-11-08 16:09:09Z raasch $
27!
28! 1257 2013-11-08 15:18:40Z raasch
29! openacc loop and loop vector clauses removed, declare create moved after
30! the FORTRAN declaration statement
31!
32! 1128 2013-04-12 06:19:32Z raasch
33! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
34! j_north
35!
36! 1036 2012-10-22 13:43:42Z raasch
37! code put under GPL (PALM 3.9)
38!
39! 1015 2012-09-27 09:23:24Z raasch
40! accelerator version (*_acc) added
41!
42! 1001 2012-09-13 14:08:46Z raasch
43! arrays comunicated by module instead of parameter list
44!
45! 978 2012-08-09 08:28:32Z fricke
46! outflow damping layer removed
47! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
48! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
49!
50! 667 2010-12-23 12:06:00Z suehring/gryschka
51! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
52!
53! 366 2009-08-25 08:06:27Z raasch
54! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
55!
56! 75 2007-03-22 09:54:05Z raasch
57! Wall functions now include diabatic conditions, call of routine wall_fluxes,
58! z0 removed from argument list
59!
60! 20 2007-02-26 00:12:32Z raasch
61! Bugfix: ddzw dimensioned 1:nzt"+1"
62!
63! RCS Log replace by Id keyword, revision history cleaned up
64!
65! Revision 1.12  2006/02/23 10:38:03  raasch
66! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
67! +z0 in argument list
68! WARNING: loops containing the MAX function are still not properly vectorized!
69!
70! Revision 1.1  1997/09/12 06:24:11  raasch
71! Initial revision
72!
73!
74! Description:
75! ------------
76! Diffusion term of the w-component
77!------------------------------------------------------------------------------!
78
79    USE wall_fluxes_mod
80
81    PRIVATE
82    PUBLIC diffusion_w, diffusion_w_acc
83
84    INTERFACE diffusion_w
85       MODULE PROCEDURE diffusion_w
86       MODULE PROCEDURE diffusion_w_ij
87    END INTERFACE diffusion_w
88
89    INTERFACE diffusion_w_acc
90       MODULE PROCEDURE diffusion_w_acc
91    END INTERFACE diffusion_w_acc
92
93 CONTAINS
94
95
96!------------------------------------------------------------------------------!
97! Call for all grid points
98!------------------------------------------------------------------------------!
99    SUBROUTINE diffusion_w
100
101       USE arrays_3d
102       USE control_parameters
103       USE grid_variables
104       USE indices
105
106       IMPLICIT NONE
107
108       INTEGER ::  i, j, k
109       REAL    ::  kmxm, kmxp, kmym, kmyp
110
111       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
112
113
114!
115!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
116!--    walls, if neccessary
117       IF ( topography /= 'flat' )  THEN
118          CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
119                            nzb_w_outer, wall_w_x )
120          CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
121                            nzb_w_outer, wall_w_y )
122       ENDIF
123
124       DO  i = nxl, nxr
125          DO  j = nys, nyn
126             DO  k = nzb_w_outer(j,i)+1, nzt-1
127!
128!--             Interpolate eddy diffusivities on staggered gridpoints
129                kmxp = 0.25 * &
130                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
131                kmxm = 0.25 * &
132                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
133                kmyp = 0.25 * &
134                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
135                kmym = 0.25 * &
136                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
137
138                tend(k,j,i) = tend(k,j,i)                                      &
139                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
140                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
141                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
142                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
143                      &   ) * ddx                                              &
144                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
145                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
146                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
147                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
148                      &   ) * ddy                                              &
149                      & + 2.0 * (                                              &
150                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
151                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
152                      &         ) * ddzu(k+1)
153             ENDDO
154
155!
156!--          Wall functions at all vertical walls, where necessary
157             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
158
159                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
160!
161!--                Interpolate eddy diffusivities on staggered gridpoints
162                   kmxp = 0.25 * &
163                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
164                   kmxm = 0.25 * &
165                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
166                   kmyp = 0.25 * &
167                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
168                   kmym = 0.25 * &
169                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
170
171                   tend(k,j,i) = tend(k,j,i)                                   &
172                                 + (   fwxp(j,i) * (                           &
173                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
174                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
175                                                   )                           &
176                                     - fwxm(j,i) * (                           &
177                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
178                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
179                                                   )                           &
180                                     + wall_w_x(j,i) * wsus(k,j,i)             &
181                                   ) * ddx                                     &
182                                 + (   fwyp(j,i) * (                           &
183                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
184                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
185                                                   )                           &
186                                     - fwym(j,i) * (                           &
187                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
188                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
189                                                   )                           &
190                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
191                                   ) * ddy                                     &
192                                 + 2.0 * (                                     &
193                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
194                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
195                                         ) * ddzu(k+1)
196                ENDDO
197             ENDIF
198
199          ENDDO
200       ENDDO
201
202    END SUBROUTINE diffusion_w
203
204
205!------------------------------------------------------------------------------!
206! Call for all grid points - accelerator version
207!------------------------------------------------------------------------------!
208    SUBROUTINE diffusion_w_acc
209
210       USE arrays_3d
211       USE control_parameters
212       USE grid_variables
213       USE indices
214
215       IMPLICIT NONE
216
217       INTEGER ::  i, j, k
218       REAL    ::  kmxm, kmxp, kmym, kmyp
219
220       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
221       !$acc declare create ( wsus, wsvs )
222
223!
224!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
225!--    walls, if neccessary
226       IF ( topography /= 'flat' )  THEN
227          CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
228                                nzb_w_outer, wall_w_x )
229          CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
230                                nzb_w_outer, wall_w_y )
231       ENDIF
232
233       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )    &
234       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
235       !$acc         present ( nzb_w_inner, nzb_w_outer )
236       DO  i = i_left, i_right
237          DO  j = j_south, j_north
238             DO  k = 1, nzt
239                IF ( k > nzb_w_outer(j,i) )  THEN
240!
241!--                Interpolate eddy diffusivities on staggered gridpoints
242                   kmxp = 0.25 * &
243                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
244                   kmxm = 0.25 * &
245                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
246                   kmyp = 0.25 * &
247                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
248                   kmym = 0.25 * &
249                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
250
251                   tend(k,j,i) = tend(k,j,i)                                     &
252                         & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx        &
253                         &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
254                         &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx          &
255                         &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)    &
256                         &   ) * ddx                                             &
257                         & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy        &
258                         &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
259                         &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy          &
260                         &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
261                         &   ) * ddy                                             &
262                         & + 2.0 * (                                             &
263                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
264                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
265                         &         ) * ddzu(k+1)
266                ENDIF
267             ENDDO
268
269!
270!--          Wall functions at all vertical walls, where necessary
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.