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

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

r1028 documented

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