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

Last change on this file since 70 was 57, checked in by raasch, 18 years ago

preliminary update of further changes, advec_particles is not running!

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