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

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

Initial repository layout and content

File size: 16.7 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: diffusion_u.f90,v $
11! Revision 1.15  2006/02/23 10:35:35  raasch
12! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
13! or nzb_diff_u, respectively, in vertical diffusion,
14! wall functions added for north and south walls, +z0 in argument list,
15! terms containing w(k-1,..) are removed from the Prandtl-layer equation
16! because they cause errors at the edges of topography
17! WARNING: loops containing the MAX function are still not properly vectorized!
18!
19! Revision 1.14  2005/03/26 20:08:21  raasch
20! Extension of horizontal loop upper bounds for non-cyclic boundary conditions,
21! additional damping layer at the outflow in case of non-cyclic lateral
22! boundaries, additional argument km_damp_y
23!
24! Revision 1.13  2004/01/30 10:21:32  raasch
25! Scalar lower k index nzb replaced by 2d-array nzb_2d
26!
27! Revision 1.12  2003/03/12 16:26:14  raasch
28! Full code replaced in the call for all gridpoints instead of calling the
29! _ij version (required by NEC, because otherwise no vectorization)
30!
31! Revision 1.11  2002/06/11 12:53:12  raasch
32! Former subroutine changed to a module which allows to be called for all grid
33! points of a single vertical column with index i,j or for all grid points by
34! using function overloading.
35!
36! Revision 1.10  2001/08/21 08:25:55  raasch
37! Special treatment at k=1 generally if momentum flux is prescribed (not only in
38! case of a Prandtl layer)
39!
40! Revision 1.9  2001/03/30 07:12:02  raasch
41! Translation of remaining German identifiers (variables, subroutines, etc.)
42!
43! Revision 1.8  2001/01/22 06:14:00  raasch
44! Module test_variables removed
45!
46! Revision 1.7  2000/07/03 12:57:46  raasch
47! dummy arguments, whose corresponding actual arguments are pointers,
48! are now also defined as pointers,
49! all comments translated into English
50!
51! Revision 1.6  2000/03/14 13:55:33  schroeter
52! vertical impulse fluxes are computed over dz at the
53! top of the Prandtl-layer
54!
55! Revision 1.5  98/07/06  12:11:42  12:11:42  raasch (Siegfried Raasch)
56! + USE test_variables
57!
58! Revision 1.4  1998/01/23 09:56:44  raasch
59! Einbau der Prandtl-Schicht
60!
61! Revision 1.3  1997/09/12 07:22:50  raasch
62! Leerzeilen mussten entfernt werden
63!
64! Revision 1.2  1997/09/12 06:43:27  raasch
65! HP-Compiler erfordert & am Beginn von Fortsetzungszeilen
66!
67! Revision 1.1  1997/09/12 06:23:51  raasch
68! Initial revision
69!
70!
71! Description:
72! ------------
73! Diffusion term of the u-component
74!------------------------------------------------------------------------------!
75
76    PRIVATE
77    PUBLIC diffusion_u
78
79    INTERFACE diffusion_u
80       MODULE PROCEDURE diffusion_u
81       MODULE PROCEDURE diffusion_u_ij
82    END INTERFACE diffusion_u
83
84 CONTAINS
85
86
87!------------------------------------------------------------------------------!
88! Call for all grid points
89!------------------------------------------------------------------------------!
90    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
91
92       USE control_parameters
93       USE grid_variables
94       USE indices
95
96       IMPLICIT NONE
97
98       INTEGER ::  i, j, k
99       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
100       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_y(nys-1:nyn+1)
101       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
102       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
103       REAL, DIMENSION(:,:),   POINTER ::  usws
104       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
105
106       DO  i = nxl, nxr+uxrp
107          DO  j = nys,nyn
108!
109!--          Compute horizontal diffusion
110             DO  k = nzb_u_outer(j,i)+1, nzt
111!
112!--             Interpolate eddy diffusivities on staggered gridpoints
113                kmyp_x = 0.25 * &
114                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
115                kmym_x = 0.25 * &
116                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
117                kmyp_y = kmyp_x
118                kmym_y = kmym_x
119!
120!--             Increase diffusion at the outflow boundary in case of
121!--             non-cyclic lateral boundaries. Damping is only needed for
122!--             velocity components parallel to the outflow boundary in
123!--             the direction normal to the outflow boundary.
124                IF ( bc_ns /= 'cyclic' )  THEN
125                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
126                   kmym_y = MAX( kmym_y, km_damp_y(j) )
127                ENDIF
128
129                tend(k,j,i) = tend(k,j,i)                                    &
130                      & + 2.0 * (                                            &
131                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
132                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
133                      &         ) * ddx2                                     &
134                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
135                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
136                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
137                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
138                      &   ) * ddy
139             ENDDO
140
141!
142!--          Wall functions at the north and south walls, respectively
143             IF ( wall_u(j,i) /= 0.0 )  THEN
144                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
145                   usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
146                   usvs   = -usvs * ABS( usvs )
147                   kmyp_x = 0.25 * &
148                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
149                   kmym_x = 0.25 * &
150                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
151                   kmyp_y = kmyp_x
152                   kmym_y = kmym_x
153!
154!--                Increase diffusion at the outflow boundary in case of
155!--                non-cyclic lateral boundaries. Damping is only needed for
156!--                velocity components parallel to the outflow boundary in
157!--                the direction normal to the outflow boundary.
158                   IF ( bc_ns /= 'cyclic' )  THEN
159                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
160                      kmym_y = MAX( kmym_y, km_damp_y(j) )
161                   ENDIF
162
163                   tend(k,j,i) = tend(k,j,i)                                   &
164                                 + 2.0 * (                                     &
165                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
166                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
167                                         ) * ddx2                              &
168                                 + (   fyp(j,i) * (                            &
169                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
170                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
171                                                  )                            &
172                                     - fym(j,i) * (                            &
173                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
174                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
175                                                  )                            &
176                                     + wall_u(j,i) * usvs                      &
177                                   ) * ddy
178                ENDDO
179             ENDIF
180
181!
182!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
183!--          index k starts at nzb_u_inner+2.
184             DO  k = nzb_diff_u(j,i), nzt
185!
186!--             Interpolate eddy diffusivities on staggered gridpoints
187                kmzp = 0.25 * &
188                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
189                kmzm = 0.25 * &
190                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
191
192                tend(k,j,i) = tend(k,j,i)                                    &
193                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
194                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
195                      &            )                                         &
196                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
197                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
198                      &            )                                         &
199                      &   ) * ddzw(k)
200             ENDDO
201
202!
203!--          Vertical diffusion at the first grid point above the surface,
204!--          if the momentum flux at the bottom is given by the Prandtl law or
205!--          if it is prescribed by the user.
206!--          Difference quotient of the momentum flux is not formed over half
207!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
208!--          with other (LES) modell showed that the values of the momentum
209!--          flux becomes too large in this case.
210!--          The term containing w(k-1,..) (see above equation) is removed here
211!--          because the vertical velocity is assumed to be zero at the surface.
212             IF ( use_surface_fluxes )  THEN
213                k = nzb_u_inner(j,i)+1
214!
215!--             Interpolate eddy diffusivities on staggered gridpoints
216                kmzp = 0.25 * &
217                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
218                kmzm = 0.25 * &
219                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
220
221                tend(k,j,i) = tend(k,j,i)                                    &
222                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
223                      &   ) * ddzw(k)                                        &
224                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
225                      &   + usws(j,i)                                        &
226                      &   ) * ddzw(k)
227             ENDIF
228
229          ENDDO
230       ENDDO
231
232    END SUBROUTINE diffusion_u
233
234
235!------------------------------------------------------------------------------!
236! Call for grid point i,j
237!------------------------------------------------------------------------------!
238    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
239                               v, w, z0 )
240
241       USE control_parameters
242       USE grid_variables
243       USE indices
244
245       IMPLICIT NONE
246
247       INTEGER ::  i, j, k
248       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
249       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_y(nys-1:nyn+1)
250       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
251       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
252       REAL, DIMENSION(:,:),   POINTER ::  usws
253       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
254
255!
256!--    Compute horizontal diffusion
257       DO  k = nzb_u_outer(j,i)+1, nzt
258!
259!--       Interpolate eddy diffusivities on staggered gridpoints
260          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
261          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
262          kmyp_y = kmyp_x
263          kmym_y = kmym_x
264
265!
266!--       Increase diffusion at the outflow boundary in case of non-cyclic
267!--       lateral boundaries. Damping is only needed for velocity components
268!--       parallel to the outflow boundary in the direction normal to the
269!--       outflow boundary.
270          IF ( bc_ns /= 'cyclic' )  THEN
271             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
272             kmym_y = MAX( kmym_y, km_damp_y(j) )
273          ENDIF
274
275          tend(k,j,i) = tend(k,j,i)                                          &
276                      & + 2.0 * (                                            &
277                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
278                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
279                      &         ) * ddx2                                     &
280                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
281                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
282                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
283                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
284                      &   ) * ddy
285       ENDDO
286
287!
288!--    Wall functions at the north and south walls, respectively
289       IF ( wall_u(j,i) .NE. 0.0 )  THEN
290          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
291             usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
292             usvs   = -usvs * ABS( usvs )
293             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
294             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
295             kmyp_y = kmyp_x
296             kmym_y = kmym_x
297!
298!--          Increase diffusion at the outflow boundary in case of
299!--          non-cyclic lateral boundaries. Damping is only needed for
300!--          velocity components parallel to the outflow boundary in
301!--          the direction normal to the outflow boundary.
302             IF ( bc_ns /= 'cyclic' )  THEN
303                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
304                kmym_y = MAX( kmym_y, km_damp_y(j) )
305             ENDIF
306
307             tend(k,j,i) = tend(k,j,i)                                         &
308                                 + 2.0 * (                                     &
309                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
310                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
311                                         ) * ddx2                              &
312                                 + (   fyp(j,i) * (                            &
313                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
314                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
315                                                  )                            &
316                                     - fym(j,i) * (                            &
317                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
318                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
319                                                  )                            &
320                                     + wall_u(j,i) * usvs                      &
321                                   ) * ddy
322          ENDDO
323       ENDIF
324
325!
326!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
327!--    index k starts at nzb_u_inner+2.
328       DO  k = nzb_diff_u(j,i), nzt
329!
330!--       Interpolate eddy diffusivities on staggered gridpoints
331          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
332          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
333
334          tend(k,j,i) = tend(k,j,i)                                          &
335                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
336                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
337                      &            )                                         &
338                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
339                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
340                      &            )                                         &
341                      &   ) * ddzw(k)
342       ENDDO
343
344!
345!--    Vertical diffusion at the first grid point above the surface, if the
346!--    momentum flux at the bottom is given by the Prandtl law or if it is
347!--    prescribed by the user.
348!--    Difference quotient of the momentum flux is not formed over half of
349!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
350!--    other (LES) modell showed that the values of the momentum flux becomes
351!--    too large in this case.
352!--    The term containing w(k-1,..) (see above equation) is removed here
353!--    because the vertical velocity is assumed to be zero at the surface.
354       IF ( use_surface_fluxes )  THEN
355          k = nzb_u_inner(j,i)+1
356!
357!--       Interpolate eddy diffusivities on staggered gridpoints
358          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
359          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
360
361          tend(k,j,i) = tend(k,j,i)                                          &
362                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
363                      &   ) * ddzw(k)                                        &
364                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
365                      &   + usws(j,i)                                        &
366                      &   ) * ddzw(k)
367       ENDIF
368
369    END SUBROUTINE diffusion_u_ij
370
371 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.