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

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

typo in file headers removed

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