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

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

last commit documented

  • Property svn:keywords set to Id
File size: 26.7 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-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 1874 2016-04-18 14:51:10Z suehring $
26!
27! 1873 2016-04-18 14:50:06Z maronga
28! Module renamed (removed _mod)
29!
30!
31! 1850 2016-04-08 13:29:27Z maronga
32! Module renamed
33!
34!
35! 1740 2016-01-13 08:19:40Z raasch
36! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
37!
38! 1682 2015-10-07 23:56:08Z knoop
39! Code annotations made doxygen readable
40!
41! 1340 2014-03-25 19:45:13Z kanani
42! REAL constants defined as wp-kind
43!
44! 1320 2014-03-20 08:40:49Z raasch
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
51!
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!
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!
60! 1036 2012-10-22 13:43:42Z raasch
61! code put under GPL (PALM 3.9)
62!
63! 1015 2012-09-27 09:23:24Z raasch
64! accelerator version (*_acc) added
65!
66! 1001 2012-09-13 14:08:46Z raasch
67! arrays comunicated by module instead of parameter list
68!
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!
73! Revision 1.1  1997/09/12 06:24:01  raasch
74! Initial revision
75!
76!
77! Description:
78! ------------
79!> Diffusion term of the v-component
80!------------------------------------------------------------------------------!
81 MODULE diffusion_v_mod
82 
83
84    USE wall_fluxes_mod
85
86    PRIVATE
87    PUBLIC diffusion_v, diffusion_v_acc
88
89    INTERFACE diffusion_v
90       MODULE PROCEDURE diffusion_v
91       MODULE PROCEDURE diffusion_v_ij
92    END INTERFACE diffusion_v
93
94    INTERFACE diffusion_v_acc
95       MODULE PROCEDURE diffusion_v_acc
96    END INTERFACE diffusion_v_acc
97
98 CONTAINS
99
100
101!------------------------------------------------------------------------------!
102! Description:
103! ------------
104!> Call for all grid points
105!------------------------------------------------------------------------------!
106    SUBROUTINE diffusion_v
107
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
123
124       IMPLICIT NONE
125
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  !<
133
134       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
135
136!
137!--    First calculate horizontal momentum flux v'u' at vertical walls,
138!--    if neccessary
139       IF ( topography /= 'flat' )  THEN
140          CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
141                            nzb_v_outer, wall_v )
142       ENDIF
143
144       DO  i = nxl, nxr
145          DO  j = nysv, nyn
146!
147!--          Compute horizontal diffusion
148             DO  k = nzb_v_outer(j,i)+1, nzt
149!
150!--             Interpolate eddy diffusivities on staggered gridpoints
151                kmxp = 0.25_wp * &
152                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
153                kmxm = 0.25_wp * &
154                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
155
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                                              &
162                      & + 2.0_wp * (                                           &
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) )      &
165                      &            ) * ddy2
166             ENDDO
167
168!
169!--          Wall functions at the left and right walls, respectively
170             IF ( wall_v(j,i) /= 0.0_wp )  THEN
171
172                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
173                   kmxp = 0.25_wp *                                            &
174                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
175                   kmxm = 0.25_wp *                                            &
176                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
177                   
178                   tend(k,j,i) = tend(k,j,i)                                   &
179                                 + 2.0_wp * (                                  &
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) ) &
182                                            ) * ddy2                           &
183                                 + (   fxp(j,i) * (                            &
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   &
186                                                  )                            &
187                                     - fxm(j,i) * (                            &
188                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
189                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
190                                                  )                            &
191                                     + wall_v(j,i) * vsus(k,j,i)               &
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.
199             DO  k = nzb_diff_v(j,i), nzt_diff
200!
201!--             Interpolate eddy diffusivities on staggered gridpoints
202                kmzp = 0.25_wp * &
203                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
204                kmzm = 0.25_wp * &
205                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
206
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                      &            )                                           &
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
223!--          comparison with other (LES) models showed that the values of
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
231                kmzp = 0.25_wp *                                               &
232                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
233
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)                                          &
239                      &   ) * ddzw(k)
240             ENDIF
241
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
245             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
246                k = nzt
247!
248!--             Interpolate eddy diffusivities on staggered gridpoints
249                kmzm = 0.25_wp *                                               &
250                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
251
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)    &
257                      &   ) * ddzw(k)
258             ENDIF
259
260          ENDDO
261       ENDDO
262
263    END SUBROUTINE diffusion_v
264
265
266!------------------------------------------------------------------------------!
267! Description:
268! ------------
269!> Call for all grid points - accelerator version
270!------------------------------------------------------------------------------!
271    SUBROUTINE diffusion_v_acc
272
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
288
289       IMPLICIT NONE
290
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  !<
298
299       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
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
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 )
308       ENDIF
309
310       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
311       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
312       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
313       DO  i = i_left, i_right
314          DO  j = j_south, j_north
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
321                   kmxp = 0.25_wp *                                            &
322                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
323                   kmxm = 0.25_wp *                                            &
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                                           &
332                         & + 2.0_wp * (                                        &
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) )   &
335                         &            ) * ddy2
336                ENDIF
337             ENDDO
338
339!
340!--          Wall functions at the left and right walls, respectively
341             DO  k = 1, nzt
342                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
343                    wall_v(j,i) /= 0.0_wp )  THEN
344
345                   kmxp = 0.25_wp *                                            &
346                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
347                   kmxm = 0.25_wp *                                            &
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)                                   &
351                                 + 2.0_wp * (                                  &
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) ) &
354                                            ) * ddy2                           &
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
375                   kmzp = 0.25_wp *                                            &
376                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
377                   kmzm = 0.25_wp *                                            &
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
400!--    comparison with other (LES) models showed that the values of
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
406          DO  i = i_left, i_right
407             DO  j = j_south, j_north
408         
409                k = nzb_v_inner(j,i)+1
410!
411!--             Interpolate eddy diffusivities on staggered gridpoints
412                kmzp = 0.25_wp *                                               &
413                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
414
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)                                          &
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
433          DO  i = i_left, i_right
434             DO  j = j_south, j_north
435
436!
437!--             Interpolate eddy diffusivities on staggered gridpoints
438                kmzm = 0.25_wp *                                               &
439                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
440
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)    &
446                      &   ) * ddzw(k)
447             ENDDO
448          ENDDO
449
450       ENDIF
451       !$acc end kernels
452
453    END SUBROUTINE diffusion_v_acc
454
455
456!------------------------------------------------------------------------------!
457! Description:
458! ------------
459!> Call for grid point i,j
460!------------------------------------------------------------------------------!
461    SUBROUTINE diffusion_v_ij( i, j )
462
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
476
477       IMPLICIT NONE
478
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  !<
486
487       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
488
489!
490!--    Compute horizontal diffusion
491       DO  k = nzb_v_outer(j,i)+1, nzt
492!
493!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
496
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                                              &
503                      & + 2.0_wp * (                                           &
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) )      &
506                      &            ) * ddy2
507       ENDDO
508
509!
510!--    Wall functions at the left and right walls, respectively
511       IF ( wall_v(j,i) /= 0.0_wp )  THEN
512
513!
514!--       Calculate the horizontal momentum flux v'u'
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 )
517
518          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
519             kmxp = 0.25_wp *                                                  &
520                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
521             kmxm = 0.25_wp *                                                  &
522                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
523
524             tend(k,j,i) = tend(k,j,i)                                         &
525                                 + 2.0_wp * (                                  &
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) ) &
528                                            ) * ddy2                           &
529                                 + (   fxp(j,i) * (                            &
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   &
532                                                  )                            &
533                                     - fxm(j,i) * (                            &
534                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
535                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
536                                                  )                            &
537                                     + wall_v(j,i) * vsus(k)                   &
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.
545       DO  k = nzb_diff_v(j,i), nzt_diff
546!
547!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
550
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                      &            )                                           &
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
567!--    other (LES) models showed that the values of the momentum flux becomes
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
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) )
576
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)                                          &
582                      &   ) * ddzw(k)
583       ENDIF
584
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
588       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
589          k = nzt
590!
591!--       Interpolate eddy diffusivities on staggered gridpoints
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.