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

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

Initial repository layout and content

File size: 16.4 KB
Line 
1 MODULE diffusion_w_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: diffusion_w.f90,v $
11! Revision 1.12  2006/02/23 10:38:03  raasch
12! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
13! +z0 in argument list
14! WARNING: loops containing the MAX function are still not properly vectorized!
15!
16! Revision 1.11  2005/03/26 20:11:05  raasch
17! Additional damping layer at the outflow in case of non-cyclic lateral
18! boundaries, additional arguments km_damp_x, km_damp_y
19!
20! Revision 1.10  2004/01/30 10:22:32  raasch
21! Scalar lower k index nzb replaced by 2d-array nzb_2d
22!
23! Revision 1.9  2003/03/12 16:26:49  raasch
24! Full code replaced in the call for all gridpoints instead of calling the
25! _ij version (required by NEC, because otherwise no vectorization)
26!
27! Revision 1.8  2002/06/11 12:54:16  raasch
28! Former subroutine changed to a module which allows to be called for all grid
29! points of a single vertical column with index i,j or for all grid points by
30! using function overloading.
31!
32! Revision 1.7  2001/03/30 07:12:35  raasch
33! Translation of remaining German identifiers (variables, subroutines, etc.)
34!
35! Revision 1.6  2001/01/22 06:20:12  raasch
36! Module test_variables removed
37!
38! Revision 1.5  2000/07/03 12:58:40  raasch
39! dummy arguments, whose corresponding actual arguments are pointers,
40! are now also defined as pointers,
41! all comments translated into English
42!
43! Revision 1.4  1998/07/06 12:12:24  raasch
44! + USE test_variables
45!
46! Revision 1.3  1997/09/12 07:24:58  raasch
47! Leerzeilen mussten entfernt werden
48!
49! Revision 1.2  1997/09/12 06:44:00  raasch
50! HP-Compiler erfordert & am Beginn von Fortsetzungszeilen
51!
52! Revision 1.1  1997/09/12 06:24:11  raasch
53! Initial revision
54!
55!
56! Description:
57! ------------
58! Diffusion term of the w-component
59!------------------------------------------------------------------------------!
60
61    PRIVATE
62    PUBLIC diffusion_w
63
64    INTERFACE diffusion_w
65       MODULE PROCEDURE diffusion_w
66       MODULE PROCEDURE diffusion_w_ij
67    END INTERFACE diffusion_w
68
69 CONTAINS
70
71
72!------------------------------------------------------------------------------!
73! Call for all grid points
74!------------------------------------------------------------------------------!
75    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
76                            w, z0 )
77
78       USE control_parameters
79       USE grid_variables
80       USE indices
81
82       IMPLICIT NONE
83
84       INTEGER ::  i, j, k
85       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
86                   kmyp_z, wsus, wsvs
87       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
88                   km_damp_y(nys-1:nyn+1)
89       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
90       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
91       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
92
93
94       DO  i = nxl, nxr
95          DO  j = nys, nyn
96             DO  k = nzb_w_outer(j,i)+1, nzt-1
97!
98!--             Interpolate eddy diffusivities on staggered gridpoints
99                kmxp_x = 0.25 * &
100                         ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
101                kmxm_x = 0.25 * &
102                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
103                kmxp_z = kmxp_x
104                kmxm_z = kmxm_x
105                kmyp_y = 0.25 * &
106                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
107                kmym_y = 0.25 * &
108                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
109                kmyp_z = kmyp_y
110                kmym_z = kmym_y
111!
112!--             Increase diffusion at the outflow boundary in case of
113!--             non-cyclic lateral boundaries. Damping is only needed for
114!--             velocity components parallel to the outflow boundary in
115!--             the direction normal to the outflow boundary.
116                IF ( bc_lr /= 'cyclic' )  THEN
117                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
118                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
119                ENDIF
120                IF ( bc_ns /= 'cyclic' )  THEN
121                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
122                   kmym_y = MAX( kmym_y, km_damp_y(j) )
123                ENDIF
124
125                tend(k,j,i) = tend(k,j,i)                                      &
126                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
127                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
128                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
129                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
130                      &   ) * ddx                                              &
131                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
132                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
133                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
134                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
135                      &   ) * ddy                                              &
136                      & + 2.0 * (                                              &
137                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
138                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
139                      &         ) * ddzu(k+1)
140             ENDDO
141
142!
143!--          Wall functions at all vertical walls, where necessary
144             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
145                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
146                   IF ( wall_w_x(j,i) /= 0.0 )  THEN
147                      wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
148                      wsus = -wsus * ABS( wsus )
149                   ELSE
150                      wsus = 0.0
151                   ENDIF
152                   IF ( wall_w_y(j,i) /= 0.0 )  THEN
153                      wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
154                      wsvs = -wsvs * ABS( wsvs )
155                   ELSE
156                      wsvs = 0.0
157                   ENDIF
158!
159!--                Interpolate eddy diffusivities on staggered gridpoints
160                   kmxp_x = 0.25 * &
161                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
162                   kmxm_x = 0.25 * &
163                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
164                   kmxp_z = kmxp_x
165                   kmxm_z = kmxm_x
166                   kmyp_y = 0.25 * &
167                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
168                   kmym_y = 0.25 * &
169                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
170                   kmyp_z = kmyp_y
171                   kmym_z = kmym_y
172!
173!--                Increase diffusion at the outflow boundary in case of
174!--                non-cyclic lateral boundaries. Damping is only needed for
175!--                velocity components parallel to the outflow boundary in
176!--                the direction normal to the outflow boundary.
177                   IF ( bc_lr /= 'cyclic' )  THEN
178                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
179                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
180                   ENDIF
181                   IF ( bc_ns /= 'cyclic' )  THEN
182                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
183                      kmym_y = MAX( kmym_y, km_damp_y(j) )
184                   ENDIF
185
186                   tend(k,j,i) = tend(k,j,i)                                   &
187                                 + (   fwxp(j,i) * (                           &
188                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
189                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
190                                                   )                           &
191                                     - fwxm(j,i) * (                           &
192                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
193                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
194                                                   )                           &
195                                     + wall_w_x(j,i) * wsus                    &
196                                   ) * ddx                                     &
197                                 + (   fwyp(j,i) * (                           &
198                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
199                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
200                                                   )                           &
201                                     - fwym(j,i) * (                           &
202                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
203                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
204                                                   )                           &
205                                     + wall_w_y(j,i) * wsvs                    &
206                                   ) * ddy                                     &
207                                 + 2.0 * (                                     &
208                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
209                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
210                                         ) * ddzu(k+1)
211                ENDDO
212             ENDIF
213
214          ENDDO
215       ENDDO
216
217    END SUBROUTINE diffusion_w
218
219
220!------------------------------------------------------------------------------!
221! Call for grid point i,j
222!------------------------------------------------------------------------------!
223    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
224                               tend, u, v, w, z0 )
225
226       USE control_parameters
227       USE grid_variables
228       USE indices
229
230       IMPLICIT NONE
231
232       INTEGER ::  i, j, k
233       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
234                   kmyp_z, wsus, wsvs
235       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
236                   km_damp_y(nys-1:nyn+1)
237       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
238       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
239       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
240
241
242       DO  k = nzb_w_outer(j,i)+1, nzt-1
243!
244!--       Interpolate eddy diffusivities on staggered gridpoints
245          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
246          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
247          kmxp_z = kmxp_x
248          kmxm_z = kmxm_x
249          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
250          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
251          kmyp_z = kmyp_y
252          kmym_z = kmym_y
253!
254!--       Increase diffusion at the outflow boundary in case of non-cyclic
255!--       lateral boundaries. Damping is only needed for velocity components
256!--       parallel to the outflow boundary in the direction normal to the
257!--       outflow boundary.
258          IF ( bc_lr /= 'cyclic' )  THEN
259             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
260             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
261          ENDIF
262          IF ( bc_ns /= 'cyclic' )  THEN
263             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
264             kmym_y = MAX( kmym_y, km_damp_y(j) )
265          ENDIF
266
267          tend(k,j,i) = tend(k,j,i)                                            &
268                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
269                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
270                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
271                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
272                      &   ) * ddx                                              &
273                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
274                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
275                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
276                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
277                      &   ) * ddy                                              &
278                      & + 2.0 * (                                              &
279                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
280                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
281                      &         ) * ddzu(k+1)
282       ENDDO
283
284!
285!--    Wall functions at all vertical walls, where necessary
286       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
287          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
288             IF ( wall_w_x(j,i) /= 0.0 )  THEN
289                wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
290                wsus = -wsus * ABS( wsus )
291             ELSE
292                wsus = 0.0
293             ENDIF
294             IF ( wall_w_y(j,i) /= 0.0 )  THEN
295                wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
296                wsvs = -wsvs * ABS( wsvs )
297             ELSE
298                wsvs = 0.0
299             ENDIF
300!
301!--          Interpolate eddy diffusivities on staggered gridpoints
302             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
303             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
304             kmxp_z = kmxp_x
305             kmxm_z = kmxm_x
306             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
307             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
308             kmyp_z = kmyp_y
309             kmym_z = kmym_y
310!
311!--          Increase diffusion at the outflow boundary in case of
312!--          non-cyclic lateral boundaries. Damping is only needed for
313!--          velocity components parallel to the outflow boundary in
314!--          the direction normal to the outflow boundary.
315             IF ( bc_lr /= 'cyclic' )  THEN
316                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
317                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
318             ENDIF
319             IF ( bc_ns /= 'cyclic' )  THEN
320                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
321                kmym_y = MAX( kmym_y, km_damp_y(j) )
322             ENDIF
323
324             tend(k,j,i) = tend(k,j,i)                                         &
325                                 + (   fwxp(j,i) * (                           &
326                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
327                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
328                                                   )                           &
329                                     - fwxm(j,i) * (                           &
330                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
331                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
332                                                   )                           &
333                                     + wall_w_x(j,i) * wsus                    &
334                                   ) * ddx                                     &
335                                 + (   fwyp(j,i) * (                           &
336                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
337                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
338                                                   )                           &
339                                     - fwym(j,i) * (                           &
340                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
341                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
342                                                   )                           &
343                                     + wall_w_y(j,i) * wsvs                    &
344                                   ) * ddy                                     &
345                                 + 2.0 * (                                     &
346                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
347                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
348                                         ) * ddzu(k+1)
349          ENDDO
350       ENDIF
351
352    END SUBROUTINE diffusion_w_ij
353
354 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.