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

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

last commit documented

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