source: palm/trunk/SOURCE/diffusion_u.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_u_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_u.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! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
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_ns replaced by bc_ns_cyc
27!
28! 106 2007-08-16 14:30:26Z raasch
29! Momentumflux at top (uswst) included as boundary condition,
30! i loop is starting from nxlu (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, uxrp 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:35:35  raasch
42! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
43! or nzb_diff_u, 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:23:51  raasch
50! Initial revision
51!
52!
53! Description:
54! ------------
55! Diffusion term of the u-component
56! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
57!        and slows down the speed on NEC about 5-10%
58!------------------------------------------------------------------------------!
59
60    USE wall_fluxes_mod
61
62    PRIVATE
63    PUBLIC diffusion_u, diffusion_u_acc
64
65    INTERFACE diffusion_u
66       MODULE PROCEDURE diffusion_u
67       MODULE PROCEDURE diffusion_u_ij
68    END INTERFACE diffusion_u
69
70    INTERFACE diffusion_u_acc
71       MODULE PROCEDURE diffusion_u_acc
72    END INTERFACE diffusion_u_acc
73
74 CONTAINS
75
76
77!------------------------------------------------------------------------------!
78! Call for all grid points
79!------------------------------------------------------------------------------!
80    SUBROUTINE diffusion_u
81
82       USE arrays_3d
83       USE control_parameters
84       USE grid_variables
85       USE indices
86
87       IMPLICIT NONE
88
89       INTEGER ::  i, j, k
90       REAL    ::  kmym, kmyp, kmzm, kmzp
91
92       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
93
94!
95!--    First calculate horizontal momentum flux u'v' at vertical walls,
96!--    if neccessary
97       IF ( topography /= 'flat' )  THEN
98          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
99                            nzb_u_outer, wall_u )
100       ENDIF
101
102       DO  i = nxlu, nxr
103          DO  j = nys, nyn
104!
105!--          Compute horizontal diffusion
106             DO  k = nzb_u_outer(j,i)+1, nzt
107!
108!--             Interpolate eddy diffusivities on staggered gridpoints
109                kmyp = 0.25 * &
110                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
111                kmym = 0.25 * &
112                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
113
114                tend(k,j,i) = tend(k,j,i)                                    &
115                      & + 2.0 * (                                            &
116                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
117                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
118                      &         ) * ddx2                                     &
119                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
120                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
121                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
122                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
123                      &   ) * ddy
124             ENDDO
125
126!
127!--          Wall functions at the north and south walls, respectively
128             IF ( wall_u(j,i) /= 0.0 )  THEN
129
130                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
131                   kmyp = 0.25 * &
132                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
133                   kmym = 0.25 * &
134                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
135
136                   tend(k,j,i) = tend(k,j,i)                                   &
137                                 + 2.0 * (                                     &
138                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
139                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
140                                         ) * ddx2                              &
141                                 + (   fyp(j,i) * (                            &
142                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
143                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
144                                                  )                            &
145                                     - fym(j,i) * (                            &
146                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
147                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
148                                                  )                            &
149                                     + wall_u(j,i) * usvs(k,j,i)               &
150                                   ) * ddy
151                ENDDO
152             ENDIF
153
154!
155!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
156!--          index k starts at nzb_u_inner+2.
157             DO  k = nzb_diff_u(j,i), nzt_diff
158!
159!--             Interpolate eddy diffusivities on staggered gridpoints
160                kmzp = 0.25 * &
161                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
162                kmzm = 0.25 * &
163                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
164
165                tend(k,j,i) = tend(k,j,i)                                    &
166                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
167                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
168                      &            )                                         &
169                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
170                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
171                      &            )                                          &
172                      &   ) * ddzw(k)
173             ENDDO
174
175!
176!--          Vertical diffusion at the first grid point above the surface,
177!--          if the momentum flux at the bottom is given by the Prandtl law or
178!--          if it is prescribed by the user.
179!--          Difference quotient of the momentum flux is not formed over half
180!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
181!--          with other (LES) modell showed that the values of the momentum
182!--          flux becomes too large in this case.
183!--          The term containing w(k-1,..) (see above equation) is removed here
184!--          because the vertical velocity is assumed to be zero at the surface.
185             IF ( use_surface_fluxes )  THEN
186                k = nzb_u_inner(j,i)+1
187!
188!--             Interpolate eddy diffusivities on staggered gridpoints
189                kmzp = 0.25 * &
190                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
191                kmzm = 0.25 * &
192                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
193
194                tend(k,j,i) = tend(k,j,i)                                    &
195                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
196                      &   ) * ddzw(k)                                        &
197                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
198                      &   + usws(j,i)                                        &
199                      &   ) * ddzw(k)
200             ENDIF
201
202!
203!--          Vertical diffusion at the first gridpoint below the top boundary,
204!--          if the momentum flux at the top is prescribed by the user
205             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
206                k = nzt
207!
208!--             Interpolate eddy diffusivities on staggered gridpoints
209                kmzp = 0.25 * &
210                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
211                kmzm = 0.25 * &
212                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
213
214                tend(k,j,i) = tend(k,j,i)                                    &
215                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
216                      &   ) * ddzw(k)                                        &
217                      & + ( -uswst(j,i)                                      &
218                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
219                      &   ) * ddzw(k)
220             ENDIF
221
222          ENDDO
223       ENDDO
224
225    END SUBROUTINE diffusion_u
226
227
228!------------------------------------------------------------------------------!
229! Call for all grid points - accelerator version
230!------------------------------------------------------------------------------!
231    SUBROUTINE diffusion_u_acc
232
233       USE arrays_3d
234       USE control_parameters
235       USE grid_variables
236       USE indices
237
238       IMPLICIT NONE
239
240       INTEGER ::  i, j, k
241       REAL    ::  kmym, kmyp, kmzm, kmzp
242
243       !$acc declare create ( usvs )
244       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
245
246!
247!--    First calculate horizontal momentum flux u'v' at vertical walls,
248!--    if neccessary
249       IF ( topography /= 'flat' )  THEN
250          CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
251                                nzb_u_outer, wall_u )
252       ENDIF
253
254       !$acc kernels present ( u, v, w, km, tend, usws, uswst )   &
255       !$acc         present ( ddzu, ddzw, fym, fyp, wall_u )           &
256       !$acc         present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
257       !$acc loop
258       DO  i = nxlu, nxr
259          DO  j = nys, nyn
260!
261!--          Compute horizontal diffusion
262             !$acc loop vector(32)
263             DO  k = 1, nzt
264                IF ( k > nzb_u_outer(j,i) )  THEN
265!
266!--                Interpolate eddy diffusivities on staggered gridpoints
267                   kmyp = 0.25 * &
268                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
269                   kmym = 0.25 * &
270                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
271
272                   tend(k,j,i) = tend(k,j,i)                                   &
273                         & + 2.0 * (                                           &
274                         &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
275                         &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
276                         &         ) * ddx2                                    &
277                         & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
278                         &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
279                         &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
280                         &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
281                         &   ) * ddy
282                ENDIF
283             ENDDO
284
285!
286!--          Wall functions at the north and south walls, respectively
287             !$acc loop vector(32)
288             DO  k = 1, nzt
289                IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND. &
290                    wall_u(j,i) /= 0.0 )  THEN
291
292                   kmyp = 0.25 * &
293                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
294                   kmym = 0.25 * &
295                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
296
297                   tend(k,j,i) = tend(k,j,i)                                   &
298                                 + 2.0 * (                                     &
299                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
300                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
301                                         ) * ddx2                              &
302                                 + (   fyp(j,i) * (                            &
303                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
304                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
305                                                  )                            &
306                                     - fym(j,i) * (                            &
307                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
308                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
309                                                  )                            &
310                                     + wall_u(j,i) * usvs(k,j,i)               &
311                                   ) * ddy
312                ENDIF
313             ENDDO
314
315!
316!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
317!--          index k starts at nzb_u_inner+2.
318             !$acc loop vector(32)
319             DO  k = 1, nzt_diff
320                IF ( k >= nzb_diff_u(j,i) )  THEN
321!
322!--                Interpolate eddy diffusivities on staggered gridpoints
323                   kmzp = 0.25 * &
324                          ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
325                   kmzm = 0.25 * &
326                          ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
327
328                   tend(k,j,i) = tend(k,j,i)                                   &
329                         & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
330                         &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
331                         &            )                                        &
332                         &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
333                         &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
334                         &            )                                        &
335                         &   ) * ddzw(k)
336                ENDIF
337             ENDDO
338
339          ENDDO
340       ENDDO
341
342!
343!--    Vertical diffusion at the first grid point above the surface,
344!--    if the momentum flux at the bottom is given by the Prandtl law or
345!--    if it is prescribed by the user.
346!--    Difference quotient of the momentum flux is not formed over half
347!--    of the grid spacing (2.0*ddzw(k)) any more, since the comparison
348!--    with other (LES) modell showed that the values of the momentum
349!--    flux becomes too large in this case.
350!--    The term containing w(k-1,..) (see above equation) is removed here
351!--    because the vertical velocity is assumed to be zero at the surface.
352       IF ( use_surface_fluxes )  THEN
353
354          !$acc loop
355          DO  i = nxlu, nxr
356             !$acc loop vector(32)
357             DO  j = nys, nyn
358         
359                k = nzb_u_inner(j,i)+1
360!
361!--             Interpolate eddy diffusivities on staggered gridpoints
362                kmzp = 0.25 * &
363                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
364                kmzm = 0.25 * &
365                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
366
367                tend(k,j,i) = tend(k,j,i)                                    &
368                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
369                      &   ) * ddzw(k)                                        &
370                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
371                      &   + usws(j,i)                                        &
372                      &   ) * ddzw(k)
373             ENDDO
374          ENDDO
375
376       ENDIF
377
378!
379!--    Vertical diffusion at the first gridpoint below the top boundary,
380!--    if the momentum flux at the top is prescribed by the user
381       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
382
383          k = nzt
384
385          !$acc loop
386          DO  i = nxlu, nxr
387             !$acc loop vector(32)
388             DO  j = nys, nyn
389
390!
391!--             Interpolate eddy diffusivities on staggered gridpoints
392                kmzp = 0.25 * &
393                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
394                kmzm = 0.25 * &
395                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
396
397                tend(k,j,i) = tend(k,j,i)                                    &
398                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
399                      &   ) * ddzw(k)                                        &
400                      & + ( -uswst(j,i)                                      &
401                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
402                      &   ) * ddzw(k)
403             ENDDO
404          ENDDO
405
406       ENDIF
407       !$acc end kernels
408
409    END SUBROUTINE diffusion_u_acc
410
411
412!------------------------------------------------------------------------------!
413! Call for grid point i,j
414!------------------------------------------------------------------------------!
415    SUBROUTINE diffusion_u_ij( i, j )
416
417       USE arrays_3d
418       USE control_parameters
419       USE grid_variables
420       USE indices
421
422       IMPLICIT NONE
423
424       INTEGER ::  i, j, k
425       REAL    ::  kmym, kmyp, kmzm, kmzp
426
427       REAL, DIMENSION(nzb:nzt+1) ::  usvs
428
429!
430!--    Compute horizontal diffusion
431       DO  k = nzb_u_outer(j,i)+1, nzt
432!
433!--       Interpolate eddy diffusivities on staggered gridpoints
434          kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
435          kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
436
437          tend(k,j,i) = tend(k,j,i)                                          &
438                      & + 2.0 * (                                            &
439                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
440                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
441                      &         ) * ddx2                                     &
442                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
443                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
444                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
445                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
446                      &   ) * ddy
447       ENDDO
448
449!
450!--    Wall functions at the north and south walls, respectively
451       IF ( wall_u(j,i) .NE. 0.0 )  THEN
452
453!
454!--       Calculate the horizontal momentum flux u'v'
455          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
456                            usvs, 1.0, 0.0, 0.0, 0.0 )
457
458          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
459             kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
460             kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
461
462             tend(k,j,i) = tend(k,j,i)                                         &
463                                 + 2.0 * (                                     &
464                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
465                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
466                                         ) * ddx2                              &
467                                 + (   fyp(j,i) * (                            &
468                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
469                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
470                                                  )                            &
471                                     - fym(j,i) * (                            &
472                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
473                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
474                                                  )                            &
475                                     + wall_u(j,i) * usvs(k)                   &
476                                   ) * ddy
477          ENDDO
478       ENDIF
479
480!
481!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
482!--    index k starts at nzb_u_inner+2.
483       DO  k = nzb_diff_u(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,i-1)+km(k+1,j,i-1) )
487          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
488
489          tend(k,j,i) = tend(k,j,i)                                          &
490                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
491                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
492                      &            )                                         &
493                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
494                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
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_u_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,i-1)+km(k+1,j,i-1) )
514          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
515
516          tend(k,j,i) = tend(k,j,i)                                          &
517                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
518                      &   ) * ddzw(k)                                        &
519                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
520                      &   + usws(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,i-1)+km(k+1,j,i-1) )
533          kmzm = 0.25 * &
534                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
535
536          tend(k,j,i) = tend(k,j,i)                                          &
537                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
538                      &   ) * ddzw(k)                                        &
539                      & + ( -uswst(j,i)                                      &
540                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
541                      &   ) * ddzw(k)
542       ENDIF
543
544    END SUBROUTINE diffusion_u_ij
545
546 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.