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

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

further preliminary changes concerning coupling

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