source: palm/trunk/SOURCE/diffusion_u.f90 @ 1258

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

last commit documented

  • Property svn:keywords set to Id
File size: 24.3 KB
Line 
1 MODULE diffusion_u_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_u.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! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
48!
49! 667 2010-12-23 12:06:00Z suehring/gryschka
50! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
51!
52! 366 2009-08-25 08:06:27Z raasch
53! bc_ns replaced by bc_ns_cyc
54!
55! 106 2007-08-16 14:30:26Z raasch
56! Momentumflux at top (uswst) included as boundary condition,
57! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
58!
59! 75 2007-03-22 09:54:05Z raasch
60! Wall functions now include diabatic conditions, call of routine wall_fluxes,
61! z0 removed from argument list, uxrp eliminated
62!
63! 20 2007-02-26 00:12:32Z raasch
64! Bugfix: ddzw dimensioned 1:nzt"+1"
65!
66! RCS Log replace by Id keyword, revision history cleaned up
67!
68! Revision 1.15  2006/02/23 10:35:35  raasch
69! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
70! or nzb_diff_u, respectively, in vertical diffusion,
71! wall functions added for north and south walls, +z0 in argument list,
72! terms containing w(k-1,..) are removed from the Prandtl-layer equation
73! because they cause errors at the edges of topography
74! WARNING: loops containing the MAX function are still not properly vectorized!
75!
76! Revision 1.1  1997/09/12 06:23:51  raasch
77! Initial revision
78!
79!
80! Description:
81! ------------
82! Diffusion term of the u-component
83! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
84!        and slows down the speed on NEC about 5-10%
85!------------------------------------------------------------------------------!
86
87    USE wall_fluxes_mod
88
89    PRIVATE
90    PUBLIC diffusion_u, diffusion_u_acc
91
92    INTERFACE diffusion_u
93       MODULE PROCEDURE diffusion_u
94       MODULE PROCEDURE diffusion_u_ij
95    END INTERFACE diffusion_u
96
97    INTERFACE diffusion_u_acc
98       MODULE PROCEDURE diffusion_u_acc
99    END INTERFACE diffusion_u_acc
100
101 CONTAINS
102
103
104!------------------------------------------------------------------------------!
105! Call for all grid points
106!------------------------------------------------------------------------------!
107    SUBROUTINE diffusion_u
108
109       USE arrays_3d
110       USE control_parameters
111       USE grid_variables
112       USE indices
113
114       IMPLICIT NONE
115
116       INTEGER ::  i, j, k
117       REAL    ::  kmym, kmyp, kmzm, kmzp
118
119       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
120
121!
122!--    First calculate horizontal momentum flux u'v' at vertical walls,
123!--    if neccessary
124       IF ( topography /= 'flat' )  THEN
125          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
126                            nzb_u_outer, wall_u )
127       ENDIF
128
129       DO  i = nxlu, nxr
130          DO  j = nys, nyn
131!
132!--          Compute horizontal diffusion
133             DO  k = nzb_u_outer(j,i)+1, nzt
134!
135!--             Interpolate eddy diffusivities on staggered gridpoints
136                kmyp = 0.25 * &
137                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
138                kmym = 0.25 * &
139                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
140
141                tend(k,j,i) = tend(k,j,i)                                    &
142                      & + 2.0 * (                                            &
143                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
144                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
145                      &         ) * ddx2                                     &
146                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
147                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
148                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
149                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
150                      &   ) * ddy
151             ENDDO
152
153!
154!--          Wall functions at the north and south walls, respectively
155             IF ( wall_u(j,i) /= 0.0 )  THEN
156
157                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
158                   kmyp = 0.25 * &
159                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
160                   kmym = 0.25 * &
161                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
162
163                   tend(k,j,i) = tend(k,j,i)                                   &
164                                 + 2.0 * (                                     &
165                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
166                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
167                                         ) * ddx2                              &
168                                 + (   fyp(j,i) * (                            &
169                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
170                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
171                                                  )                            &
172                                     - fym(j,i) * (                            &
173                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
174                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
175                                                  )                            &
176                                     + wall_u(j,i) * usvs(k,j,i)               &
177                                   ) * ddy
178                ENDDO
179             ENDIF
180
181!
182!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
183!--          index k starts at nzb_u_inner+2.
184             DO  k = nzb_diff_u(j,i), nzt_diff
185!
186!--             Interpolate eddy diffusivities on staggered gridpoints
187                kmzp = 0.25 * &
188                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
189                kmzm = 0.25 * &
190                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
191
192                tend(k,j,i) = tend(k,j,i)                                    &
193                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
194                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
195                      &            )                                         &
196                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
197                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
198                      &            )                                          &
199                      &   ) * ddzw(k)
200             ENDDO
201
202!
203!--          Vertical diffusion at the first grid point above the surface,
204!--          if the momentum flux at the bottom is given by the Prandtl law or
205!--          if it is prescribed by the user.
206!--          Difference quotient of the momentum flux is not formed over half
207!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
208!--          with other (LES) modell showed that the values of the momentum
209!--          flux becomes too large in this case.
210!--          The term containing w(k-1,..) (see above equation) is removed here
211!--          because the vertical velocity is assumed to be zero at the surface.
212             IF ( use_surface_fluxes )  THEN
213                k = nzb_u_inner(j,i)+1
214!
215!--             Interpolate eddy diffusivities on staggered gridpoints
216                kmzp = 0.25 * &
217                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
218                kmzm = 0.25 * &
219                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
220
221                tend(k,j,i) = tend(k,j,i)                                    &
222                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
223                      &   ) * ddzw(k)                                        &
224                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
225                      &   + usws(j,i)                                        &
226                      &   ) * ddzw(k)
227             ENDIF
228
229!
230!--          Vertical diffusion at the first gridpoint below the top boundary,
231!--          if the momentum flux at the top is prescribed by the user
232             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
233                k = nzt
234!
235!--             Interpolate eddy diffusivities on staggered gridpoints
236                kmzp = 0.25 * &
237                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
238                kmzm = 0.25 * &
239                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
240
241                tend(k,j,i) = tend(k,j,i)                                    &
242                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
243                      &   ) * ddzw(k)                                        &
244                      & + ( -uswst(j,i)                                      &
245                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
246                      &   ) * ddzw(k)
247             ENDIF
248
249          ENDDO
250       ENDDO
251
252    END SUBROUTINE diffusion_u
253
254
255!------------------------------------------------------------------------------!
256! Call for all grid points - accelerator version
257!------------------------------------------------------------------------------!
258    SUBROUTINE diffusion_u_acc
259
260       USE arrays_3d
261       USE control_parameters
262       USE grid_variables
263       USE indices
264
265       IMPLICIT NONE
266
267       INTEGER ::  i, j, k
268       REAL    ::  kmym, kmyp, kmzm, kmzp
269
270       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
271       !$acc declare create ( usvs )
272
273!
274!--    First calculate horizontal momentum flux u'v' at vertical walls,
275!--    if neccessary
276       IF ( topography /= 'flat' )  THEN
277          CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
278                                nzb_u_outer, wall_u )
279       ENDIF
280
281       !$acc kernels present ( u, v, w, km, tend, usws, uswst )   &
282       !$acc         present ( ddzu, ddzw, fym, fyp, wall_u )           &
283       !$acc         present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
284       DO  i = i_left, i_right
285          DO  j = j_south, j_north
286!
287!--          Compute horizontal diffusion
288             DO  k = 1, nzt
289                IF ( k > nzb_u_outer(j,i) )  THEN
290!
291!--                Interpolate eddy diffusivities on staggered gridpoints
292                   kmyp = 0.25 * &
293                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
294                   kmym = 0.25 * &
295                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
296
297                   tend(k,j,i) = tend(k,j,i)                                   &
298                         & + 2.0 * (                                           &
299                         &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
300                         &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
301                         &         ) * ddx2                                    &
302                         & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
303                         &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
304                         &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
305                         &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
306                         &   ) * ddy
307                ENDIF
308             ENDDO
309
310!
311!--          Wall functions at the north and south walls, respectively
312             DO  k = 1, nzt
313                IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND. &
314                    wall_u(j,i) /= 0.0 )  THEN
315
316                   kmyp = 0.25 * &
317                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
318                   kmym = 0.25 * &
319                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
320
321                   tend(k,j,i) = tend(k,j,i)                                   &
322                                 + 2.0 * (                                     &
323                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
324                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
325                                         ) * ddx2                              &
326                                 + (   fyp(j,i) * (                            &
327                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
328                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
329                                                  )                            &
330                                     - fym(j,i) * (                            &
331                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
332                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
333                                                  )                            &
334                                     + wall_u(j,i) * usvs(k,j,i)               &
335                                   ) * ddy
336                ENDIF
337             ENDDO
338
339!
340!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
341!--          index k starts at nzb_u_inner+2.
342             DO  k = 1, nzt_diff
343                IF ( k >= nzb_diff_u(j,i) )  THEN
344!
345!--                Interpolate eddy diffusivities on staggered gridpoints
346                   kmzp = 0.25 * &
347                          ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
348                   kmzm = 0.25 * &
349                          ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
350
351                   tend(k,j,i) = tend(k,j,i)                                   &
352                         & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
353                         &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
354                         &            )                                        &
355                         &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
356                         &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
357                         &            )                                        &
358                         &   ) * ddzw(k)
359                ENDIF
360             ENDDO
361
362          ENDDO
363       ENDDO
364
365!
366!--    Vertical diffusion at the first grid point above the surface,
367!--    if the momentum flux at the bottom is given by the Prandtl law or
368!--    if it is prescribed by the user.
369!--    Difference quotient of the momentum flux is not formed over half
370!--    of the grid spacing (2.0*ddzw(k)) any more, since the comparison
371!--    with other (LES) modell showed that the values of the momentum
372!--    flux becomes too large in this case.
373!--    The term containing w(k-1,..) (see above equation) is removed here
374!--    because the vertical velocity is assumed to be zero at the surface.
375       IF ( use_surface_fluxes )  THEN
376
377          DO  i = i_left, i_right
378             DO  j = j_south, j_north
379         
380                k = nzb_u_inner(j,i)+1
381!
382!--             Interpolate eddy diffusivities on staggered gridpoints
383                kmzp = 0.25 * &
384                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
385                kmzm = 0.25 * &
386                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
387
388                tend(k,j,i) = tend(k,j,i)                                    &
389                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
390                      &   ) * ddzw(k)                                        &
391                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
392                      &   + usws(j,i)                                        &
393                      &   ) * ddzw(k)
394             ENDDO
395          ENDDO
396
397       ENDIF
398
399!
400!--    Vertical diffusion at the first gridpoint below the top boundary,
401!--    if the momentum flux at the top is prescribed by the user
402       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
403
404          k = nzt
405
406          DO  i = i_left, i_right
407             DO  j = j_south, j_north
408
409!
410!--             Interpolate eddy diffusivities on staggered gridpoints
411                kmzp = 0.25 * &
412                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
413                kmzm = 0.25 * &
414                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
415
416                tend(k,j,i) = tend(k,j,i)                                    &
417                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
418                      &   ) * ddzw(k)                                        &
419                      & + ( -uswst(j,i)                                      &
420                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
421                      &   ) * ddzw(k)
422             ENDDO
423          ENDDO
424
425       ENDIF
426       !$acc end kernels
427
428    END SUBROUTINE diffusion_u_acc
429
430
431!------------------------------------------------------------------------------!
432! Call for grid point i,j
433!------------------------------------------------------------------------------!
434    SUBROUTINE diffusion_u_ij( i, j )
435
436       USE arrays_3d
437       USE control_parameters
438       USE grid_variables
439       USE indices
440
441       IMPLICIT NONE
442
443       INTEGER ::  i, j, k
444       REAL    ::  kmym, kmyp, kmzm, kmzp
445
446       REAL, DIMENSION(nzb:nzt+1) ::  usvs
447
448!
449!--    Compute horizontal diffusion
450       DO  k = nzb_u_outer(j,i)+1, nzt
451!
452!--       Interpolate eddy diffusivities on staggered gridpoints
453          kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
454          kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
455
456          tend(k,j,i) = tend(k,j,i)                                          &
457                      & + 2.0 * (                                            &
458                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
459                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
460                      &         ) * ddx2                                     &
461                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
462                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
463                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
464                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
465                      &   ) * ddy
466       ENDDO
467
468!
469!--    Wall functions at the north and south walls, respectively
470       IF ( wall_u(j,i) .NE. 0.0 )  THEN
471
472!
473!--       Calculate the horizontal momentum flux u'v'
474          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
475                            usvs, 1.0, 0.0, 0.0, 0.0 )
476
477          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
478             kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
479             kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
480
481             tend(k,j,i) = tend(k,j,i)                                         &
482                                 + 2.0 * (                                     &
483                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
484                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
485                                         ) * ddx2                              &
486                                 + (   fyp(j,i) * (                            &
487                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
488                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
489                                                  )                            &
490                                     - fym(j,i) * (                            &
491                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
492                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
493                                                  )                            &
494                                     + wall_u(j,i) * usvs(k)                   &
495                                   ) * ddy
496          ENDDO
497       ENDIF
498
499!
500!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
501!--    index k starts at nzb_u_inner+2.
502       DO  k = nzb_diff_u(j,i), nzt_diff
503!
504!--       Interpolate eddy diffusivities on staggered gridpoints
505          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
506          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
507
508          tend(k,j,i) = tend(k,j,i)                                          &
509                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
510                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
511                      &            )                                         &
512                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
513                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
514                      &            )                                         &
515                      &   ) * ddzw(k)
516       ENDDO
517
518!
519!--    Vertical diffusion at the first grid point above the surface, if the
520!--    momentum flux at the bottom is given by the Prandtl law or if it is
521!--    prescribed by the user.
522!--    Difference quotient of the momentum flux is not formed over half of
523!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
524!--    other (LES) modell showed that the values of the momentum flux becomes
525!--    too large in this case.
526!--    The term containing w(k-1,..) (see above equation) is removed here
527!--    because the vertical velocity is assumed to be zero at the surface.
528       IF ( use_surface_fluxes )  THEN
529          k = nzb_u_inner(j,i)+1
530!
531!--       Interpolate eddy diffusivities on staggered gridpoints
532          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
533          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
534
535          tend(k,j,i) = tend(k,j,i)                                          &
536                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
537                      &   ) * ddzw(k)                                        &
538                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
539                      &   + usws(j,i)                                        &
540                      &   ) * ddzw(k)
541       ENDIF
542
543!
544!--    Vertical diffusion at the first gridpoint below the top boundary,
545!--    if the momentum flux at the top is prescribed by the user
546       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
547          k = nzt
548!
549!--       Interpolate eddy diffusivities on staggered gridpoints
550          kmzp = 0.25 * &
551                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
552          kmzm = 0.25 * &
553                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
554
555          tend(k,j,i) = tend(k,j,i)                                          &
556                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
557                      &   ) * ddzw(k)                                        &
558                      & + ( -uswst(j,i)                                      &
559                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
560                      &   ) * ddzw(k)
561       ENDIF
562
563    END SUBROUTINE diffusion_u_ij
564
565 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.