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

Last change on this file since 1258 was 1258, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.3 KB
Line 
1 MODULE diffusion_v_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_v.f90 1258 2013-11-08 16:09:09Z raasch $
27!
28! 1257 2013-11-08 15:18:40Z raasch
29! openacc loop and loop vector clauses removed, declare create moved after
30! the FORTRAN declaration statement
31!
32! 1128 2013-04-12 06:19:32Z raasch
33! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
34! j_north
35!
36! 1036 2012-10-22 13:43:42Z raasch
37! code put under GPL (PALM 3.9)
38!
39! 1015 2012-09-27 09:23:24Z raasch
40! accelerator version (*_acc) added
41!
42! 1001 2012-09-13 14:08:46Z raasch
43! arrays comunicated by module instead of parameter list
44!
45! 978 2012-08-09 08:28:32Z fricke
46! outflow damping layer removed
47! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
48!
49! 667 2010-12-23 12:06:00Z suehring/gryschka
50! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
51!
52! 366 2009-08-25 08:06:27Z raasch
53! bc_lr replaced by bc_lr_cyc
54!
55! 106 2007-08-16 14:30:26Z raasch
56! Momentumflux at top (vswst) included as boundary condition,
57! j loop is starting from nysv (needed for non-cyclic boundary conditions)
58!
59! 75 2007-03-22 09:54:05Z raasch
60! Wall functions now include diabatic conditions, call of routine wall_fluxes,
61! z0 removed from argument list, vynp eliminated
62!
63! 20 2007-02-26 00:12:32Z raasch
64! Bugfix: ddzw dimensioned 1:nzt"+1"
65!
66! RCS Log replace by Id keyword, revision history cleaned up
67!
68! Revision 1.15  2006/02/23 10:36:00  raasch
69! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner
70! or nzb_diff_v, respectively, in vertical diffusion,
71! wall functions added for north and south walls, +z0 in argument list,
72! terms containing w(k-1,..) are removed from the Prandtl-layer equation
73! because they cause errors at the edges of topography
74! WARNING: loops containing the MAX function are still not properly vectorized!
75!
76! Revision 1.1  1997/09/12 06:24:01  raasch
77! Initial revision
78!
79!
80! Description:
81! ------------
82! Diffusion term of the v-component
83!------------------------------------------------------------------------------!
84
85    USE wall_fluxes_mod
86
87    PRIVATE
88    PUBLIC diffusion_v, diffusion_v_acc
89
90    INTERFACE diffusion_v
91       MODULE PROCEDURE diffusion_v
92       MODULE PROCEDURE diffusion_v_ij
93    END INTERFACE diffusion_v
94
95    INTERFACE diffusion_v_acc
96       MODULE PROCEDURE diffusion_v_acc
97    END INTERFACE diffusion_v_acc
98
99 CONTAINS
100
101
102!------------------------------------------------------------------------------!
103! Call for all grid points
104!------------------------------------------------------------------------------!
105    SUBROUTINE diffusion_v
106
107       USE arrays_3d
108       USE control_parameters
109       USE grid_variables
110       USE indices
111
112       IMPLICIT NONE
113
114       INTEGER ::  i, j, k
115       REAL    ::  kmxm, kmxp, kmzm, kmzp
116
117       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
118
119!
120!--    First calculate horizontal momentum flux v'u' at vertical walls,
121!--    if neccessary
122       IF ( topography /= 'flat' )  THEN
123          CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
124                            nzb_v_outer, wall_v )
125       ENDIF
126
127       DO  i = nxl, nxr
128          DO  j = nysv, nyn
129!
130!--          Compute horizontal diffusion
131             DO  k = nzb_v_outer(j,i)+1, nzt
132!
133!--             Interpolate eddy diffusivities on staggered gridpoints
134                kmxp = 0.25 * &
135                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
136                kmxm = 0.25 * &
137                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
138
139                tend(k,j,i) = tend(k,j,i)                                    &
140                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
141                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
142                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
143                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
144                      &   ) * ddx                                            &
145                      & + 2.0 * (                                            &
146                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
147                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
148                      &         ) * ddy2
149             ENDDO
150
151!
152!--          Wall functions at the left and right walls, respectively
153             IF ( wall_v(j,i) /= 0.0 )  THEN
154
155                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
156                   kmxp = 0.25 * &
157                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
158                   kmxm = 0.25 * &
159                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
160                   
161                   tend(k,j,i) = tend(k,j,i)                                   &
162                                 + 2.0 * (                                     &
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                                 + (   fxp(j,i) * (                            &
167                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
168                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
169                                                  )                            &
170                                     - fxm(j,i) * (                            &
171                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
172                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
173                                                  )                            &
174                                     + wall_v(j,i) * vsus(k,j,i)               &
175                                   ) * ddx
176                ENDDO
177             ENDIF
178
179!
180!--          Compute vertical diffusion. In case of simulating a Prandtl
181!--          layer, index k starts at nzb_v_inner+2.
182             DO  k = nzb_diff_v(j,i), nzt_diff
183!
184!--             Interpolate eddy diffusivities on staggered gridpoints
185                kmzp = 0.25 * &
186                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
187                kmzm = 0.25 * &
188                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
189
190                tend(k,j,i) = tend(k,j,i)                                    &
191                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
192                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
193                      &            )                                         &
194                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
195                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
196                      &            )                                         &
197                      &   ) * ddzw(k)
198             ENDDO
199
200!
201!--          Vertical diffusion at the first grid point above the surface,
202!--          if the momentum flux at the bottom is given by the Prandtl law
203!--          or if it is prescribed by the user.
204!--          Difference quotient of the momentum flux is not formed over
205!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
206!--          comparison with other (LES) modell showed that the values of
207!--          the momentum flux becomes too large in this case.
208!--          The term containing w(k-1,..) (see above equation) is removed here
209!--          because the vertical velocity is assumed to be zero at the surface.
210             IF ( use_surface_fluxes )  THEN
211                k = nzb_v_inner(j,i)+1
212!
213!--             Interpolate eddy diffusivities on staggered gridpoints
214                kmzp = 0.25 * &
215                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
216                kmzm = 0.25 * &
217                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
218
219                tend(k,j,i) = tend(k,j,i)                                    &
220                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
221                      &   ) * ddzw(k)                                        &
222                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
223                      &   + vsws(j,i)                                        &
224                      &   ) * ddzw(k)
225             ENDIF
226
227!
228!--          Vertical diffusion at the first gridpoint below the top boundary,
229!--          if the momentum flux at the top is prescribed by the user
230             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
231                k = nzt
232!
233!--             Interpolate eddy diffusivities on staggered gridpoints
234                kmzp = 0.25 * &
235                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
236                kmzm = 0.25 * &
237                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
238
239                tend(k,j,i) = tend(k,j,i)                                    &
240                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
241                      &   ) * ddzw(k)                                        &
242                      & + ( -vswst(j,i)                                      &
243                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
244                      &   ) * ddzw(k)
245             ENDIF
246
247          ENDDO
248       ENDDO
249
250    END SUBROUTINE diffusion_v
251
252
253!------------------------------------------------------------------------------!
254! Call for all grid points - accelerator version
255!------------------------------------------------------------------------------!
256    SUBROUTINE diffusion_v_acc
257
258       USE arrays_3d
259       USE control_parameters
260       USE grid_variables
261       USE indices
262
263       IMPLICIT NONE
264
265       INTEGER ::  i, j, k
266       REAL    ::  kmxm, kmxp, kmzm, kmzp
267
268       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
269       !$acc declare create ( vsus )
270
271!
272!--    First calculate horizontal momentum flux v'u' at vertical walls,
273!--    if neccessary
274       IF ( topography /= 'flat' )  THEN
275          CALL wall_fluxes_acc( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
276                                nzb_v_outer, wall_v )
277       ENDIF
278
279       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )   &
280       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )           &
281       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
282       DO  i = i_left, i_right
283          DO  j = j_south, j_north
284!
285!--          Compute horizontal diffusion
286             DO  k = 1, nzt
287                IF ( k > nzb_v_outer(j,i) )  THEN
288!
289!--                Interpolate eddy diffusivities on staggered gridpoints
290                   kmxp = 0.25 * &
291                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
292                   kmxm = 0.25 * &
293                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
294
295                   tend(k,j,i) = tend(k,j,i)                                   &
296                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
297                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
298                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
299                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
300                         &   ) * ddx                                           &
301                         & + 2.0 * (                                           &
302                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
303                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
304                         &         ) * ddy2
305                ENDIF
306             ENDDO
307
308!
309!--          Wall functions at the left and right walls, respectively
310             DO  k = 1, nzt
311                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND. &
312                    wall_v(j,i) /= 0.0 )  THEN
313
314                   kmxp = 0.25 * &
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 * &
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                                 + 2.0 * (                                     &
321                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
322                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
323                                         ) * ddy2                              &
324                                 + (   fxp(j,i) * (                            &
325                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
326                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
327                                                  )                            &
328                                     - fxm(j,i) * (                            &
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                                                  )                            &
332                                     + wall_v(j,i) * vsus(k,j,i)               &
333                                   ) * ddx
334                ENDIF
335             ENDDO
336
337!
338!--          Compute vertical diffusion. In case of simulating a Prandtl
339!--          layer, index k starts at nzb_v_inner+2.
340             DO  k = 1, nzt_diff
341                IF ( k >= nzb_diff_v(j,i) )  THEN
342!
343!--                Interpolate eddy diffusivities on staggered gridpoints
344                   kmzp = 0.25 * &
345                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
346                   kmzm = 0.25 * &
347                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
348
349                   tend(k,j,i) = tend(k,j,i)                                   &
350                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
351                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
352                         &            )                                        &
353                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
354                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
355                         &            )                                        &
356                         &   ) * ddzw(k)
357                ENDIF
358             ENDDO
359
360          ENDDO
361       ENDDO
362
363!
364!--    Vertical diffusion at the first grid point above the surface,
365!--    if the momentum flux at the bottom is given by the Prandtl law
366!--    or if it is prescribed by the user.
367!--    Difference quotient of the momentum flux is not formed over
368!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
369!--    comparison with other (LES) modell showed that the values of
370!--    the momentum flux becomes too large in this case.
371!--    The term containing w(k-1,..) (see above equation) is removed here
372!--    because the vertical velocity is assumed to be zero at the surface.
373       IF ( use_surface_fluxes )  THEN
374
375          DO  i = i_left, i_right
376             DO  j = j_south, j_north
377         
378                k = nzb_v_inner(j,i)+1
379!
380!--             Interpolate eddy diffusivities on staggered gridpoints
381                kmzp = 0.25 * &
382                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
383                kmzm = 0.25 * &
384                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
385
386                tend(k,j,i) = tend(k,j,i)                                    &
387                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
388                      &   ) * ddzw(k)                                        &
389                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
390                      &   + vsws(j,i)                                        &
391                      &   ) * ddzw(k)
392             ENDDO
393          ENDDO
394
395       ENDIF
396
397!
398!--    Vertical diffusion at the first gridpoint below the top boundary,
399!--    if the momentum flux at the top is prescribed by the user
400       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
401
402          k = nzt
403
404          DO  i = i_left, i_right
405             DO  j = j_south, j_north
406
407!
408!--             Interpolate eddy diffusivities on staggered gridpoints
409                kmzp = 0.25 * &
410                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
411                kmzm = 0.25 * &
412                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
413
414                tend(k,j,i) = tend(k,j,i)                                    &
415                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
416                      &   ) * ddzw(k)                                        &
417                      & + ( -vswst(j,i)                                      &
418                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
419                      &   ) * ddzw(k)
420             ENDDO
421          ENDDO
422
423       ENDIF
424       !$acc end kernels
425
426    END SUBROUTINE diffusion_v_acc
427
428
429!------------------------------------------------------------------------------!
430! Call for grid point i,j
431!------------------------------------------------------------------------------!
432    SUBROUTINE diffusion_v_ij( i, j )
433
434       USE arrays_3d
435       USE control_parameters
436       USE grid_variables
437       USE indices
438
439       IMPLICIT NONE
440
441       INTEGER ::  i, j, k
442       REAL    ::  kmxm, kmxp, kmzm, kmzp
443
444       REAL, DIMENSION(nzb:nzt+1) ::  vsus
445
446!
447!--    Compute horizontal diffusion
448       DO  k = nzb_v_outer(j,i)+1, nzt
449!
450!--       Interpolate eddy diffusivities on staggered gridpoints
451          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
452          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
453
454          tend(k,j,i) = tend(k,j,i)                                          &
455                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
456                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
457                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
458                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
459                      &   ) * ddx                                            &
460                      & + 2.0 * (                                            &
461                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
462                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
463                      &         ) * ddy2
464       ENDDO
465
466!
467!--    Wall functions at the left and right walls, respectively
468       IF ( wall_v(j,i) /= 0.0 )  THEN
469
470!
471!--       Calculate the horizontal momentum flux v'u'
472          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
473                            vsus, 0.0, 1.0, 0.0, 0.0 )
474
475          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
476             kmxp = 0.25 * &
477                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
478             kmxm = 0.25 * &
479                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
480
481             tend(k,j,i) = tend(k,j,i)                                         &
482                                 + 2.0 * (                                     &
483                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
484                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
485                                         ) * ddy2                              &
486                                 + (   fxp(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                                                  )                            &
490                                     - fxm(j,i) * (                            &
491                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
492                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
493                                                  )                            &
494                                     + wall_v(j,i) * vsus(k)                   &
495                                   ) * ddx
496          ENDDO
497       ENDIF
498
499!
500!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
501!--    index k starts at nzb_v_inner+2.
502       DO  k = nzb_diff_v(j,i), nzt_diff
503!
504!--       Interpolate eddy diffusivities on staggered gridpoints
505          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
506          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
507
508          tend(k,j,i) = tend(k,j,i)                                          &
509                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
510                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
511                      &            )                                         &
512                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
513                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
514                      &            )                                         &
515                      &   ) * ddzw(k)
516       ENDDO
517
518!
519!--    Vertical diffusion at the first grid point above the surface, if the
520!--    momentum flux at the bottom is given by the Prandtl law or if it is
521!--    prescribed by the user.
522!--    Difference quotient of the momentum flux is not formed over half of
523!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
524!--    other (LES) modell showed that the values of the momentum flux becomes
525!--    too large in this case.
526!--    The term containing w(k-1,..) (see above equation) is removed here
527!--    because the vertical velocity is assumed to be zero at the surface.
528       IF ( use_surface_fluxes )  THEN
529          k = nzb_v_inner(j,i)+1
530!
531!--       Interpolate eddy diffusivities on staggered gridpoints
532          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
533          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
534
535          tend(k,j,i) = tend(k,j,i)                                          &
536                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
537                      &   ) * ddzw(k)                                        &
538                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
539                      &   + vsws(j,i)                                        &
540                      &   ) * ddzw(k)
541       ENDIF
542
543!
544!--    Vertical diffusion at the first gridpoint below the top boundary,
545!--    if the momentum flux at the top is prescribed by the user
546       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
547          k = nzt
548!
549!--       Interpolate eddy diffusivities on staggered gridpoints
550          kmzp = 0.25 * &
551                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
552          kmzm = 0.25 * &
553                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
554
555          tend(k,j,i) = tend(k,j,i)                                          &
556                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
557                      &   ) * ddzw(k)                                        &
558                      & + ( -vswst(j,i)                                      &
559                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
560                      &   ) * ddzw(k)
561       ENDIF
562
563    END SUBROUTINE diffusion_v_ij
564
565 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.