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

Last change on this file since 2091 was 2038, checked in by knoop, 7 years ago

last commit documented

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