source: palm/trunk/SOURCE/diffusion_u_mod.f90 @ 1851

Last change on this file since 1851 was 1851, checked in by maronga, 5 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 26.8 KB
Line 
1!> @file diffusion_u_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_u_mod.f90 1851 2016-04-08 13:32:50Z maronga $
26!
27! 1850 2016-04-08 13:29:27Z maronga
28! Module renamed
29!
30!
31! 1740 2016-01-13 08:19:40Z raasch
32! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
33!
34! 1691 2015-10-26 16:17:44Z maronga
35! Formatting corrections.
36!
37! 1682 2015-10-07 23:56:08Z knoop
38! Code annotations made doxygen readable
39!
40! 1340 2014-03-25 19:45:13Z kanani
41! REAL constants defined as wp-kind
42!
43! 1320 2014-03-20 08:40:49Z raasch
44! ONLY-attribute added to USE-statements,
45! kind-parameters added to all INTEGER and REAL declaration statements,
46! kinds are defined in new module kinds,
47! revision history before 2012 removed,
48! comment fields (!:) to be used for variable explanations added to
49! all variable declaration statements
50!
51! 1257 2013-11-08 15:18:40Z raasch
52! openacc loop and loop vector clauses removed, declare create moved after
53! the FORTRAN declaration statement
54!
55! 1128 2013-04-12 06:19:32Z raasch
56! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
57! j_north
58!
59! 1036 2012-10-22 13:43:42Z raasch
60! code put under GPL (PALM 3.9)
61!
62! 1015 2012-09-27 09:23:24Z raasch
63! accelerator version (*_acc) added
64!
65! 1001 2012-09-13 14:08:46Z raasch
66! arrays comunicated by module instead of parameter list
67!
68! 978 2012-08-09 08:28:32Z fricke
69! outflow damping layer removed
70! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
71!
72! Revision 1.1  1997/09/12 06:23:51  raasch
73! Initial revision
74!
75!
76! Description:
77! ------------
78!> Diffusion term of the u-component
79!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization
80!>       and slows down the speed on NEC about 5-10%
81!------------------------------------------------------------------------------!
82 MODULE diffusion_u_mod
83 
84
85    USE wall_fluxes_mod
86
87    PRIVATE
88    PUBLIC diffusion_u, diffusion_u_acc
89
90    INTERFACE diffusion_u
91       MODULE PROCEDURE diffusion_u
92       MODULE PROCEDURE diffusion_u_ij
93    END INTERFACE diffusion_u
94
95    INTERFACE diffusion_u_acc
96       MODULE PROCEDURE diffusion_u_acc
97    END INTERFACE diffusion_u_acc
98
99 CONTAINS
100
101
102!------------------------------------------------------------------------------!
103! Description:
104! ------------
105!> Call for all grid points
106!------------------------------------------------------------------------------!
107    SUBROUTINE diffusion_u
108
109       USE arrays_3d,                                                          &
110           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
111       
112       USE control_parameters,                                                 &
113           ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
114                  use_top_fluxes
115       
116       USE grid_variables,                                                     &
117           ONLY:  ddx, ddx2, ddy, fym, fyp, wall_u
118       
119       USE indices,                                                            &
120           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzb_diff_u, nzb_u_inner,      &
121                  nzb_u_outer, nzt, nzt_diff
122       
123       USE kinds
124
125       IMPLICIT NONE
126
127       INTEGER(iwp) ::  i     !<
128       INTEGER(iwp) ::  j     !<
129       INTEGER(iwp) ::  k     !<
130       REAL(wp)     ::  kmym  !<
131       REAL(wp)     ::  kmyp  !<
132       REAL(wp)     ::  kmzm  !<
133       REAL(wp)     ::  kmzp  !<
134
135       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs  !<
136
137!
138!--    First calculate horizontal momentum flux u'v' at vertical walls,
139!--    if neccessary
140       IF ( topography /= 'flat' )  THEN
141          CALL wall_fluxes( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, nzb_u_inner, &
142                            nzb_u_outer, wall_u )
143       ENDIF
144
145       DO  i = nxlu, nxr
146          DO  j = nys, nyn
147!
148!--          Compute horizontal diffusion
149             DO  k = nzb_u_outer(j,i)+1, nzt
150!
151!--             Interpolate eddy diffusivities on staggered gridpoints
152                kmyp = 0.25_wp *                                               &
153                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
154                kmym = 0.25_wp *                                               &
155                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
156
157                tend(k,j,i) = tend(k,j,i)                                      &
158                      & + 2.0_wp * (                                           &
159                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
160                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
161                      &            ) * ddx2                                    &
162                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
163                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
164                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
165                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
166                      &   ) * ddy
167             ENDDO
168
169!
170!--          Wall functions at the north and south walls, respectively
171             IF ( wall_u(j,i) /= 0.0_wp )  THEN
172
173                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
174                   kmyp = 0.25_wp *                                            &
175                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
176                   kmym = 0.25_wp *                                            &
177                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
178
179                   tend(k,j,i) = tend(k,j,i)                                   &
180                                 + 2.0_wp * (                                  &
181                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
182                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
183                                            ) * ddx2                           &
184                                 + (   fyp(j,i) * (                            &
185                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
186                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
187                                                  )                            &
188                                     - fym(j,i) * (                            &
189                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
190                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
191                                                  )                            &
192                                     + wall_u(j,i) * usvs(k,j,i)               &
193                                   ) * ddy
194                ENDDO
195             ENDIF
196
197!
198!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
199!--          index k starts at nzb_u_inner+2.
200             DO  k = nzb_diff_u(j,i), nzt_diff
201!
202!--             Interpolate eddy diffusivities on staggered gridpoints
203                kmzp = 0.25_wp *                                               &
204                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
205                kmzm = 0.25_wp *                                               &
206                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
207
208                tend(k,j,i) = tend(k,j,i)                                      &
209                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
210                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
211                      &            )                                           &
212                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
213                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
214                      &            )                                           &
215                      &   ) * ddzw(k)
216             ENDDO
217
218!
219!--          Vertical diffusion at the first grid point above the surface,
220!--          if the momentum flux at the bottom is given by the Prandtl law or
221!--          if it is prescribed by the user.
222!--          Difference quotient of the momentum flux is not formed over half
223!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
224!--          with other (LES) models showed that the values of the momentum
225!--          flux becomes too large in this case.
226!--          The term containing w(k-1,..) (see above equation) is removed here
227!--          because the vertical velocity is assumed to be zero at the surface.
228             IF ( use_surface_fluxes )  THEN
229                k = nzb_u_inner(j,i)+1
230!
231!--             Interpolate eddy diffusivities on staggered gridpoints
232                kmzp = 0.25_wp *                                               &
233                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
234
235                tend(k,j,i) = tend(k,j,i)                                      &
236                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
237                      &   ) * ddzw(k)                                          &
238                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
239                      &   + usws(j,i)                                          &
240                      &   ) * ddzw(k)
241             ENDIF
242
243!
244!--          Vertical diffusion at the first gridpoint below the top boundary,
245!--          if the momentum flux at the top is prescribed by the user
246             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
247                k = nzt
248!
249!--             Interpolate eddy diffusivities on staggered gridpoints
250                kmzm = 0.25_wp *                                               &
251                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
252
253                tend(k,j,i) = tend(k,j,i)                                      &
254                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
255                      &   ) * ddzw(k)                                          &
256                      & + ( -uswst(j,i)                                        &
257                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
258                      &   ) * ddzw(k)
259             ENDIF
260
261          ENDDO
262       ENDDO
263
264    END SUBROUTINE diffusion_u
265
266
267!------------------------------------------------------------------------------!
268! Description:
269! ------------
270!> Call for all grid points - accelerator version
271!------------------------------------------------------------------------------!
272    SUBROUTINE diffusion_u_acc
273
274       USE arrays_3d,                                                          &
275           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
276       
277       USE control_parameters,                                                 &
278           ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
279                  use_top_fluxes
280       
281       USE grid_variables,                                                     &
282           ONLY:  ddx, ddx2, ddy, fym, fyp, wall_u
283       
284       USE indices,                                                            &
285           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
286                  nzb_diff_u, nzb_u_inner, nzb_u_outer, nzt, nzt_diff
287       
288       USE kinds
289
290       IMPLICIT NONE
291
292       INTEGER(iwp) ::  i     !<
293       INTEGER(iwp) ::  j     !<
294       INTEGER(iwp) ::  k     !<
295       REAL(wp)     ::  kmym  !<
296       REAL(wp)     ::  kmyp  !<
297       REAL(wp)     ::  kmzm  !<
298       REAL(wp)     ::  kmzp  !<
299
300       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs  !<
301       !$acc declare create ( usvs )
302
303!
304!--    First calculate horizontal momentum flux u'v' at vertical walls,
305!--    if neccessary
306       IF ( topography /= 'flat' )  THEN
307          CALL wall_fluxes_acc( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
308                                nzb_u_inner, nzb_u_outer, wall_u )
309       ENDIF
310
311       !$acc kernels present ( u, v, w, km, tend, usws, uswst )                &
312       !$acc         present ( ddzu, ddzw, fym, fyp, wall_u )                  &
313       !$acc         present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
314       DO  i = i_left, i_right
315          DO  j = j_south, j_north
316!
317!--          Compute horizontal diffusion
318             DO  k = 1, nzt
319                IF ( k > nzb_u_outer(j,i) )  THEN
320!
321!--                Interpolate eddy diffusivities on staggered gridpoints
322                   kmyp = 0.25_wp *                                            &
323                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
324                   kmym = 0.25_wp *                                            &
325                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
326
327                   tend(k,j,i) = tend(k,j,i)                                   &
328                         & + 2.0_wp * (                                        &
329                         &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
330                         &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
331                         &            ) * ddx2                                 &
332                         & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
333                         &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
334                         &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
335                         &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
336                         &   ) * ddy
337                ENDIF
338             ENDDO
339
340!
341!--          Wall functions at the north and south walls, respectively
342             DO  k = 1, nzt
343                IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND.  &
344                    wall_u(j,i) /= 0.0_wp )  THEN
345
346                   kmyp = 0.25_wp *                                            &
347                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
348                   kmym = 0.25_wp *                                            &
349                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
350
351                   tend(k,j,i) = tend(k,j,i)                                   &
352                                 + 2.0_wp * (                                  &
353                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
354                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
355                                            ) * ddx2                           &
356                                 + (   fyp(j,i) * (                            &
357                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
358                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
359                                                  )                            &
360                                     - fym(j,i) * (                            &
361                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
362                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
363                                                  )                            &
364                                     + wall_u(j,i) * usvs(k,j,i)               &
365                                   ) * ddy
366                ENDIF
367             ENDDO
368
369!
370!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
371!--          index k starts at nzb_u_inner+2.
372             DO  k = 1, nzt_diff
373                IF ( k >= nzb_diff_u(j,i) )  THEN
374!
375!--                Interpolate eddy diffusivities on staggered gridpoints
376                   kmzp = 0.25_wp *                                            &
377                          ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
378                   kmzm = 0.25_wp *                                            &
379                          ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
380
381                   tend(k,j,i) = tend(k,j,i)                                   &
382                         & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
383                         &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
384                         &            )                                        &
385                         &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
386                         &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
387                         &            )                                        &
388                         &   ) * ddzw(k)
389                ENDIF
390             ENDDO
391
392          ENDDO
393       ENDDO
394
395!
396!--    Vertical diffusion at the first grid point above the surface,
397!--    if the momentum flux at the bottom is given by the Prandtl law or
398!--    if it is prescribed by the user.
399!--    Difference quotient of the momentum flux is not formed over half
400!--    of the grid spacing (2.0*ddzw(k)) any more, since the comparison
401!--    with other (LES) models showed that the values of the momentum
402!--    flux becomes too large in this case.
403!--    The term containing w(k-1,..) (see above equation) is removed here
404!--    because the vertical velocity is assumed to be zero at the surface.
405       IF ( use_surface_fluxes )  THEN
406
407          DO  i = i_left, i_right
408             DO  j = j_south, j_north
409         
410                k = nzb_u_inner(j,i)+1
411!
412!--             Interpolate eddy diffusivities on staggered gridpoints
413                kmzp = 0.25_wp *                                               &
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                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
418                      &   ) * ddzw(k)                                          &
419                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
420                      &   + usws(j,i)                                          &
421                      &   ) * ddzw(k)
422             ENDDO
423          ENDDO
424
425       ENDIF
426
427!
428!--    Vertical diffusion at the first gridpoint below the top boundary,
429!--    if the momentum flux at the top is prescribed by the user
430       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
431
432          k = nzt
433
434          DO  i = i_left, i_right
435             DO  j = j_south, j_north
436
437!
438!--             Interpolate eddy diffusivities on staggered gridpoints
439                kmzm = 0.25_wp *                                               &
440                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
441
442                tend(k,j,i) = tend(k,j,i)                                      &
443                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
444                      &   ) * ddzw(k)                                          &
445                      & + ( -uswst(j,i)                                        &
446                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
447                      &   ) * ddzw(k)
448             ENDDO
449          ENDDO
450
451       ENDIF
452       !$acc end kernels
453
454    END SUBROUTINE diffusion_u_acc
455
456
457!------------------------------------------------------------------------------!
458! Description:
459! ------------
460!> Call for grid point i,j
461!------------------------------------------------------------------------------!
462    SUBROUTINE diffusion_u_ij( i, j )
463
464       USE arrays_3d,                                                          &
465           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
466       
467       USE control_parameters,                                                 &
468           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
469       
470       USE grid_variables,                                                     &
471           ONLY:  ddx, ddx2, ddy, fym, fyp, wall_u
472       
473       USE indices,                                                            &
474           ONLY:  nzb, nzb_diff_u, nzb_u_inner, nzb_u_outer, nzt, nzt_diff
475       
476       USE kinds
477
478       IMPLICIT NONE
479
480       INTEGER(iwp) ::  i     !<
481       INTEGER(iwp) ::  j     !<
482       INTEGER(iwp) ::  k     !<
483       REAL(wp)     ::  kmym  !<
484       REAL(wp)     ::  kmyp  !<
485       REAL(wp)     ::  kmzm  !<
486       REAL(wp)     ::  kmzp  !<
487
488       REAL(wp), DIMENSION(nzb:nzt+1) ::  usvs  !<
489
490!
491!--    Compute horizontal diffusion
492       DO  k = nzb_u_outer(j,i)+1, nzt
493!
494!--       Interpolate eddy diffusivities on staggered gridpoints
495          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
496          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
497
498          tend(k,j,i) = tend(k,j,i)                                            &
499                      & + 2.0_wp * (                                           &
500                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
501                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
502                      &            ) * ddx2                                    &
503                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
504                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
505                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
506                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
507                      &   ) * ddy
508       ENDDO
509
510!
511!--    Wall functions at the north and south walls, respectively
512       IF ( wall_u(j,i) /= 0.0_wp )  THEN
513
514!
515!--       Calculate the horizontal momentum flux u'v'
516          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),        &
517                            usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
518
519          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
520             kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
521             kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
522
523             tend(k,j,i) = tend(k,j,i)                                         &
524                                 + 2.0_wp * (                                  &
525                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
526                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
527                                            ) * ddx2                           &
528                                 + (   fyp(j,i) * (                            &
529                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
530                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
531                                                  )                            &
532                                     - fym(j,i) * (                            &
533                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
534                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
535                                                  )                            &
536                                     + wall_u(j,i) * usvs(k)                   &
537                                   ) * ddy
538          ENDDO
539       ENDIF
540
541!
542!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
543!--    index k starts at nzb_u_inner+2.
544       DO  k = nzb_diff_u(j,i), nzt_diff
545!
546!--       Interpolate eddy diffusivities on staggered gridpoints
547          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
548          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
549
550          tend(k,j,i) = tend(k,j,i)                                            &
551                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
552                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
553                      &            )                                           &
554                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
555                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
556                      &            )                                           &
557                      &   ) * ddzw(k)
558       ENDDO
559
560!
561!--    Vertical diffusion at the first grid point above the surface, if the
562!--    momentum flux at the bottom is given by the Prandtl law or if it is
563!--    prescribed by the user.
564!--    Difference quotient of the momentum flux is not formed over half of
565!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
566!--    other (LES) models showed that the values of the momentum flux becomes
567!--    too large in this case.
568!--    The term containing w(k-1,..) (see above equation) is removed here
569!--    because the vertical velocity is assumed to be zero at the surface.
570       IF ( use_surface_fluxes )  THEN
571          k = nzb_u_inner(j,i)+1
572!
573!--       Interpolate eddy diffusivities on staggered gridpoints
574          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
575
576          tend(k,j,i) = tend(k,j,i)                                            &
577                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
578                      &   ) * ddzw(k)                                          &
579                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
580                      &   + usws(j,i)                                          &
581                      &   ) * ddzw(k)
582       ENDIF
583
584!
585!--    Vertical diffusion at the first gridpoint below the top boundary,
586!--    if the momentum flux at the top is prescribed by the user
587       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
588          k = nzt
589!
590!--       Interpolate eddy diffusivities on staggered gridpoints
591          kmzm = 0.25_wp *                                                     &
592                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
593
594          tend(k,j,i) = tend(k,j,i)                                            &
595                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
596                      &   ) * ddzw(k)                                          &
597                      & + ( -uswst(j,i)                                        &
598                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
599                      &   ) * ddzw(k)
600       ENDIF
601
602    END SUBROUTINE diffusion_u_ij
603
604 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.