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

Last change on this file since 102 was 102, checked in by raasch, 17 years ago

preliminary version for coupled runs

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