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

Last change on this file since 1037 was 1037, checked in by raasch, 9 years ago

last commit documented

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