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

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

comments prepared for 3.1c

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