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

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

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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