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

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

revised renaming of modules

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