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

Last change on this file since 53 was 53, checked in by raasch, 15 years ago

preliminary version, several changes to be explained later

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