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

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

preliminary update for changes concerning non-cyclic boundary conditions

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