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

Last change on this file since 978 was 978, checked in by fricke, 9 years ago

merge fricke branch back into trunk

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