source: palm/trunk/SOURCE/diffusion_v.f90 @ 106

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

preliminary update of bugfixes and extensions for non-cyclic BCs

  • Property svn:keywords set to Id
File size: 17.6 KB
Line 
1 MODULE diffusion_v_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Momentumflux at top (vswst) included as boundary condition,
7! j loop is starting from nysv (needed for non-cyclic boundary conditions)
8!
9! Former revisions:
10! -----------------
11! $Id: diffusion_v.f90 106 2007-08-16 14:30:26Z raasch $
12!
13! 75 2007-03-22 09:54:05Z raasch
14! Wall functions now include diabatic conditions, call of routine wall_fluxes,
15! z0 removed from argument list, vynp eliminated
16!
17! 20 2007-02-26 00:12:32Z raasch
18! Bugfix: ddzw dimensioned 1:nzt"+1"
19!
20! RCS Log replace by Id keyword, revision history cleaned up
21!
22! Revision 1.15  2006/02/23 10:36:00  raasch
23! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner
24! or nzb_diff_v, respectively, in vertical diffusion,
25! wall functions added for north and south walls, +z0 in argument list,
26! terms containing w(k-1,..) are removed from the Prandtl-layer equation
27! because they cause errors at the edges of topography
28! WARNING: loops containing the MAX function are still not properly vectorized!
29!
30! Revision 1.1  1997/09/12 06:24:01  raasch
31! Initial revision
32!
33!
34! Description:
35! ------------
36! Diffusion term of the v-component
37!------------------------------------------------------------------------------!
38
39    USE wall_fluxes_mod
40
41    PRIVATE
42    PUBLIC diffusion_v
43
44    INTERFACE diffusion_v
45       MODULE PROCEDURE diffusion_v
46       MODULE PROCEDURE diffusion_v_ij
47    END INTERFACE diffusion_v
48
49 CONTAINS
50
51
52!------------------------------------------------------------------------------!
53! Call for all grid points
54!------------------------------------------------------------------------------!
55    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, &
56                            vswst, w )
57
58       USE control_parameters
59       USE grid_variables
60       USE indices
61
62       IMPLICIT NONE
63
64       INTEGER ::  i, j, k
65       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
66       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
67       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
68       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
69       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
70       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
71
72!
73!--    First calculate horizontal momentum flux v'u' at vertical walls,
74!--    if neccessary
75       IF ( topography /= 'flat' )  THEN
76          CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
77                            nzb_v_outer, wall_v )
78       ENDIF
79
80       DO  i = nxl, nxr
81          DO  j = nysv, nyn
82!
83!--          Compute horizontal diffusion
84             DO  k = nzb_v_outer(j,i)+1, nzt
85!
86!--             Interpolate eddy diffusivities on staggered gridpoints
87                kmxp_x = 0.25 * &
88                         ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
89                kmxm_x = 0.25 * &
90                         ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
91                kmxp_y = kmxp_x
92                kmxm_y = kmxm_x
93!
94!--             Increase diffusion at the outflow boundary in case of
95!--             non-cyclic lateral boundaries. Damping is only needed for
96!--             velocity components parallel to the outflow boundary in
97!--             the direction normal to the outflow boundary.
98                IF ( bc_lr /= 'cyclic' )  THEN
99                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
100                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
101                ENDIF
102
103                tend(k,j,i) = tend(k,j,i)                                    &
104                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
105                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
106                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
107                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
108                      &   ) * ddx                                            &
109                      & + 2.0 * (                                            &
110                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
111                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
112                      &         ) * ddy2
113             ENDDO
114
115!
116!--          Wall functions at the left and right walls, respectively
117             IF ( wall_v(j,i) /= 0.0 )  THEN
118
119                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
120                   kmxp_x = 0.25 * &
121                            ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
122                   kmxm_x = 0.25 * &
123                            ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
124                   kmxp_y = kmxp_x
125                   kmxm_y = kmxm_x
126!
127!--                Increase diffusion at the outflow boundary in case of
128!--                non-cyclic lateral boundaries. Damping is only needed for
129!--                velocity components parallel to the outflow boundary in
130!--                the direction normal to the outflow boundary.
131                   IF ( bc_lr /= 'cyclic' )  THEN
132                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
133                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
134                   ENDIF
135
136                   tend(k,j,i) = tend(k,j,i)                                   &
137                                 + 2.0 * (                                     &
138                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
139                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
140                                         ) * ddy2                              &
141                                 + (   fxp(j,i) * (                            &
142                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
143                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
144                                                  )                            &
145                                     - fxm(j,i) * (                            &
146                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
147                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
148                                                  )                            &
149                                     + wall_v(j,i) * vsus(k,j,i)               &
150                                   ) * ddx
151                ENDDO
152             ENDIF
153
154!
155!--          Compute vertical diffusion. In case of simulating a Prandtl
156!--          layer, index k starts at nzb_v_inner+2.
157             DO  k = nzb_diff_v(j,i), nzt_diff
158!
159!--             Interpolate eddy diffusivities on staggered gridpoints
160                kmzp = 0.25 * &
161                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
162                kmzm = 0.25 * &
163                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
164
165                tend(k,j,i) = tend(k,j,i)                                    &
166                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
167                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
168                      &            )                                         &
169                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
170                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
171                      &            )                                         &
172                      &   ) * ddzw(k)
173             ENDDO
174
175!
176!--          Vertical diffusion at the first grid point above the surface,
177!--          if the momentum flux at the bottom is given by the Prandtl law
178!--          or if it is prescribed by the user.
179!--          Difference quotient of the momentum flux is not formed over
180!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
181!--          comparison with other (LES) modell showed that the values of
182!--          the momentum flux becomes too large in this case.
183!--          The term containing w(k-1,..) (see above equation) is removed here
184!--          because the vertical velocity is assumed to be zero at the surface.
185             IF ( use_surface_fluxes )  THEN
186                k = nzb_v_inner(j,i)+1
187!
188!--             Interpolate eddy diffusivities on staggered gridpoints
189                kmzp = 0.25 * &
190                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
191                kmzm = 0.25 * &
192                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
193
194                tend(k,j,i) = tend(k,j,i)                                    &
195                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
196                      &   ) * ddzw(k)                                        &
197                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
198                      &   + vsws(j,i)                                        &
199                      &   ) * ddzw(k)
200             ENDIF
201
202!
203!--          Vertical diffusion at the first gridpoint below the top boundary,
204!--          if the momentum flux at the top is prescribed by the user
205             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
206                k = nzt
207!
208!--             Interpolate eddy diffusivities on staggered gridpoints
209                kmzp = 0.25 * &
210                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
211                kmzm = 0.25 * &
212                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
213
214                tend(k,j,i) = tend(k,j,i)                                    &
215                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
216                      &   ) * ddzw(k)                                        &
217                      & + ( -vswst(j,i)                                      &
218                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
219                      &   ) * ddzw(k)
220             ENDIF
221
222          ENDDO
223       ENDDO
224
225    END SUBROUTINE diffusion_v
226
227
228!------------------------------------------------------------------------------!
229! Call for grid point i,j
230!------------------------------------------------------------------------------!
231    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
232                               vsws, vswst, w )
233
234       USE control_parameters
235       USE grid_variables
236       USE indices
237
238       IMPLICIT NONE
239
240       INTEGER ::  i, j, k
241       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
242       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
243       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
244       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
245       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
246       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
247
248!
249!--    Compute horizontal diffusion
250       DO  k = nzb_v_outer(j,i)+1, nzt
251!
252!--       Interpolate eddy diffusivities on staggered gridpoints
253          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
254          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
255          kmxp_y = kmxp_x
256          kmxm_y = kmxm_x
257!
258!--       Increase diffusion at the outflow boundary in case of non-cyclic
259!--       lateral boundaries. Damping is only needed for velocity components
260!--       parallel to the outflow boundary in the direction normal to the
261!--       outflow boundary.
262          IF ( bc_lr /= 'cyclic' )  THEN
263             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
264             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
265          ENDIF
266
267          tend(k,j,i) = tend(k,j,i)                                          &
268                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
269                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
270                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
271                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
272                      &   ) * ddx                                            &
273                      & + 2.0 * (                                            &
274                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
275                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
276                      &         ) * ddy2
277       ENDDO
278
279!
280!--    Wall functions at the left and right walls, respectively
281       IF ( wall_v(j,i) /= 0.0 )  THEN
282
283!
284!--       Calculate the horizontal momentum flux v'u'
285          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
286                            vsus, 0.0, 1.0, 0.0, 0.0 )
287
288          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
289             kmxp_x = 0.25 * &
290                      ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
291             kmxm_x = 0.25 * &
292                      ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
293             kmxp_y = kmxp_x
294             kmxm_y = kmxm_x
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 ( bc_lr /= 'cyclic' )  THEN
301                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
302                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
303             ENDIF
304
305             tend(k,j,i) = tend(k,j,i)                                         &
306                                 + 2.0 * (                                     &
307                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
308                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
309                                         ) * ddy2                              &
310                                 + (   fxp(j,i) * (                            &
311                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
312                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
313                                                  )                            &
314                                     - fxm(j,i) * (                            &
315                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
316                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
317                                                  )                            &
318                                     + wall_v(j,i) * vsus(k)                   &
319                                   ) * ddx
320          ENDDO
321       ENDIF
322
323!
324!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
325!--    index k starts at nzb_v_inner+2.
326       DO  k = nzb_diff_v(j,i), nzt_diff
327!
328!--       Interpolate eddy diffusivities on staggered gridpoints
329          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
330          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
331
332          tend(k,j,i) = tend(k,j,i)                                          &
333                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
334                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
335                      &            )                                         &
336                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
337                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
338                      &            )                                         &
339                      &   ) * ddzw(k)
340       ENDDO
341
342!
343!--    Vertical diffusion at the first grid point above the surface, if the
344!--    momentum flux at the bottom is given by the Prandtl law or if it is
345!--    prescribed by the user.
346!--    Difference quotient of the momentum flux is not formed over half of
347!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
348!--    other (LES) modell showed that the values of the momentum flux becomes
349!--    too large in this case.
350!--    The term containing w(k-1,..) (see above equation) is removed here
351!--    because the vertical velocity is assumed to be zero at the surface.
352       IF ( use_surface_fluxes )  THEN
353          k = nzb_v_inner(j,i)+1
354!
355!--       Interpolate eddy diffusivities on staggered gridpoints
356          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
357          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
358
359          tend(k,j,i) = tend(k,j,i)                                          &
360                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
361                      &   ) * ddzw(k)                                        &
362                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
363                      &   + vsws(j,i)                                        &
364                      &   ) * ddzw(k)
365       ENDIF
366
367!
368!--    Vertical diffusion at the first gridpoint below the top boundary,
369!--    if the momentum flux at the top is prescribed by the user
370       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
371          k = nzt
372!
373!--       Interpolate eddy diffusivities on staggered gridpoints
374          kmzp = 0.25 * &
375                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
376          kmzm = 0.25 * &
377                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
378
379          tend(k,j,i) = tend(k,j,i)                                          &
380                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
381                      &   ) * ddzw(k)                                        &
382                      & + ( -vswst(j,i)                                      &
383                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
384                      &   ) * ddzw(k)
385       ENDIF
386
387    END SUBROUTINE diffusion_v_ij
388
389 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.