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

Last change on this file since 1111 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.2 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 1037 2012-10-22 14:10:22Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 1015 2012-09-27 09:23:24Z raasch
32! accelerator version (*_acc) added
33!
34! 1001 2012-09-13 14:08:46Z raasch
35! arrays comunicated by module instead of parameter list
36!
37! 978 2012-08-09 08:28:32Z fricke
38! outflow damping layer removed
39! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
40!
41! 667 2010-12-23 12:06:00Z suehring/gryschka
42! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
43!
44! 366 2009-08-25 08:06:27Z raasch
45! bc_lr replaced by bc_lr_cyc
46!
47! 106 2007-08-16 14:30:26Z raasch
48! Momentumflux at top (vswst) included as boundary condition,
49! j loop is starting from nysv (needed for non-cyclic boundary conditions)
50!
51! 75 2007-03-22 09:54:05Z raasch
52! Wall functions now include diabatic conditions, call of routine wall_fluxes,
53! z0 removed from argument list, vynp eliminated
54!
55! 20 2007-02-26 00:12:32Z raasch
56! Bugfix: ddzw dimensioned 1:nzt"+1"
57!
58! RCS Log replace by Id keyword, revision history cleaned up
59!
60! Revision 1.15  2006/02/23 10:36:00  raasch
61! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner
62! or nzb_diff_v, respectively, in vertical diffusion,
63! wall functions added for north and south walls, +z0 in argument list,
64! terms containing w(k-1,..) are removed from the Prandtl-layer equation
65! because they cause errors at the edges of topography
66! WARNING: loops containing the MAX function are still not properly vectorized!
67!
68! Revision 1.1  1997/09/12 06:24:01  raasch
69! Initial revision
70!
71!
72! Description:
73! ------------
74! Diffusion term of the v-component
75!------------------------------------------------------------------------------!
76
77    USE wall_fluxes_mod
78
79    PRIVATE
80    PUBLIC diffusion_v, diffusion_v_acc
81
82    INTERFACE diffusion_v
83       MODULE PROCEDURE diffusion_v
84       MODULE PROCEDURE diffusion_v_ij
85    END INTERFACE diffusion_v
86
87    INTERFACE diffusion_v_acc
88       MODULE PROCEDURE diffusion_v_acc
89    END INTERFACE diffusion_v_acc
90
91 CONTAINS
92
93
94!------------------------------------------------------------------------------!
95! Call for all grid points
96!------------------------------------------------------------------------------!
97    SUBROUTINE diffusion_v
98
99       USE arrays_3d
100       USE control_parameters
101       USE grid_variables
102       USE indices
103
104       IMPLICIT NONE
105
106       INTEGER ::  i, j, k
107       REAL    ::  kmxm, kmxp, kmzm, kmzp
108
109       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
110
111!
112!--    First calculate horizontal momentum flux v'u' at vertical walls,
113!--    if neccessary
114       IF ( topography /= 'flat' )  THEN
115          CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
116                            nzb_v_outer, wall_v )
117       ENDIF
118
119       DO  i = nxl, nxr
120          DO  j = nysv, nyn
121!
122!--          Compute horizontal diffusion
123             DO  k = nzb_v_outer(j,i)+1, nzt
124!
125!--             Interpolate eddy diffusivities on staggered gridpoints
126                kmxp = 0.25 * &
127                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
128                kmxm = 0.25 * &
129                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
130
131                tend(k,j,i) = tend(k,j,i)                                    &
132                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
133                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
134                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
135                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
136                      &   ) * ddx                                            &
137                      & + 2.0 * (                                            &
138                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
139                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
140                      &         ) * ddy2
141             ENDDO
142
143!
144!--          Wall functions at the left and right walls, respectively
145             IF ( wall_v(j,i) /= 0.0 )  THEN
146
147                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
148                   kmxp = 0.25 * &
149                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
150                   kmxm = 0.25 * &
151                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
152                   
153                   tend(k,j,i) = tend(k,j,i)                                   &
154                                 + 2.0 * (                                     &
155                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
156                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
157                                         ) * ddy2                              &
158                                 + (   fxp(j,i) * (                            &
159                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
160                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
161                                                  )                            &
162                                     - fxm(j,i) * (                            &
163                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
164                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
165                                                  )                            &
166                                     + wall_v(j,i) * vsus(k,j,i)               &
167                                   ) * ddx
168                ENDDO
169             ENDIF
170
171!
172!--          Compute vertical diffusion. In case of simulating a Prandtl
173!--          layer, index k starts at nzb_v_inner+2.
174             DO  k = nzb_diff_v(j,i), nzt_diff
175!
176!--             Interpolate eddy diffusivities on staggered gridpoints
177                kmzp = 0.25 * &
178                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
179                kmzm = 0.25 * &
180                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
181
182                tend(k,j,i) = tend(k,j,i)                                    &
183                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
184                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
185                      &            )                                         &
186                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
187                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
188                      &            )                                         &
189                      &   ) * ddzw(k)
190             ENDDO
191
192!
193!--          Vertical diffusion at the first grid point above the surface,
194!--          if the momentum flux at the bottom is given by the Prandtl law
195!--          or if it is prescribed by the user.
196!--          Difference quotient of the momentum flux is not formed over
197!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
198!--          comparison with other (LES) modell showed that the values of
199!--          the momentum flux becomes too large in this case.
200!--          The term containing w(k-1,..) (see above equation) is removed here
201!--          because the vertical velocity is assumed to be zero at the surface.
202             IF ( use_surface_fluxes )  THEN
203                k = nzb_v_inner(j,i)+1
204!
205!--             Interpolate eddy diffusivities on staggered gridpoints
206                kmzp = 0.25 * &
207                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
208                kmzm = 0.25 * &
209                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
210
211                tend(k,j,i) = tend(k,j,i)                                    &
212                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
213                      &   ) * ddzw(k)                                        &
214                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
215                      &   + vsws(j,i)                                        &
216                      &   ) * ddzw(k)
217             ENDIF
218
219!
220!--          Vertical diffusion at the first gridpoint below the top boundary,
221!--          if the momentum flux at the top is prescribed by the user
222             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
223                k = nzt
224!
225!--             Interpolate eddy diffusivities on staggered gridpoints
226                kmzp = 0.25 * &
227                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
228                kmzm = 0.25 * &
229                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
230
231                tend(k,j,i) = tend(k,j,i)                                    &
232                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
233                      &   ) * ddzw(k)                                        &
234                      & + ( -vswst(j,i)                                      &
235                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
236                      &   ) * ddzw(k)
237             ENDIF
238
239          ENDDO
240       ENDDO
241
242    END SUBROUTINE diffusion_v
243
244
245!------------------------------------------------------------------------------!
246! Call for all grid points - accelerator version
247!------------------------------------------------------------------------------!
248    SUBROUTINE diffusion_v_acc
249
250       USE arrays_3d
251       USE control_parameters
252       USE grid_variables
253       USE indices
254
255       IMPLICIT NONE
256
257       INTEGER ::  i, j, k
258       REAL    ::  kmxm, kmxp, kmzm, kmzp
259
260       !$acc declare create ( vsus )
261       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
262
263!
264!--    First calculate horizontal momentum flux v'u' at vertical walls,
265!--    if neccessary
266       IF ( topography /= 'flat' )  THEN
267          CALL wall_fluxes_acc( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
268                                nzb_v_outer, wall_v )
269       ENDIF
270
271       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )   &
272       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )           &
273       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
274       !$acc loop
275       DO  i = nxl, nxr
276          DO  j = nysv, nyn
277!
278!--          Compute horizontal diffusion
279             !$acc loop vector(32)
280             DO  k = 1, nzt
281                IF ( k > nzb_v_outer(j,i) )  THEN
282!
283!--                Interpolate eddy diffusivities on staggered gridpoints
284                   kmxp = 0.25 * &
285                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
286                   kmxm = 0.25 * &
287                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
288
289                   tend(k,j,i) = tend(k,j,i)                                   &
290                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
291                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
292                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
293                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
294                         &   ) * ddx                                           &
295                         & + 2.0 * (                                           &
296                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
297                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
298                         &         ) * ddy2
299                ENDIF
300             ENDDO
301
302!
303!--          Wall functions at the left and right walls, respectively
304             !$acc loop vector(32)
305             DO  k = 1, nzt
306                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND. &
307                    wall_v(j,i) /= 0.0 )  THEN
308
309                   kmxp = 0.25 * &
310                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
311                   kmxm = 0.25 * &
312                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
313                   
314                   tend(k,j,i) = tend(k,j,i)                                   &
315                                 + 2.0 * (                                     &
316                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
317                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
318                                         ) * ddy2                              &
319                                 + (   fxp(j,i) * (                            &
320                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
321                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
322                                                  )                            &
323                                     - fxm(j,i) * (                            &
324                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
325                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
326                                                  )                            &
327                                     + wall_v(j,i) * vsus(k,j,i)               &
328                                   ) * ddx
329                ENDIF
330             ENDDO
331
332!
333!--          Compute vertical diffusion. In case of simulating a Prandtl
334!--          layer, index k starts at nzb_v_inner+2.
335             !$acc loop vector(32)
336             DO  k = 1, nzt_diff
337                IF ( k >= nzb_diff_v(j,i) )  THEN
338!
339!--                Interpolate eddy diffusivities on staggered gridpoints
340                   kmzp = 0.25 * &
341                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
342                   kmzm = 0.25 * &
343                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
344
345                   tend(k,j,i) = tend(k,j,i)                                   &
346                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
347                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
348                         &            )                                        &
349                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
350                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
351                         &            )                                        &
352                         &   ) * ddzw(k)
353                ENDIF
354             ENDDO
355
356          ENDDO
357       ENDDO
358
359!
360!--    Vertical diffusion at the first grid point above the surface,
361!--    if the momentum flux at the bottom is given by the Prandtl law
362!--    or if it is prescribed by the user.
363!--    Difference quotient of the momentum flux is not formed over
364!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
365!--    comparison with other (LES) modell showed that the values of
366!--    the momentum flux becomes too large in this case.
367!--    The term containing w(k-1,..) (see above equation) is removed here
368!--    because the vertical velocity is assumed to be zero at the surface.
369       IF ( use_surface_fluxes )  THEN
370
371          !$acc loop
372          DO  i = nxl, nxr
373             !$acc loop vector(32)
374             DO  j = nysv, nyn
375         
376                k = nzb_v_inner(j,i)+1
377!
378!--             Interpolate eddy diffusivities on staggered gridpoints
379                kmzp = 0.25 * &
380                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
381                kmzm = 0.25 * &
382                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
383
384                tend(k,j,i) = tend(k,j,i)                                    &
385                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
386                      &   ) * ddzw(k)                                        &
387                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
388                      &   + vsws(j,i)                                        &
389                      &   ) * ddzw(k)
390             ENDDO
391          ENDDO
392
393       ENDIF
394
395!
396!--    Vertical diffusion at the first gridpoint below the top boundary,
397!--    if the momentum flux at the top is prescribed by the user
398       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
399
400          k = nzt
401
402          !$acc loop
403          DO  i = nxl, nxr
404             !$acc loop vector(32)
405             DO  j = nysv, nyn
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.