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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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