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

Last change on this file since 686 was 668, checked in by suehring, 13 years ago

last commit documented

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