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

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

added _mod string to several filenames to meet the naming convection for modules

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