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

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

revised renaming of modules

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