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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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