source: palm/tags/release-3.1c/SOURCE/diffusion_u.f90 @ 366

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

comments prepared for 3.1c

  • Property svn:keywords set to Id
File size: 15.0 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_u.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:35:35  raasch
18! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
19! or nzb_diff_u, 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:23:51  raasch
26! Initial revision
27!
28!
29! Description:
30! ------------
31! Diffusion term of the u-component
32!------------------------------------------------------------------------------!
33
34    PRIVATE
35    PUBLIC diffusion_u
36
37    INTERFACE diffusion_u
38       MODULE PROCEDURE diffusion_u
39       MODULE PROCEDURE diffusion_u_ij
40    END INTERFACE diffusion_u
41
42 CONTAINS
43
44
45!------------------------------------------------------------------------------!
46! Call for all grid points
47!------------------------------------------------------------------------------!
48    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, 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    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
58       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+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 ::  usws
62       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
63
64       DO  i = nxl, nxr+uxrp
65          DO  j = nys,nyn
66!
67!--          Compute horizontal diffusion
68             DO  k = nzb_u_outer(j,i)+1, nzt
69!
70!--             Interpolate eddy diffusivities on staggered gridpoints
71                kmyp_x = 0.25 * &
72                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
73                kmym_x = 0.25 * &
74                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
75                kmyp_y = kmyp_x
76                kmym_y = kmym_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_ns /= 'cyclic' )  THEN
83                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
84                   kmym_y = MAX( kmym_y, km_damp_y(j) )
85                ENDIF
86
87                tend(k,j,i) = tend(k,j,i)                                    &
88                      & + 2.0 * (                                            &
89                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
90                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
91                      &         ) * ddx2                                     &
92                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
93                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
94                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
95                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
96                      &   ) * ddy
97             ENDDO
98
99!
100!--          Wall functions at the north and south walls, respectively
101             IF ( wall_u(j,i) /= 0.0 )  THEN
102                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
103                   usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
104                   usvs   = -usvs * ABS( usvs )
105                   kmyp_x = 0.25 * &
106                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
107                   kmym_x = 0.25 * &
108                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
109                   kmyp_y = kmyp_x
110                   kmym_y = kmym_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_ns /= 'cyclic' )  THEN
117                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
118                      kmym_y = MAX( kmym_y, km_damp_y(j) )
119                   ENDIF
120
121                   tend(k,j,i) = tend(k,j,i)                                   &
122                                 + 2.0 * (                                     &
123                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
124                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
125                                         ) * ddx2                              &
126                                 + (   fyp(j,i) * (                            &
127                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
128                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
129                                                  )                            &
130                                     - fym(j,i) * (                            &
131                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
132                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
133                                                  )                            &
134                                     + wall_u(j,i) * usvs                      &
135                                   ) * ddy
136                ENDDO
137             ENDIF
138
139!
140!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
141!--          index k starts at nzb_u_inner+2.
142             DO  k = nzb_diff_u(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,i-1)+km(k+1,j,i-1) )
147                kmzm = 0.25 * &
148                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
149
150                tend(k,j,i) = tend(k,j,i)                                    &
151                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
152                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
153                      &            )                                         &
154                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
155                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
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 or
163!--          if it is prescribed by the user.
164!--          Difference quotient of the momentum flux is not formed over half
165!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
166!--          with other (LES) modell showed that the values of the momentum
167!--          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_u_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,i-1)+km(k+1,j,i-1) )
176                kmzm = 0.25 * &
177                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
178
179                tend(k,j,i) = tend(k,j,i)                                    &
180                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
181                      &   ) * ddzw(k)                                        &
182                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
183                      &   + usws(j,i)                                        &
184                      &   ) * ddzw(k)
185             ENDIF
186
187          ENDDO
188       ENDDO
189
190    END SUBROUTINE diffusion_u
191
192
193!------------------------------------------------------------------------------!
194! Call for grid point i,j
195!------------------------------------------------------------------------------!
196    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
197                               v, 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    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
207       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+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 ::  usws
211       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
212
213!
214!--    Compute horizontal diffusion
215       DO  k = nzb_u_outer(j,i)+1, nzt
216!
217!--       Interpolate eddy diffusivities on staggered gridpoints
218          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
219          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
220          kmyp_y = kmyp_x
221          kmym_y = kmym_x
222
223!
224!--       Increase diffusion at the outflow boundary in case of non-cyclic
225!--       lateral boundaries. Damping is only needed for velocity components
226!--       parallel to the outflow boundary in the direction normal to the
227!--       outflow boundary.
228          IF ( bc_ns /= 'cyclic' )  THEN
229             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
230             kmym_y = MAX( kmym_y, km_damp_y(j) )
231          ENDIF
232
233          tend(k,j,i) = tend(k,j,i)                                          &
234                      & + 2.0 * (                                            &
235                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
236                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
237                      &         ) * ddx2                                     &
238                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
239                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
240                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
241                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
242                      &   ) * ddy
243       ENDDO
244
245!
246!--    Wall functions at the north and south walls, respectively
247       IF ( wall_u(j,i) .NE. 0.0 )  THEN
248          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
249             usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
250             usvs   = -usvs * ABS( usvs )
251             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
252             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
253             kmyp_y = kmyp_x
254             kmym_y = kmym_x
255!
256!--          Increase diffusion at the outflow boundary in case of
257!--          non-cyclic lateral boundaries. Damping is only needed for
258!--          velocity components parallel to the outflow boundary in
259!--          the direction normal to the outflow boundary.
260             IF ( bc_ns /= 'cyclic' )  THEN
261                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
262                kmym_y = MAX( kmym_y, km_damp_y(j) )
263             ENDIF
264
265             tend(k,j,i) = tend(k,j,i)                                         &
266                                 + 2.0 * (                                     &
267                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
268                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
269                                         ) * ddx2                              &
270                                 + (   fyp(j,i) * (                            &
271                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
272                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
273                                                  )                            &
274                                     - fym(j,i) * (                            &
275                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
276                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
277                                                  )                            &
278                                     + wall_u(j,i) * usvs                      &
279                                   ) * ddy
280          ENDDO
281       ENDIF
282
283!
284!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
285!--    index k starts at nzb_u_inner+2.
286       DO  k = nzb_diff_u(j,i), nzt
287!
288!--       Interpolate eddy diffusivities on staggered gridpoints
289          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
290          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
291
292          tend(k,j,i) = tend(k,j,i)                                          &
293                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
294                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
295                      &            )                                         &
296                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
297                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
298                      &            )                                         &
299                      &   ) * ddzw(k)
300       ENDDO
301
302!
303!--    Vertical diffusion at the first grid point above the surface, if the
304!--    momentum flux at the bottom is given by the Prandtl law or if it is
305!--    prescribed by the user.
306!--    Difference quotient of the momentum flux is not formed over half of
307!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
308!--    other (LES) modell showed that the values of the momentum flux becomes
309!--    too large in this case.
310!--    The term containing w(k-1,..) (see above equation) is removed here
311!--    because the vertical velocity is assumed to be zero at the surface.
312       IF ( use_surface_fluxes )  THEN
313          k = nzb_u_inner(j,i)+1
314!
315!--       Interpolate eddy diffusivities on staggered gridpoints
316          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
317          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
318
319          tend(k,j,i) = tend(k,j,i)                                          &
320                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
321                      &   ) * ddzw(k)                                        &
322                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
323                      &   + usws(j,i)                                        &
324                      &   ) * ddzw(k)
325       ENDIF
326
327    END SUBROUTINE diffusion_u_ij
328
329 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.