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

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

RCS Log replace by Id keyword, revision history cleaned up

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