source: palm/trunk/SOURCE/diffusion_w.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.1 KB
Line 
1 MODULE diffusion_w_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.12  2006/02/23 10:38:03  raasch
14! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
15! +z0 in argument list
16! WARNING: loops containing the MAX function are still not properly vectorized!
17!
18! Revision 1.1  1997/09/12 06:24:11  raasch
19! Initial revision
20!
21!
22! Description:
23! ------------
24! Diffusion term of the w-component
25!------------------------------------------------------------------------------!
26
27    PRIVATE
28    PUBLIC diffusion_w
29
30    INTERFACE diffusion_w
31       MODULE PROCEDURE diffusion_w
32       MODULE PROCEDURE diffusion_w_ij
33    END INTERFACE diffusion_w
34
35 CONTAINS
36
37
38!------------------------------------------------------------------------------!
39! Call for all grid points
40!------------------------------------------------------------------------------!
41    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
42                            w, z0 )
43
44       USE control_parameters
45       USE grid_variables
46       USE indices
47
48       IMPLICIT NONE
49
50       INTEGER ::  i, j, k
51       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
52                   kmyp_z, wsus, wsvs
53       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
54                   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 ::  km, u, v, w
58
59
60       DO  i = nxl, nxr
61          DO  j = nys, nyn
62             DO  k = nzb_w_outer(j,i)+1, nzt-1
63!
64!--             Interpolate eddy diffusivities on staggered gridpoints
65                kmxp_x = 0.25 * &
66                         ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
67                kmxm_x = 0.25 * &
68                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
69                kmxp_z = kmxp_x
70                kmxm_z = kmxm_x
71                kmyp_y = 0.25 * &
72                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
73                kmym_y = 0.25 * &
74                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
75                kmyp_z = kmyp_y
76                kmym_z = kmym_y
77!
78!--             Increase diffusion at the outflow boundary in case of
79!--             non-cyclic lateral boundaries. Damping is only needed for
80!--             velocity components parallel to the outflow boundary in
81!--             the direction normal to the outflow boundary.
82                IF ( bc_lr /= 'cyclic' )  THEN
83                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
84                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
85                ENDIF
86                IF ( bc_ns /= 'cyclic' )  THEN
87                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
88                   kmym_y = MAX( kmym_y, km_damp_y(j) )
89                ENDIF
90
91                tend(k,j,i) = tend(k,j,i)                                      &
92                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
93                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
94                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
95                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
96                      &   ) * ddx                                              &
97                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
98                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
99                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
100                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
101                      &   ) * ddy                                              &
102                      & + 2.0 * (                                              &
103                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
104                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
105                      &         ) * ddzu(k+1)
106             ENDDO
107
108!
109!--          Wall functions at all vertical walls, where necessary
110             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
111                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
112                   IF ( wall_w_x(j,i) /= 0.0 )  THEN
113                      wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
114                      wsus = -wsus * ABS( wsus )
115                   ELSE
116                      wsus = 0.0
117                   ENDIF
118                   IF ( wall_w_y(j,i) /= 0.0 )  THEN
119                      wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
120                      wsvs = -wsvs * ABS( wsvs )
121                   ELSE
122                      wsvs = 0.0
123                   ENDIF
124!
125!--                Interpolate eddy diffusivities on staggered gridpoints
126                   kmxp_x = 0.25 * &
127                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
128                   kmxm_x = 0.25 * &
129                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
130                   kmxp_z = kmxp_x
131                   kmxm_z = kmxm_x
132                   kmyp_y = 0.25 * &
133                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
134                   kmym_y = 0.25 * &
135                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
136                   kmyp_z = kmyp_y
137                   kmym_z = kmym_y
138!
139!--                Increase diffusion at the outflow boundary in case of
140!--                non-cyclic lateral boundaries. Damping is only needed for
141!--                velocity components parallel to the outflow boundary in
142!--                the direction normal to the outflow boundary.
143                   IF ( bc_lr /= 'cyclic' )  THEN
144                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
145                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
146                   ENDIF
147                   IF ( bc_ns /= 'cyclic' )  THEN
148                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
149                      kmym_y = MAX( kmym_y, km_damp_y(j) )
150                   ENDIF
151
152                   tend(k,j,i) = tend(k,j,i)                                   &
153                                 + (   fwxp(j,i) * (                           &
154                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
155                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
156                                                   )                           &
157                                     - fwxm(j,i) * (                           &
158                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
159                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
160                                                   )                           &
161                                     + wall_w_x(j,i) * wsus                    &
162                                   ) * ddx                                     &
163                                 + (   fwyp(j,i) * (                           &
164                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
165                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
166                                                   )                           &
167                                     - fwym(j,i) * (                           &
168                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
169                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
170                                                   )                           &
171                                     + wall_w_y(j,i) * wsvs                    &
172                                   ) * ddy                                     &
173                                 + 2.0 * (                                     &
174                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
175                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
176                                         ) * ddzu(k+1)
177                ENDDO
178             ENDIF
179
180          ENDDO
181       ENDDO
182
183    END SUBROUTINE diffusion_w
184
185
186!------------------------------------------------------------------------------!
187! Call for grid point i,j
188!------------------------------------------------------------------------------!
189    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
190                               tend, u, v, w, z0 )
191
192       USE control_parameters
193       USE grid_variables
194       USE indices
195
196       IMPLICIT NONE
197
198       INTEGER ::  i, j, k
199       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
200                   kmyp_z, wsus, wsvs
201       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
202                   km_damp_y(nys-1:nyn+1)
203       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
204       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
205       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
206
207
208       DO  k = nzb_w_outer(j,i)+1, nzt-1
209!
210!--       Interpolate eddy diffusivities on staggered gridpoints
211          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
212          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
213          kmxp_z = kmxp_x
214          kmxm_z = kmxm_x
215          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
216          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
217          kmyp_z = kmyp_y
218          kmym_z = kmym_y
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_lr /= 'cyclic' )  THEN
225             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
226             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
227          ENDIF
228          IF ( bc_ns /= 'cyclic' )  THEN
229             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
230             kmym_y = MAX( kmym_y, km_damp_y(j) )
231          ENDIF
232
233          tend(k,j,i) = tend(k,j,i)                                            &
234                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
235                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
236                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
237                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
238                      &   ) * ddx                                              &
239                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
240                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
241                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
242                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
243                      &   ) * ddy                                              &
244                      & + 2.0 * (                                              &
245                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
246                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
247                      &         ) * ddzu(k+1)
248       ENDDO
249
250!
251!--    Wall functions at all vertical walls, where necessary
252       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
253          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
254             IF ( wall_w_x(j,i) /= 0.0 )  THEN
255                wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
256                wsus = -wsus * ABS( wsus )
257             ELSE
258                wsus = 0.0
259             ENDIF
260             IF ( wall_w_y(j,i) /= 0.0 )  THEN
261                wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
262                wsvs = -wsvs * ABS( wsvs )
263             ELSE
264                wsvs = 0.0
265             ENDIF
266!
267!--          Interpolate eddy diffusivities on staggered gridpoints
268             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
269             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
270             kmxp_z = kmxp_x
271             kmxm_z = kmxm_x
272             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
273             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
274             kmyp_z = kmyp_y
275             kmym_z = kmym_y
276!
277!--          Increase diffusion at the outflow boundary in case of
278!--          non-cyclic lateral boundaries. Damping is only needed for
279!--          velocity components parallel to the outflow boundary in
280!--          the direction normal to the outflow boundary.
281             IF ( bc_lr /= 'cyclic' )  THEN
282                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
283                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
284             ENDIF
285             IF ( bc_ns /= 'cyclic' )  THEN
286                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
287                kmym_y = MAX( kmym_y, km_damp_y(j) )
288             ENDIF
289
290             tend(k,j,i) = tend(k,j,i)                                         &
291                                 + (   fwxp(j,i) * (                           &
292                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
293                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
294                                                   )                           &
295                                     - fwxm(j,i) * (                           &
296                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
297                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
298                                                   )                           &
299                                     + wall_w_x(j,i) * wsus                    &
300                                   ) * ddx                                     &
301                                 + (   fwyp(j,i) * (                           &
302                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
303                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
304                                                   )                           &
305                                     - fwym(j,i) * (                           &
306                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
307                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
308                                                   )                           &
309                                     + wall_w_y(j,i) * wsvs                    &
310                                   ) * ddy                                     &
311                                 + 2.0 * (                                     &
312                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
313                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
314                                         ) * ddzu(k+1)
315          ENDDO
316       ENDIF
317
318    END SUBROUTINE diffusion_w_ij
319
320 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.