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

Last change on this file since 1841 was 1818, checked in by maronga, 8 years ago

last commit documented / copyright update

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