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

Last change on this file since 106 was 106, checked in by raasch, 14 years ago

preliminary update of bugfixes and extensions for non-cyclic BCs

  • Property svn:keywords set to Id
File size: 17.6 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Momentumflux at top (uswst) included as boundary condition,
7! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
8!
9! Former revisions:
10! -----------------
11! $Id: diffusion_u.f90 106 2007-08-16 14:30:26Z raasch $
12!
13! 75 2007-03-22 09:54:05Z raasch
14! Wall functions now include diabatic conditions, call of routine wall_fluxes,
15! z0 removed from argument list, uxrp eliminated
16!
17! 20 2007-02-26 00:12:32Z raasch
18! Bugfix: ddzw dimensioned 1:nzt"+1"
19!
20! RCS Log replace by Id keyword, revision history cleaned up
21!
22! Revision 1.15  2006/02/23 10:35:35  raasch
23! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
24! or nzb_diff_u, respectively, in vertical diffusion,
25! wall functions added for north and south walls, +z0 in argument list,
26! terms containing w(k-1,..) are removed from the Prandtl-layer equation
27! because they cause errors at the edges of topography
28! WARNING: loops containing the MAX function are still not properly vectorized!
29!
30! Revision 1.1  1997/09/12 06:23:51  raasch
31! Initial revision
32!
33!
34! Description:
35! ------------
36! Diffusion term of the u-component
37! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
38!        and slows down the speed on NEC about 5-10%
39!------------------------------------------------------------------------------!
40
41    USE wall_fluxes_mod
42
43    PRIVATE
44    PUBLIC diffusion_u
45
46    INTERFACE diffusion_u
47       MODULE PROCEDURE diffusion_u
48       MODULE PROCEDURE diffusion_u_ij
49    END INTERFACE diffusion_u
50
51 CONTAINS
52
53
54!------------------------------------------------------------------------------!
55! Call for all grid points
56!------------------------------------------------------------------------------!
57    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
58                            v, w )
59
60       USE control_parameters
61       USE grid_variables
62       USE indices
63
64       IMPLICIT NONE
65
66       INTEGER ::  i, j, k
67       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
68       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
69       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
70       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
71       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
72       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
73
74!
75!--    First calculate horizontal momentum flux u'v' at vertical walls,
76!--    if neccessary
77       IF ( topography /= 'flat' )  THEN
78          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
79                            nzb_u_outer, wall_u )
80       ENDIF
81
82       DO  i = nxlu, nxr
83          DO  j = nys,nyn
84!
85!--          Compute horizontal diffusion
86             DO  k = nzb_u_outer(j,i)+1, nzt
87!
88!--             Interpolate eddy diffusivities on staggered gridpoints
89                kmyp_x = 0.25 * &
90                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
91                kmym_x = 0.25 * &
92                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
93                kmyp_y = kmyp_x
94                kmym_y = kmym_x
95!
96!--             Increase diffusion at the outflow boundary in case of
97!--             non-cyclic lateral boundaries. Damping is only needed for
98!--             velocity components parallel to the outflow boundary in
99!--             the direction normal to the outflow boundary.
100                IF ( bc_ns /= 'cyclic' )  THEN
101                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
102                   kmym_y = MAX( kmym_y, km_damp_y(j) )
103                ENDIF
104
105                tend(k,j,i) = tend(k,j,i)                                    &
106                      & + 2.0 * (                                            &
107                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
108                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
109                      &         ) * ddx2                                     &
110                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
111                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
112                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
113                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
114                      &   ) * ddy
115             ENDDO
116
117!
118!--          Wall functions at the north and south walls, respectively
119             IF ( wall_u(j,i) /= 0.0 )  THEN
120
121                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
122                   kmyp_x = 0.25 * &
123                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
124                   kmym_x = 0.25 * &
125                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
126                   kmyp_y = kmyp_x
127                   kmym_y = kmym_x
128!
129!--                Increase diffusion at the outflow boundary in case of
130!--                non-cyclic lateral boundaries. Damping is only needed for
131!--                velocity components parallel to the outflow boundary in
132!--                the direction normal to the outflow boundary.
133                   IF ( bc_ns /= 'cyclic' )  THEN
134                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
135                      kmym_y = MAX( kmym_y, km_damp_y(j) )
136                   ENDIF
137
138                   tend(k,j,i) = tend(k,j,i)                                   &
139                                 + 2.0 * (                                     &
140                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
141                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
142                                         ) * ddx2                              &
143                                 + (   fyp(j,i) * (                            &
144                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
145                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
146                                                  )                            &
147                                     - fym(j,i) * (                            &
148                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
149                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
150                                                  )                            &
151                                     + wall_u(j,i) * usvs(k,j,i)               &
152                                   ) * ddy
153                ENDDO
154             ENDIF
155
156!
157!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
158!--          index k starts at nzb_u_inner+2.
159             DO  k = nzb_diff_u(j,i), nzt_diff
160!
161!--             Interpolate eddy diffusivities on staggered gridpoints
162                kmzp = 0.25 * &
163                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
164                kmzm = 0.25 * &
165                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
166
167                tend(k,j,i) = tend(k,j,i)                                    &
168                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
169                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
170                      &            )                                         &
171                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
172                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
173                      &            )                                         &
174                      &   ) * ddzw(k)
175             ENDDO
176
177!
178!--          Vertical diffusion at the first grid point above the surface,
179!--          if the momentum flux at the bottom is given by the Prandtl law or
180!--          if it is prescribed by the user.
181!--          Difference quotient of the momentum flux is not formed over half
182!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
183!--          with other (LES) modell showed that the values of the momentum
184!--          flux becomes too large in this case.
185!--          The term containing w(k-1,..) (see above equation) is removed here
186!--          because the vertical velocity is assumed to be zero at the surface.
187             IF ( use_surface_fluxes )  THEN
188                k = nzb_u_inner(j,i)+1
189!
190!--             Interpolate eddy diffusivities on staggered gridpoints
191                kmzp = 0.25 * &
192                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
193                kmzm = 0.25 * &
194                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
195
196                tend(k,j,i) = tend(k,j,i)                                    &
197                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
198                      &   ) * ddzw(k)                                        &
199                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
200                      &   + usws(j,i)                                        &
201                      &   ) * ddzw(k)
202             ENDIF
203
204!
205!--          Vertical diffusion at the first gridpoint below the top boundary,
206!--          if the momentum flux at the top is prescribed by the user
207             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
208                k = nzt
209!
210!--             Interpolate eddy diffusivities on staggered gridpoints
211                kmzp = 0.25 * &
212                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
213                kmzm = 0.25 * &
214                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
215
216                tend(k,j,i) = tend(k,j,i)                                    &
217                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
218                      &   ) * ddzw(k)                                        &
219                      & + ( -uswst(j,i)                                      &
220                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
221                      &   ) * ddzw(k)
222             ENDIF
223
224          ENDDO
225       ENDDO
226
227    END SUBROUTINE diffusion_u
228
229
230!------------------------------------------------------------------------------!
231! Call for grid point i,j
232!------------------------------------------------------------------------------!
233    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
234                               uswst, v, w )
235
236       USE control_parameters
237       USE grid_variables
238       USE indices
239
240       IMPLICIT NONE
241
242       INTEGER ::  i, j, k
243       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
244       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
245       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
246       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
247       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
248       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
249
250!
251!--    Compute horizontal diffusion
252       DO  k = nzb_u_outer(j,i)+1, nzt
253!
254!--       Interpolate eddy diffusivities on staggered gridpoints
255          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
256          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
257          kmyp_y = kmyp_x
258          kmym_y = kmym_x
259
260!
261!--       Increase diffusion at the outflow boundary in case of non-cyclic
262!--       lateral boundaries. Damping is only needed for velocity components
263!--       parallel to the outflow boundary in the direction normal to the
264!--       outflow boundary.
265          IF ( bc_ns /= 'cyclic' )  THEN
266             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
267             kmym_y = MAX( kmym_y, km_damp_y(j) )
268          ENDIF
269
270          tend(k,j,i) = tend(k,j,i)                                          &
271                      & + 2.0 * (                                            &
272                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
273                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
274                      &         ) * ddx2                                     &
275                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
276                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
277                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
278                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
279                      &   ) * ddy
280       ENDDO
281
282!
283!--    Wall functions at the north and south walls, respectively
284       IF ( wall_u(j,i) .NE. 0.0 )  THEN
285
286!
287!--       Calculate the horizontal momentum flux u'v'
288          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
289                            usvs, 1.0, 0.0, 0.0, 0.0 )
290
291          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
292             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
293             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
294             kmyp_y = kmyp_x
295             kmym_y = kmym_x
296!
297!--          Increase diffusion at the outflow boundary in case of
298!--          non-cyclic lateral boundaries. Damping is only needed for
299!--          velocity components parallel to the outflow boundary in
300!--          the direction normal to the outflow boundary.
301             IF ( bc_ns /= 'cyclic' )  THEN
302                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
303                kmym_y = MAX( kmym_y, km_damp_y(j) )
304             ENDIF
305
306             tend(k,j,i) = tend(k,j,i)                                         &
307                                 + 2.0 * (                                     &
308                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
309                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
310                                         ) * ddx2                              &
311                                 + (   fyp(j,i) * (                            &
312                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
313                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
314                                                  )                            &
315                                     - fym(j,i) * (                            &
316                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
317                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
318                                                  )                            &
319                                     + wall_u(j,i) * usvs(k)                   &
320                                   ) * ddy
321          ENDDO
322       ENDIF
323
324!
325!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
326!--    index k starts at nzb_u_inner+2.
327       DO  k = nzb_diff_u(j,i), nzt_diff
328!
329!--       Interpolate eddy diffusivities on staggered gridpoints
330          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
331          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
332
333          tend(k,j,i) = tend(k,j,i)                                          &
334                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
335                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
336                      &            )                                         &
337                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
338                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
339                      &            )                                         &
340                      &   ) * ddzw(k)
341       ENDDO
342
343!
344!--    Vertical diffusion at the first grid point above the surface, if the
345!--    momentum flux at the bottom is given by the Prandtl law or if it is
346!--    prescribed by the user.
347!--    Difference quotient of the momentum flux is not formed over half of
348!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
349!--    other (LES) modell showed that the values of the momentum flux becomes
350!--    too large in this case.
351!--    The term containing w(k-1,..) (see above equation) is removed here
352!--    because the vertical velocity is assumed to be zero at the surface.
353       IF ( use_surface_fluxes )  THEN
354          k = nzb_u_inner(j,i)+1
355!
356!--       Interpolate eddy diffusivities on staggered gridpoints
357          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
358          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
359
360          tend(k,j,i) = tend(k,j,i)                                          &
361                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
362                      &   ) * ddzw(k)                                        &
363                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
364                      &   + usws(j,i)                                        &
365                      &   ) * ddzw(k)
366       ENDIF
367
368!
369!--    Vertical diffusion at the first gridpoint below the top boundary,
370!--    if the momentum flux at the top is prescribed by the user
371       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
372          k = nzt
373!
374!--       Interpolate eddy diffusivities on staggered gridpoints
375          kmzp = 0.25 * &
376                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
377          kmzm = 0.25 * &
378                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
379
380          tend(k,j,i) = tend(k,j,i)                                          &
381                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
382                      &   ) * ddzw(k)                                        &
383                      & + ( -uswst(j,i)                                      &
384                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
385                      &   ) * ddzw(k)
386       ENDIF
387
388    END SUBROUTINE diffusion_u_ij
389
390 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.