source: palm/trunk/SOURCE/diffusion_u.f90 @ 635

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

typo in file headers removed

  • Property svn:keywords set to Id
File size: 17.7 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_u.f90 484 2010-02-05 07:36:54Z raasch $
11!
12! 366 2009-08-25 08:06:27Z raasch
13! bc_ns replaced by bc_ns_cyc
14!
15! 106 2007-08-16 14:30:26Z raasch
16! Momentumflux at top (uswst) included as boundary condition,
17! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
18!
19! 75 2007-03-22 09:54:05Z raasch
20! Wall functions now include diabatic conditions, call of routine wall_fluxes,
21! z0 removed from argument list, uxrp eliminated
22!
23! 20 2007-02-26 00:12:32Z raasch
24! Bugfix: ddzw dimensioned 1:nzt"+1"
25!
26! RCS Log replace by Id keyword, revision history cleaned up
27!
28! Revision 1.15  2006/02/23 10:35:35  raasch
29! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
30! or nzb_diff_u, respectively, in vertical diffusion,
31! wall functions added for north and south walls, +z0 in argument list,
32! terms containing w(k-1,..) are removed from the Prandtl-layer equation
33! because they cause errors at the edges of topography
34! WARNING: loops containing the MAX function are still not properly vectorized!
35!
36! Revision 1.1  1997/09/12 06:23:51  raasch
37! Initial revision
38!
39!
40! Description:
41! ------------
42! Diffusion term of the u-component
43! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
44!        and slows down the speed on NEC about 5-10%
45!------------------------------------------------------------------------------!
46
47    USE wall_fluxes_mod
48
49    PRIVATE
50    PUBLIC diffusion_u
51
52    INTERFACE diffusion_u
53       MODULE PROCEDURE diffusion_u
54       MODULE PROCEDURE diffusion_u_ij
55    END INTERFACE diffusion_u
56
57 CONTAINS
58
59
60!------------------------------------------------------------------------------!
61! Call for all grid points
62!------------------------------------------------------------------------------!
63    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
64                            v, w )
65
66       USE control_parameters
67       USE grid_variables
68       USE indices
69
70       IMPLICIT NONE
71
72       INTEGER ::  i, j, k
73       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
74       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
75       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
76       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
77       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
78       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
79
80!
81!--    First calculate horizontal momentum flux u'v' at vertical walls,
82!--    if neccessary
83       IF ( topography /= 'flat' )  THEN
84          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
85                            nzb_u_outer, wall_u )
86       ENDIF
87
88       DO  i = nxlu, nxr
89          DO  j = nys,nyn
90!
91!--          Compute horizontal diffusion
92             DO  k = nzb_u_outer(j,i)+1, nzt
93!
94!--             Interpolate eddy diffusivities on staggered gridpoints
95                kmyp_x = 0.25 * &
96                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
97                kmym_x = 0.25 * &
98                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
99                kmyp_y = kmyp_x
100                kmym_y = kmym_x
101!
102!--             Increase diffusion at the outflow boundary in case of
103!--             non-cyclic lateral boundaries. Damping is only needed for
104!--             velocity components parallel to the outflow boundary in
105!--             the direction normal to the outflow boundary.
106                IF ( .NOT. bc_ns_cyc )  THEN
107                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
108                   kmym_y = MAX( kmym_y, km_damp_y(j) )
109                ENDIF
110
111                tend(k,j,i) = tend(k,j,i)                                    &
112                      & + 2.0 * (                                            &
113                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
114                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
115                      &         ) * ddx2                                     &
116                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
117                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
118                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
119                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
120                      &   ) * ddy
121             ENDDO
122
123!
124!--          Wall functions at the north and south walls, respectively
125             IF ( wall_u(j,i) /= 0.0 )  THEN
126
127                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
128                   kmyp_x = 0.25 * &
129                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
130                   kmym_x = 0.25 * &
131                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
132                   kmyp_y = kmyp_x
133                   kmym_y = kmym_x
134!
135!--                Increase diffusion at the outflow boundary in case of
136!--                non-cyclic lateral boundaries. Damping is only needed for
137!--                velocity components parallel to the outflow boundary in
138!--                the direction normal to the outflow boundary.
139                   IF ( .NOT. bc_ns_cyc )  THEN
140                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
141                      kmym_y = MAX( kmym_y, km_damp_y(j) )
142                   ENDIF
143
144                   tend(k,j,i) = tend(k,j,i)                                   &
145                                 + 2.0 * (                                     &
146                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
147                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
148                                         ) * ddx2                              &
149                                 + (   fyp(j,i) * (                            &
150                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
151                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
152                                                  )                            &
153                                     - fym(j,i) * (                            &
154                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
155                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
156                                                  )                            &
157                                     + wall_u(j,i) * usvs(k,j,i)               &
158                                   ) * ddy
159                ENDDO
160             ENDIF
161
162!
163!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
164!--          index k starts at nzb_u_inner+2.
165             DO  k = nzb_diff_u(j,i), nzt_diff
166!
167!--             Interpolate eddy diffusivities on staggered gridpoints
168                kmzp = 0.25 * &
169                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
170                kmzm = 0.25 * &
171                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
172
173                tend(k,j,i) = tend(k,j,i)                                    &
174                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
175                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
176                      &            )                                         &
177                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
178                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
179                      &            )                                         &
180                      &   ) * ddzw(k)
181             ENDDO
182
183!
184!--          Vertical diffusion at the first grid point above the surface,
185!--          if the momentum flux at the bottom is given by the Prandtl law or
186!--          if it is prescribed by the user.
187!--          Difference quotient of the momentum flux is not formed over half
188!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
189!--          with other (LES) modell showed that the values of the momentum
190!--          flux becomes too large in this case.
191!--          The term containing w(k-1,..) (see above equation) is removed here
192!--          because the vertical velocity is assumed to be zero at the surface.
193             IF ( use_surface_fluxes )  THEN
194                k = nzb_u_inner(j,i)+1
195!
196!--             Interpolate eddy diffusivities on staggered gridpoints
197                kmzp = 0.25 * &
198                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
199                kmzm = 0.25 * &
200                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
201
202                tend(k,j,i) = tend(k,j,i)                                    &
203                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
204                      &   ) * ddzw(k)                                        &
205                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
206                      &   + usws(j,i)                                        &
207                      &   ) * ddzw(k)
208             ENDIF
209
210!
211!--          Vertical diffusion at the first gridpoint below the top boundary,
212!--          if the momentum flux at the top is prescribed by the user
213             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
214                k = nzt
215!
216!--             Interpolate eddy diffusivities on staggered gridpoints
217                kmzp = 0.25 * &
218                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
219                kmzm = 0.25 * &
220                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
221
222                tend(k,j,i) = tend(k,j,i)                                    &
223                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
224                      &   ) * ddzw(k)                                        &
225                      & + ( -uswst(j,i)                                      &
226                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
227                      &   ) * ddzw(k)
228             ENDIF
229
230          ENDDO
231       ENDDO
232
233    END SUBROUTINE diffusion_u
234
235
236!------------------------------------------------------------------------------!
237! Call for grid point i,j
238!------------------------------------------------------------------------------!
239    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
240                               uswst, v, w )
241
242       USE control_parameters
243       USE grid_variables
244       USE indices
245
246       IMPLICIT NONE
247
248       INTEGER ::  i, j, k
249       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
250       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
251       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
252       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
253       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
254       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
255
256!
257!--    Compute horizontal diffusion
258       DO  k = nzb_u_outer(j,i)+1, nzt
259!
260!--       Interpolate eddy diffusivities on staggered gridpoints
261          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
262          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
263          kmyp_y = kmyp_x
264          kmym_y = kmym_x
265
266!
267!--       Increase diffusion at the outflow boundary in case of non-cyclic
268!--       lateral boundaries. Damping is only needed for velocity components
269!--       parallel to the outflow boundary in the direction normal to the
270!--       outflow boundary.
271          IF ( .NOT. bc_ns_cyc )  THEN
272             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
273             kmym_y = MAX( kmym_y, km_damp_y(j) )
274          ENDIF
275
276          tend(k,j,i) = tend(k,j,i)                                          &
277                      & + 2.0 * (                                            &
278                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
279                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
280                      &         ) * ddx2                                     &
281                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
282                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
283                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
284                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
285                      &   ) * ddy
286       ENDDO
287
288!
289!--    Wall functions at the north and south walls, respectively
290       IF ( wall_u(j,i) .NE. 0.0 )  THEN
291
292!
293!--       Calculate the horizontal momentum flux u'v'
294          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
295                            usvs, 1.0, 0.0, 0.0, 0.0 )
296
297          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
298             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
299             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
300             kmyp_y = kmyp_x
301             kmym_y = kmym_x
302!
303!--          Increase diffusion at the outflow boundary in case of
304!--          non-cyclic lateral boundaries. Damping is only needed for
305!--          velocity components parallel to the outflow boundary in
306!--          the direction normal to the outflow boundary.
307             IF ( .NOT. bc_ns_cyc )  THEN
308                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
309                kmym_y = MAX( kmym_y, km_damp_y(j) )
310             ENDIF
311
312             tend(k,j,i) = tend(k,j,i)                                         &
313                                 + 2.0 * (                                     &
314                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
315                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
316                                         ) * ddx2                              &
317                                 + (   fyp(j,i) * (                            &
318                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
319                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
320                                                  )                            &
321                                     - fym(j,i) * (                            &
322                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
323                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
324                                                  )                            &
325                                     + wall_u(j,i) * usvs(k)                   &
326                                   ) * ddy
327          ENDDO
328       ENDIF
329
330!
331!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
332!--    index k starts at nzb_u_inner+2.
333       DO  k = nzb_diff_u(j,i), nzt_diff
334!
335!--       Interpolate eddy diffusivities on staggered gridpoints
336          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
337          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
338
339          tend(k,j,i) = tend(k,j,i)                                          &
340                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
341                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
342                      &            )                                         &
343                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
344                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
345                      &            )                                         &
346                      &   ) * ddzw(k)
347       ENDDO
348
349!
350!--    Vertical diffusion at the first grid point above the surface, if the
351!--    momentum flux at the bottom is given by the Prandtl law or if it is
352!--    prescribed by the user.
353!--    Difference quotient of the momentum flux is not formed over half of
354!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
355!--    other (LES) modell showed that the values of the momentum flux becomes
356!--    too large in this case.
357!--    The term containing w(k-1,..) (see above equation) is removed here
358!--    because the vertical velocity is assumed to be zero at the surface.
359       IF ( use_surface_fluxes )  THEN
360          k = nzb_u_inner(j,i)+1
361!
362!--       Interpolate eddy diffusivities on staggered gridpoints
363          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
364          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
365
366          tend(k,j,i) = tend(k,j,i)                                          &
367                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
368                      &   ) * ddzw(k)                                        &
369                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
370                      &   + usws(j,i)                                        &
371                      &   ) * ddzw(k)
372       ENDIF
373
374!
375!--    Vertical diffusion at the first gridpoint below the top boundary,
376!--    if the momentum flux at the top is prescribed by the user
377       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
378          k = nzt
379!
380!--       Interpolate eddy diffusivities on staggered gridpoints
381          kmzp = 0.25 * &
382                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
383          kmzm = 0.25 * &
384                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
385
386          tend(k,j,i) = tend(k,j,i)                                          &
387                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
388                      &   ) * ddzw(k)                                        &
389                      & + ( -uswst(j,i)                                      &
390                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
391                      &   ) * ddzw(k)
392       ENDIF
393
394    END SUBROUTINE diffusion_u_ij
395
396 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.