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

Last change on this file since 1740 was 1740, checked in by raasch, 6 years ago

unnecessary calculations of kmzm and kmzp removed from wall bounded parts

  • Property svn:keywords set to Id
File size: 26.5 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! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 1740 2016-01-13 08:19:40Z raasch $
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
223                tend(k,j,i) = tend(k,j,i)                                      &
224                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
225                      &   ) * ddzw(k)                                          &
226                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
227                      &   + vsws(j,i)                                          &
228                      &   ) * ddzw(k)
229             ENDIF
230
231!
232!--          Vertical diffusion at the first gridpoint below the top boundary,
233!--          if the momentum flux at the top is prescribed by the user
234             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
235                k = nzt
236!
237!--             Interpolate eddy diffusivities on staggered gridpoints
238                kmzm = 0.25_wp *                                               &
239                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
240
241                tend(k,j,i) = tend(k,j,i)                                      &
242                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
243                      &   ) * ddzw(k)                                          &
244                      & + ( -vswst(j,i)                                        &
245                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
246                      &   ) * ddzw(k)
247             ENDIF
248
249          ENDDO
250       ENDDO
251
252    END SUBROUTINE diffusion_v
253
254
255!------------------------------------------------------------------------------!
256! Description:
257! ------------
258!> Call for all grid points - accelerator version
259!------------------------------------------------------------------------------!
260    SUBROUTINE diffusion_v_acc
261
262       USE arrays_3d,                                                          &
263           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
264       
265       USE control_parameters,                                                 &
266           ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
267                  use_top_fluxes
268       
269       USE grid_variables,                                                     &
270           ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
271       
272       USE indices,                                                            &
273           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
274                  nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
275       
276       USE kinds
277
278       IMPLICIT NONE
279
280       INTEGER(iwp) ::  i     !<
281       INTEGER(iwp) ::  j     !<
282       INTEGER(iwp) ::  k     !<
283       REAL(wp)     ::  kmxm  !<
284       REAL(wp)     ::  kmxp  !<
285       REAL(wp)     ::  kmzm  !<
286       REAL(wp)     ::  kmzp  !<
287
288       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
289       !$acc declare create ( vsus )
290
291!
292!--    First calculate horizontal momentum flux v'u' at vertical walls,
293!--    if neccessary
294       IF ( topography /= 'flat' )  THEN
295          CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp,          &
296                                nzb_v_inner, nzb_v_outer, wall_v )
297       ENDIF
298
299       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
300       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
301       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
302       DO  i = i_left, i_right
303          DO  j = j_south, j_north
304!
305!--          Compute horizontal diffusion
306             DO  k = 1, nzt
307                IF ( k > nzb_v_outer(j,i) )  THEN
308!
309!--                Interpolate eddy diffusivities on staggered gridpoints
310                   kmxp = 0.25_wp *                                            &
311                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
312                   kmxm = 0.25_wp *                                            &
313                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
314
315                   tend(k,j,i) = tend(k,j,i)                                   &
316                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
317                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
318                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
319                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
320                         &   ) * ddx                                           &
321                         & + 2.0_wp * (                                        &
322                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
323                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
324                         &            ) * ddy2
325                ENDIF
326             ENDDO
327
328!
329!--          Wall functions at the left and right walls, respectively
330             DO  k = 1, nzt
331                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
332                    wall_v(j,i) /= 0.0_wp )  THEN
333
334                   kmxp = 0.25_wp *                                            &
335                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
336                   kmxm = 0.25_wp *                                            &
337                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
338                   
339                   tend(k,j,i) = tend(k,j,i)                                   &
340                                 + 2.0_wp * (                                  &
341                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
342                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
343                                            ) * ddy2                           &
344                                 + (   fxp(j,i) * (                            &
345                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
346                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
347                                                  )                            &
348                                     - fxm(j,i) * (                            &
349                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
350                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
351                                                  )                            &
352                                     + wall_v(j,i) * vsus(k,j,i)               &
353                                   ) * ddx
354                ENDIF
355             ENDDO
356
357!
358!--          Compute vertical diffusion. In case of simulating a Prandtl
359!--          layer, index k starts at nzb_v_inner+2.
360             DO  k = 1, nzt_diff
361                IF ( k >= nzb_diff_v(j,i) )  THEN
362!
363!--                Interpolate eddy diffusivities on staggered gridpoints
364                   kmzp = 0.25_wp *                                            &
365                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
366                   kmzm = 0.25_wp *                                            &
367                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
368
369                   tend(k,j,i) = tend(k,j,i)                                   &
370                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
371                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
372                         &            )                                        &
373                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
374                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
375                         &            )                                        &
376                         &   ) * ddzw(k)
377                ENDIF
378             ENDDO
379
380          ENDDO
381       ENDDO
382
383!
384!--    Vertical diffusion at the first grid point above the surface,
385!--    if the momentum flux at the bottom is given by the Prandtl law
386!--    or if it is prescribed by the user.
387!--    Difference quotient of the momentum flux is not formed over
388!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
389!--    comparison with other (LES) models showed that the values of
390!--    the momentum flux becomes too large in this case.
391!--    The term containing w(k-1,..) (see above equation) is removed here
392!--    because the vertical velocity is assumed to be zero at the surface.
393       IF ( use_surface_fluxes )  THEN
394
395          DO  i = i_left, i_right
396             DO  j = j_south, j_north
397         
398                k = nzb_v_inner(j,i)+1
399!
400!--             Interpolate eddy diffusivities on staggered gridpoints
401                kmzp = 0.25_wp *                                               &
402                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
403
404                tend(k,j,i) = tend(k,j,i)                                      &
405                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
406                      &   ) * ddzw(k)                                          &
407                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
408                      &   + vsws(j,i)                                          &
409                      &   ) * ddzw(k)
410             ENDDO
411          ENDDO
412
413       ENDIF
414
415!
416!--    Vertical diffusion at the first gridpoint below the top boundary,
417!--    if the momentum flux at the top is prescribed by the user
418       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
419
420          k = nzt
421
422          DO  i = i_left, i_right
423             DO  j = j_south, j_north
424
425!
426!--             Interpolate eddy diffusivities on staggered gridpoints
427                kmzm = 0.25_wp *                                               &
428                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
429
430                tend(k,j,i) = tend(k,j,i)                                      &
431                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
432                      &   ) * ddzw(k)                                          &
433                      & + ( -vswst(j,i)                                        &
434                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
435                      &   ) * ddzw(k)
436             ENDDO
437          ENDDO
438
439       ENDIF
440       !$acc end kernels
441
442    END SUBROUTINE diffusion_v_acc
443
444
445!------------------------------------------------------------------------------!
446! Description:
447! ------------
448!> Call for grid point i,j
449!------------------------------------------------------------------------------!
450    SUBROUTINE diffusion_v_ij( i, j )
451
452       USE arrays_3d,                                                          &
453           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
454       
455       USE control_parameters,                                                 &
456           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
457       
458       USE grid_variables,                                                     &
459           ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
460       
461       USE indices,                                                            &
462           ONLY:  nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
463       
464       USE kinds
465
466       IMPLICIT NONE
467
468       INTEGER(iwp) ::  i     !<
469       INTEGER(iwp) ::  j     !<
470       INTEGER(iwp) ::  k     !<
471       REAL(wp)     ::  kmxm  !<
472       REAL(wp)     ::  kmxp  !<
473       REAL(wp)     ::  kmzm  !<
474       REAL(wp)     ::  kmzp  !<
475
476       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
477
478!
479!--    Compute horizontal diffusion
480       DO  k = nzb_v_outer(j,i)+1, nzt
481!
482!--       Interpolate eddy diffusivities on staggered gridpoints
483          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
484          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
485
486          tend(k,j,i) = tend(k,j,i)                                            &
487                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
488                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
489                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
490                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
491                      &   ) * ddx                                              &
492                      & + 2.0_wp * (                                           &
493                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
494                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
495                      &            ) * ddy2
496       ENDDO
497
498!
499!--    Wall functions at the left and right walls, respectively
500       IF ( wall_v(j,i) /= 0.0_wp )  THEN
501
502!
503!--       Calculate the horizontal momentum flux v'u'
504          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i),        &
505                            vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
506
507          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
508             kmxp = 0.25_wp *                                                  &
509                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
510             kmxm = 0.25_wp *                                                  &
511                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
512
513             tend(k,j,i) = tend(k,j,i)                                         &
514                                 + 2.0_wp * (                                  &
515                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
516                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
517                                            ) * ddy2                           &
518                                 + (   fxp(j,i) * (                            &
519                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
520                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
521                                                  )                            &
522                                     - fxm(j,i) * (                            &
523                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
524                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
525                                                  )                            &
526                                     + wall_v(j,i) * vsus(k)                   &
527                                   ) * ddx
528          ENDDO
529       ENDIF
530
531!
532!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
533!--    index k starts at nzb_v_inner+2.
534       DO  k = nzb_diff_v(j,i), nzt_diff
535!
536!--       Interpolate eddy diffusivities on staggered gridpoints
537          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
538          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
539
540          tend(k,j,i) = tend(k,j,i)                                            &
541                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
542                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
543                      &            )                                           &
544                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
545                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
546                      &            )                                           &
547                      &   ) * ddzw(k)
548       ENDDO
549
550!
551!--    Vertical diffusion at the first grid point above the surface, if the
552!--    momentum flux at the bottom is given by the Prandtl law or if it is
553!--    prescribed by the user.
554!--    Difference quotient of the momentum flux is not formed over half of
555!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
556!--    other (LES) models showed that the values of the momentum flux becomes
557!--    too large in this case.
558!--    The term containing w(k-1,..) (see above equation) is removed here
559!--    because the vertical velocity is assumed to be zero at the surface.
560       IF ( use_surface_fluxes )  THEN
561          k = nzb_v_inner(j,i)+1
562!
563!--       Interpolate eddy diffusivities on staggered gridpoints
564          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
565
566          tend(k,j,i) = tend(k,j,i)                                            &
567                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
568                      &   ) * ddzw(k)                                          &
569                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
570                      &   + vsws(j,i)                                          &
571                      &   ) * ddzw(k)
572       ENDIF
573
574!
575!--    Vertical diffusion at the first gridpoint below the top boundary,
576!--    if the momentum flux at the top is prescribed by the user
577       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
578          k = nzt
579!
580!--       Interpolate eddy diffusivities on staggered gridpoints
581          kmzm = 0.25_wp * &
582                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
583
584          tend(k,j,i) = tend(k,j,i)                                            &
585                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
586                      &   ) * ddzw(k)                                          &
587                      & + ( -vswst(j,i)                                        &
588                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
589                      &   ) * ddzw(k)
590       ENDIF
591
592    END SUBROUTINE diffusion_v_ij
593
594 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.