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

Last change on this file since 791 was 668, checked in by suehring, 14 years ago

last commit documented

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