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

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

RCS Log replace by Id keyword, revision history cleaned up

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