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

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

preliminary version, several changes to be explained later

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