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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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