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

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

last commit documented

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