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

Last change on this file since 1858 was 1851, checked in by maronga, 9 years ago

last commit documented

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