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

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

last commit documented

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