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

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

Initial repository layout and content

File size: 16.7 KB
Line 
1 MODULE diffusion_v_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: diffusion_v.f90,v $
11! Revision 1.15  2006/02/23 10:36:00  raasch
12! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner
13! or nzb_diff_v, 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:09:56  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_x
23!
24! Revision 1.13  2004/01/30 10:21:56  raasch
25! Scalar lower k index nzb replaced by 2d-array nzb_2d
26!
27! Revision 1.12  2003/03/12 16:26:31  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:46  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:29:31  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:19  raasch
41! Translation of remaining German identifiers (variables, subroutines, etc.)
42!
43! Revision 1.8  2001/01/22 06:18:41  raasch
44! Module test_variables removed
45!
46! Revision 1.7  2000/07/03 12:58:14  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:58:18  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:12:05  12:12:05  raasch (Siegfried Raasch)
56! + USE test_variables
57!
58! Revision 1.4  1998/01/23 09:57:18  raasch
59! Einbau der Prandtl-Schicht
60!
61! Revision 1.3  1997/09/12 07:24:20  raasch
62! Leerzeilen mussten entfernt werden
63!
64! Revision 1.2  1997/09/12 06:43:46  raasch
65! HP-Compiler erfordert & am Beginn von Fortsetzungszeilen
66!
67! Revision 1.1  1997/09/12 06:24:01  raasch
68! Initial revision
69!
70!
71! Description:
72! ------------
73! Diffusion term of the v-component
74!------------------------------------------------------------------------------!
75
76    PRIVATE
77    PUBLIC diffusion_v
78
79    INTERFACE diffusion_v
80       MODULE PROCEDURE diffusion_v
81       MODULE PROCEDURE diffusion_v_ij
82    END INTERFACE diffusion_v
83
84 CONTAINS
85
86
87!------------------------------------------------------------------------------!
88! Call for all grid points
89!------------------------------------------------------------------------------!
90    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, 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    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp, vsus
100       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+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 ::  vsws
104       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
105
106       DO  i = nxl, nxr
107          DO  j = nys, nyn+vynp
108!
109!--          Compute horizontal diffusion
110             DO  k = nzb_v_outer(j,i)+1, nzt
111!
112!--             Interpolate eddy diffusivities on staggered gridpoints
113                kmxp_x = 0.25 * &
114                         ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
115                kmxm_x = 0.25 * &
116                         ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
117                kmxp_y = kmxp_x
118                kmxm_y = kmxm_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_lr /= 'cyclic' )  THEN
125                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
126                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
127                ENDIF
128
129                tend(k,j,i) = tend(k,j,i)                                    &
130                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
131                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
132                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
133                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
134                      &   ) * ddx                                            &
135                      & + 2.0 * (                                            &
136                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
137                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
138                      &         ) * ddy2
139             ENDDO
140
141!
142!--          Wall functions at the left and right walls, respectively
143             IF ( wall_v(j,i) /= 0.0 )  THEN
144                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
145                   vsus   = kappa * v(k,j,i) / LOG( 0.5 * dx / z0(j,i))
146                   vsus   = -vsus * ABS( vsus )
147                   kmxp_x = 0.25 * &
148                            ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
149                   kmxm_x = 0.25 * &
150                            ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
151                   kmxp_y = kmxp_x
152                   kmxm_y = kmxm_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_lr /= 'cyclic' )  THEN
159                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
160                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
161                   ENDIF
162
163                   tend(k,j,i) = tend(k,j,i)                                   &
164                                 + 2.0 * (                                     &
165                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
166                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
167                                         ) * ddy2                              &
168                                 + (   fxp(j,i) * (                            &
169                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
170                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
171                                                  )                            &
172                                     - fxm(j,i) * (                            &
173                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
174                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
175                                                  )                            &
176                                     + wall_v(j,i) * vsus                      &
177                                   ) * ddx
178                ENDDO
179             ENDIF
180
181!
182!--          Compute vertical diffusion. In case of simulating a Prandtl
183!--          layer, index k starts at nzb_v_inner+2.
184             DO  k = nzb_diff_v(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-1,i)+km(k+1,j-1,i) )
189                kmzm = 0.25 * &
190                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
191
192                tend(k,j,i) = tend(k,j,i)                                    &
193                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
194                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
195                      &            )                                         &
196                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
197                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
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
205!--          or if it is prescribed by the user.
206!--          Difference quotient of the momentum flux is not formed over
207!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
208!--          comparison with other (LES) modell showed that the values of
209!--          the momentum 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_v_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-1,i)+km(k+1,j-1,i) )
218                kmzm = 0.25 * &
219                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
220
221                tend(k,j,i) = tend(k,j,i)                                    &
222                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
223                      &   ) * ddzw(k)                                        &
224                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
225                      &   + vsws(j,i)                                        &
226                      &   ) * ddzw(k)
227             ENDIF
228
229          ENDDO
230       ENDDO
231
232    END SUBROUTINE diffusion_v
233
234
235!------------------------------------------------------------------------------!
236! Call for grid point i,j
237!------------------------------------------------------------------------------!
238    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
239                               vsws, 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    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp, vsus
249       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+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 ::  vsws
253       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
254
255!
256!--    Compute horizontal diffusion
257       DO  k = nzb_v_outer(j,i)+1, nzt
258!
259!--       Interpolate eddy diffusivities on staggered gridpoints
260          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
261          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
262          kmxp_y = kmxp_x
263          kmxm_y = kmxm_x
264!
265!--       Increase diffusion at the outflow boundary in case of non-cyclic
266!--       lateral boundaries. Damping is only needed for velocity components
267!--       parallel to the outflow boundary in the direction normal to the
268!--       outflow boundary.
269          IF ( bc_lr /= 'cyclic' )  THEN
270             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
271             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
272          ENDIF
273
274          tend(k,j,i) = tend(k,j,i)                                          &
275                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
276                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
277                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
278                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
279                      &   ) * ddx                                            &
280                      & + 2.0 * (                                            &
281                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
282                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
283                      &         ) * ddy2
284       ENDDO
285
286!
287!--    Wall functions at the left and right walls, respectively
288       IF ( wall_v(j,i) /= 0.0 )  THEN
289          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
290             vsus   = kappa * v(k,j,i) / LOG( 0.5 * dx / z0(j,i))
291             vsus   = -vsus * ABS( vsus )
292             kmxp_x = 0.25 * &
293                      ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
294             kmxm_x = 0.25 * &
295                      ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
296             kmxp_y = kmxp_x
297             kmxm_y = kmxm_x
298!
299!--          Increase diffusion at the outflow boundary in case of
300!--          non-cyclic lateral boundaries. Damping is only needed for
301!--          velocity components parallel to the outflow boundary in
302!--          the direction normal to the outflow boundary.
303             IF ( bc_lr /= 'cyclic' )  THEN
304                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
305                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
306             ENDIF
307
308             tend(k,j,i) = tend(k,j,i)                                         &
309                                 + 2.0 * (                                     &
310                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
311                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
312                                         ) * ddy2                              &
313                                 + (   fxp(j,i) * (                            &
314                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
315                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
316                                                  )                            &
317                                     - fxm(j,i) * (                            &
318                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
319                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
320                                                  )                            &
321                                     + wall_v(j,i) * vsus                      &
322                                   ) * ddx
323          ENDDO
324       ENDIF
325
326!
327!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
328!--    index k starts at nzb_v_inner+2.
329       DO  k = nzb_diff_v(j,i), nzt
330!
331!--       Interpolate eddy diffusivities on staggered gridpoints
332          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
333          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
334
335          tend(k,j,i) = tend(k,j,i)                                          &
336                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
337                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
338                      &            )                                         &
339                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
340                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
341                      &            )                                         &
342                      &   ) * ddzw(k)
343       ENDDO
344
345!
346!--    Vertical diffusion at the first grid point above the surface, if the
347!--    momentum flux at the bottom is given by the Prandtl law or if it is
348!--    prescribed by the user.
349!--    Difference quotient of the momentum flux is not formed over half of
350!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
351!--    other (LES) modell showed that the values of the momentum flux becomes
352!--    too large in this case.
353!--    The term containing w(k-1,..) (see above equation) is removed here
354!--    because the vertical velocity is assumed to be zero at the surface.
355       IF ( use_surface_fluxes )  THEN
356          k = nzb_v_inner(j,i)+1
357!
358!--       Interpolate eddy diffusivities on staggered gridpoints
359          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
360          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
361
362          tend(k,j,i) = tend(k,j,i)                                          &
363                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
364                      &   ) * ddzw(k)                                        &
365                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
366                      &   + vsws(j,i)                                        &
367                      &   ) * ddzw(k)
368       ENDIF
369
370    END SUBROUTINE diffusion_v_ij
371
372 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.