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

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