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

Last change on this file since 1128 was 1128, checked in by raasch, 11 years ago

asynchronous transfer of ghost point data for acc-optimized version

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