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

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

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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