source: palm/trunk/SOURCE/diffusion_w.f90 @ 39

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

comments prepared for 3.1c

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