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

Last change on this file since 1002 was 1002, checked in by raasch, 9 years ago

last commit documented

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