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

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

further checkin of preliminary changes

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