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

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

code has been put under the GNU General Public License (v3)

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