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

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

further checkin of preliminary changes

  • Property svn:keywords set to Id
File size: 15.5 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!
8! Former revisions:
9! -----------------
10! $Id: diffusion_v.f90 56 2007-03-08 13:57:07Z raasch $
11!
12! 20 2007-02-26 00:12:32Z raasch
13! Bugfix: ddzw dimensioned 1:nzt"+1"
14!
15! RCS Log replace by Id keyword, revision history cleaned up
16!
17! Revision 1.15  2006/02/23 10:36:00  raasch
18! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner
19! or nzb_diff_v, respectively, in vertical diffusion,
20! wall functions added for north and south walls, +z0 in argument list,
21! terms containing w(k-1,..) are removed from the Prandtl-layer equation
22! because they cause errors at the edges of topography
23! WARNING: loops containing the MAX function are still not properly vectorized!
24!
25! Revision 1.1  1997/09/12 06:24:01  raasch
26! Initial revision
27!
28!
29! Description:
30! ------------
31! Diffusion term of the v-component
32!------------------------------------------------------------------------------!
33
34    USE wall_fluxes_mod
35
36    PRIVATE
37    PUBLIC diffusion_v
38
39    INTERFACE diffusion_v
40       MODULE PROCEDURE diffusion_v
41       MODULE PROCEDURE diffusion_v_ij
42    END INTERFACE diffusion_v
43
44 CONTAINS
45
46
47!------------------------------------------------------------------------------!
48! Call for all grid points
49!------------------------------------------------------------------------------!
50    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
51
52       USE control_parameters
53       USE grid_variables
54       USE indices
55
56       IMPLICIT NONE
57
58       INTEGER ::  i, j, k
59       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
60       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
61       REAL    ::  z0(nys-1:nyn+1,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, z0 )
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    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
219       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
220       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
221       REAL, DIMENSION(:,:),   POINTER ::  vsws
222       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
223
224!
225!--    Compute horizontal diffusion
226       DO  k = nzb_v_outer(j,i)+1, nzt
227!
228!--       Interpolate eddy diffusivities on staggered gridpoints
229          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
230          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
231          kmxp_y = kmxp_x
232          kmxm_y = kmxm_x
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 ( bc_lr /= 'cyclic' )  THEN
239             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
240             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
241          ENDIF
242
243          tend(k,j,i) = tend(k,j,i)                                          &
244                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
245                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
246                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
247                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
248                      &   ) * ddx                                            &
249                      & + 2.0 * (                                            &
250                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
251                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
252                      &         ) * ddy2
253       ENDDO
254
255!
256!--    Wall functions at the left and right walls, respectively
257       IF ( wall_v(j,i) /= 0.0 )  THEN
258
259!
260!--       Calculate the horizontal momentum flux v'u'
261          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
262                            vsus, 0.0, 1.0, 0.0, 0.0 )
263
264          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
265             kmxp_x = 0.25 * &
266                      ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
267             kmxm_x = 0.25 * &
268                      ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
269             kmxp_y = kmxp_x
270             kmxm_y = kmxm_x
271!
272!--          Increase diffusion at the outflow boundary in case of
273!--          non-cyclic lateral boundaries. Damping is only needed for
274!--          velocity components parallel to the outflow boundary in
275!--          the direction normal to the outflow boundary.
276             IF ( bc_lr /= 'cyclic' )  THEN
277                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
278                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
279             ENDIF
280
281             tend(k,j,i) = tend(k,j,i)                                         &
282                                 + 2.0 * (                                     &
283                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
284                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
285                                         ) * ddy2                              &
286                                 + (   fxp(j,i) * (                            &
287                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
288                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
289                                                  )                            &
290                                     - fxm(j,i) * (                            &
291                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
292                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
293                                                  )                            &
294                                     + wall_v(j,i) * vsus(k)                   &
295                                   ) * ddx
296          ENDDO
297       ENDIF
298
299!
300!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
301!--    index k starts at nzb_v_inner+2.
302       DO  k = nzb_diff_v(j,i), nzt
303!
304!--       Interpolate eddy diffusivities on staggered gridpoints
305          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
306          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
307
308          tend(k,j,i) = tend(k,j,i)                                          &
309                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
310                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
311                      &            )                                         &
312                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
313                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
314                      &            )                                         &
315                      &   ) * ddzw(k)
316       ENDDO
317
318!
319!--    Vertical diffusion at the first grid point above the surface, if the
320!--    momentum flux at the bottom is given by the Prandtl law or if it is
321!--    prescribed by the user.
322!--    Difference quotient of the momentum flux is not formed over half of
323!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
324!--    other (LES) modell showed that the values of the momentum flux becomes
325!--    too large in this case.
326!--    The term containing w(k-1,..) (see above equation) is removed here
327!--    because the vertical velocity is assumed to be zero at the surface.
328       IF ( use_surface_fluxes )  THEN
329          k = nzb_v_inner(j,i)+1
330!
331!--       Interpolate eddy diffusivities on staggered gridpoints
332          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
333          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
334
335          tend(k,j,i) = tend(k,j,i)                                          &
336                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
337                      &   ) * ddzw(k)                                        &
338                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
339                      &   + vsws(j,i)                                        &
340                      &   ) * ddzw(k)
341       ENDIF
342
343    END SUBROUTINE diffusion_v_ij
344
345 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.