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

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

typo in file headers removed

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