source: palm/trunk/SOURCE/diffusion_v.f90 @ 2000

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

  • Property svn:keywords set to Id
File size: 26.7 KB
RevLine 
[1873]1!> @file diffusion_v.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
4!
[2000]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.
[1036]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!
[1818]17! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[2000]22! Forced header and separation lines into 80 columns
[1341]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_v.f90 2000 2016-08-20 18:09:15Z knoop $
27!
[1874]28! 1873 2016-04-18 14:50:06Z maronga
29! Module renamed (removed _mod)
30!
31!
[1851]32! 1850 2016-04-08 13:29:27Z maronga
33! Module renamed
34!
35!
[1741]36! 1740 2016-01-13 08:19:40Z raasch
37! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
38!
[1683]39! 1682 2015-10-07 23:56:08Z knoop
40! Code annotations made doxygen readable
41!
[1341]42! 1340 2014-03-25 19:45:13Z kanani
43! REAL constants defined as wp-kind
44!
[1321]45! 1320 2014-03-20 08:40:49Z raasch
[1320]46! ONLY-attribute added to USE-statements,
47! kind-parameters added to all INTEGER and REAL declaration statements,
48! kinds are defined in new module kinds,
49! revision history before 2012 removed,
50! comment fields (!:) to be used for variable explanations added to
51! all variable declaration statements
[1321]52!
[1258]53! 1257 2013-11-08 15:18:40Z raasch
54! openacc loop and loop vector clauses removed, declare create moved after
55! the FORTRAN declaration statement
56!
[1132]57! 1128 2013-04-12 06:19:32Z raasch
58! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
59! j_north
60!
[1037]61! 1036 2012-10-22 13:43:42Z raasch
62! code put under GPL (PALM 3.9)
63!
[1017]64! 1015 2012-09-27 09:23:24Z raasch
65! accelerator version (*_acc) added
66!
[1002]67! 1001 2012-09-13 14:08:46Z raasch
68! arrays comunicated by module instead of parameter list
69!
[979]70! 978 2012-08-09 08:28:32Z fricke
71! outflow damping layer removed
72! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
73!
[1]74! Revision 1.1  1997/09/12 06:24:01  raasch
75! Initial revision
76!
77!
78! Description:
79! ------------
[1682]80!> Diffusion term of the v-component
[1]81!------------------------------------------------------------------------------!
[1682]82 MODULE diffusion_v_mod
83 
[1]84
[56]85    USE wall_fluxes_mod
86
[1]87    PRIVATE
[1015]88    PUBLIC diffusion_v, diffusion_v_acc
[1]89
90    INTERFACE diffusion_v
91       MODULE PROCEDURE diffusion_v
92       MODULE PROCEDURE diffusion_v_ij
93    END INTERFACE diffusion_v
94
[1015]95    INTERFACE diffusion_v_acc
96       MODULE PROCEDURE diffusion_v_acc
97    END INTERFACE diffusion_v_acc
98
[1]99 CONTAINS
100
101
102!------------------------------------------------------------------------------!
[1682]103! Description:
104! ------------
105!> Call for all grid points
[1]106!------------------------------------------------------------------------------!
[1001]107    SUBROUTINE diffusion_v
[1]108
[1320]109       USE arrays_3d,                                                          &
110           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, 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, ddy, ddy2, fxm, fxp, wall_v
118       
119       USE indices,                                                            &
120           ONLY:  nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner,      &
121                  nzb_v_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)     ::  kmxm  !<
131       REAL(wp)     ::  kmxp  !<
132       REAL(wp)     ::  kmzm  !<
133       REAL(wp)     ::  kmzp  !<
[1001]134
[1682]135       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
[1]136
[56]137!
138!--    First calculate horizontal momentum flux v'u' at vertical walls,
139!--    if neccessary
140       IF ( topography /= 'flat' )  THEN
[1320]141          CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
[56]142                            nzb_v_outer, wall_v )
143       ENDIF
144
[1]145       DO  i = nxl, nxr
[106]146          DO  j = nysv, nyn
[1]147!
148!--          Compute horizontal diffusion
149             DO  k = nzb_v_outer(j,i)+1, nzt
150!
151!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]152                kmxp = 0.25_wp * &
[978]153                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]154                kmxm = 0.25_wp * &
[978]155                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]156
[1320]157                tend(k,j,i) = tend(k,j,i)                                      &
158                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
159                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
160                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
161                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
162                      &   ) * ddx                                              &
[1340]163                      & + 2.0_wp * (                                           &
[1320]164                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
165                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
[1340]166                      &            ) * ddy2
[1]167             ENDDO
168
169!
170!--          Wall functions at the left and right walls, respectively
[1340]171             IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]172
[1]173                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]174                   kmxp = 0.25_wp *                                            &
[978]175                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]176                   kmxm = 0.25_wp *                                            &
[978]177                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
178                   
[1]179                   tend(k,j,i) = tend(k,j,i)                                   &
[1340]180                                 + 2.0_wp * (                                  &
[1]181                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
182                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]183                                            ) * ddy2                           &
[1]184                                 + (   fxp(j,i) * (                            &
[978]185                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
186                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
[1]187                                                  )                            &
188                                     - fxm(j,i) * (                            &
[978]189                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
190                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]191                                                  )                            &
[56]192                                     + wall_v(j,i) * vsus(k,j,i)               &
[1]193                                   ) * ddx
194                ENDDO
195             ENDIF
196
197!
198!--          Compute vertical diffusion. In case of simulating a Prandtl
199!--          layer, index k starts at nzb_v_inner+2.
[102]200             DO  k = nzb_diff_v(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-1,i)+km(k+1,j-1,i) )
[1340]205                kmzm = 0.25_wp * &
[1]206                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
207
[1320]208                tend(k,j,i) = tend(k,j,i)                                      &
209                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
210                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
211                      &            )                                           &
212                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
213                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
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
221!--          or if it is prescribed by the user.
222!--          Difference quotient of the momentum flux is not formed over
223!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]224!--          comparison with other (LES) models showed that the values of
[1]225!--          the momentum 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_v_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-1,i)+km(k+1,j-1,i) )
234
[1320]235                tend(k,j,i) = tend(k,j,i)                                      &
236                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
237                      &   ) * ddzw(k)                                          &
238                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
239                      &   + vsws(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-1,i)+km(k-1,j-1,i) )
252
[1320]253                tend(k,j,i) = tend(k,j,i)                                      &
254                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
255                      &   ) * ddzw(k)                                          &
256                      & + ( -vswst(j,i)                                        &
257                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
[102]258                      &   ) * ddzw(k)
259             ENDIF
260
[1]261          ENDDO
262       ENDDO
263
264    END SUBROUTINE diffusion_v
265
266
267!------------------------------------------------------------------------------!
[1682]268! Description:
269! ------------
270!> Call for all grid points - accelerator version
[1015]271!------------------------------------------------------------------------------!
272    SUBROUTINE diffusion_v_acc
273
[1320]274       USE arrays_3d,                                                          &
275           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, 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, ddy, ddy2, fxm, fxp, wall_v
283       
284       USE indices,                                                            &
285           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
286                  nzb_diff_v, nzb_v_inner, nzb_v_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)     ::  kmxm  !<
296       REAL(wp)     ::  kmxp  !<
297       REAL(wp)     ::  kmzm  !<
298       REAL(wp)     ::  kmzp  !<
[1015]299
[1682]300       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
[1015]301       !$acc declare create ( vsus )
302
303!
304!--    First calculate horizontal momentum flux v'u' at vertical walls,
305!--    if neccessary
306       IF ( topography /= 'flat' )  THEN
[1320]307          CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp,          &
308                                nzb_v_inner, nzb_v_outer, wall_v )
[1015]309       ENDIF
310
[1320]311       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
312       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
[1015]313       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
[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_v_outer(j,i) )  THEN
320!
321!--                Interpolate eddy diffusivities on staggered gridpoints
[1340]322                   kmxp = 0.25_wp *                                            &
[1015]323                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]324                   kmxm = 0.25_wp *                                            &
[1015]325                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
326
327                   tend(k,j,i) = tend(k,j,i)                                   &
328                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
329                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
330                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
331                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
332                         &   ) * ddx                                           &
[1340]333                         & + 2.0_wp * (                                        &
[1015]334                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
335                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
[1340]336                         &            ) * ddy2
[1015]337                ENDIF
338             ENDDO
339
340!
341!--          Wall functions at the left and right walls, respectively
342             DO  k = 1, nzt
[1320]343                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
[1340]344                    wall_v(j,i) /= 0.0_wp )  THEN
[1015]345
[1340]346                   kmxp = 0.25_wp *                                            &
[1015]347                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]348                   kmxm = 0.25_wp *                                            &
[1015]349                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+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)   * ( v(k,j+1,i) - v(k,j,i) ) &
354                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]355                                            ) * ddy2                           &
[1015]356                                 + (   fxp(j,i) * (                            &
357                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
358                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
359                                                  )                            &
360                                     - fxm(j,i) * (                            &
361                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
362                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
363                                                  )                            &
364                                     + wall_v(j,i) * vsus(k,j,i)               &
365                                   ) * ddx
366                ENDIF
367             ENDDO
368
369!
370!--          Compute vertical diffusion. In case of simulating a Prandtl
371!--          layer, index k starts at nzb_v_inner+2.
372             DO  k = 1, nzt_diff
373                IF ( k >= nzb_diff_v(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-1,i)+km(k+1,j-1,i) )
[1340]378                   kmzm = 0.25_wp *                                            &
[1015]379                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
380
381                   tend(k,j,i) = tend(k,j,i)                                   &
382                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
383                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
384                         &            )                                        &
385                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
386                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
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
398!--    or if it is prescribed by the user.
399!--    Difference quotient of the momentum flux is not formed over
400!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]401!--    comparison with other (LES) models showed that the values of
[1015]402!--    the momentum 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_v_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-1,i)+km(k+1,j-1,i) )
415
[1320]416                tend(k,j,i) = tend(k,j,i)                                      &
417                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
418                      &   ) * ddzw(k)                                          &
419                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
420                      &   + vsws(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-1,i)+km(k-1,j-1,i) )
441
[1320]442                tend(k,j,i) = tend(k,j,i)                                      &
443                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
444                      &   ) * ddzw(k)                                          &
445                      & + ( -vswst(j,i)                                        &
446                      &   - kmzm * ( v(k,j,i)   - v(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_v_acc
455
456
457!------------------------------------------------------------------------------!
[1682]458! Description:
459! ------------
460!> Call for grid point i,j
[1]461!------------------------------------------------------------------------------!
[1001]462    SUBROUTINE diffusion_v_ij( i, j )
[1]463
[1320]464       USE arrays_3d,                                                          &
465           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, 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, ddy, ddy2, fxm, fxp, wall_v
472       
473       USE indices,                                                            &
474           ONLY:  nzb, nzb_diff_v, nzb_v_inner, nzb_v_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)     ::  kmxm  !<
484       REAL(wp)     ::  kmxp  !<
485       REAL(wp)     ::  kmzm  !<
486       REAL(wp)     ::  kmzp  !<
[1]487
[1682]488       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
[1001]489
[1]490!
491!--    Compute horizontal diffusion
492       DO  k = nzb_v_outer(j,i)+1, nzt
493!
494!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]495          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
496          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]497
[1320]498          tend(k,j,i) = tend(k,j,i)                                            &
499                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
500                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
501                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
502                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
503                      &   ) * ddx                                              &
[1340]504                      & + 2.0_wp * (                                           &
[1320]505                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
506                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
[1340]507                      &            ) * ddy2
[1]508       ENDDO
509
510!
511!--    Wall functions at the left and right walls, respectively
[1340]512       IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]513
514!
515!--       Calculate the horizontal momentum flux v'u'
[1320]516          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i),        &
517                            vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[51]518
[1]519          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]520             kmxp = 0.25_wp *                                                  &
[978]521                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]522             kmxm = 0.25_wp *                                                  &
[978]523                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]524
525             tend(k,j,i) = tend(k,j,i)                                         &
[1340]526                                 + 2.0_wp * (                                  &
[1]527                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
528                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]529                                            ) * ddy2                           &
[1]530                                 + (   fxp(j,i) * (                            &
[978]531                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
532                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
[1]533                                                  )                            &
534                                     - fxm(j,i) * (                            &
[978]535                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
536                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]537                                                  )                            &
[51]538                                     + wall_v(j,i) * vsus(k)                   &
[1]539                                   ) * ddx
540          ENDDO
541       ENDIF
542
543!
544!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
545!--    index k starts at nzb_v_inner+2.
[102]546       DO  k = nzb_diff_v(j,i), nzt_diff
[1]547!
548!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]549          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
550          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]551
[1320]552          tend(k,j,i) = tend(k,j,i)                                            &
553                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
554                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
555                      &            )                                           &
556                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
557                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
558                      &            )                                           &
[1]559                      &   ) * ddzw(k)
560       ENDDO
561
562!
563!--    Vertical diffusion at the first grid point above the surface, if the
564!--    momentum flux at the bottom is given by the Prandtl law or if it is
565!--    prescribed by the user.
566!--    Difference quotient of the momentum flux is not formed over half of
567!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]568!--    other (LES) models showed that the values of the momentum flux becomes
[1]569!--    too large in this case.
570!--    The term containing w(k-1,..) (see above equation) is removed here
571!--    because the vertical velocity is assumed to be zero at the surface.
572       IF ( use_surface_fluxes )  THEN
573          k = nzb_v_inner(j,i)+1
574!
575!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]576          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]577
[1320]578          tend(k,j,i) = tend(k,j,i)                                            &
579                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
580                      &   ) * ddzw(k)                                          &
581                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
582                      &   + vsws(j,i)                                          &
[1]583                      &   ) * ddzw(k)
584       ENDIF
585
[102]586!
587!--    Vertical diffusion at the first gridpoint below the top boundary,
588!--    if the momentum flux at the top is prescribed by the user
[103]589       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[102]590          k = nzt
591!
592!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]593          kmzm = 0.25_wp * &
[102]594                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
595
[1320]596          tend(k,j,i) = tend(k,j,i)                                            &
597                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
598                      &   ) * ddzw(k)                                          &
599                      & + ( -vswst(j,i)                                        &
600                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
[102]601                      &   ) * ddzw(k)
602       ENDIF
603
[1]604    END SUBROUTINE diffusion_v_ij
605
[1321]606 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.