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

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

last commit documented

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