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

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

last commit documented

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