source: palm/trunk/SOURCE/diffusion_v.f90 @ 987

Last change on this file since 987 was 979, checked in by fricke, 12 years ago

last commit documented

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