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

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

leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code,
reading of dt_fixed from restart file removed

  • Property svn:keywords set to Id
File size: 15.3 KB
Line 
1 MODULE diffusion_v_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! arrays comunicated by module instead of parameter list
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_v.f90 1001 2012-09-13 14:08:46Z raasch $
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
69
70       USE arrays_3d
71       USE control_parameters
72       USE grid_variables
73       USE indices
74
75       IMPLICIT NONE
76
77       INTEGER ::  i, j, k
78       REAL    ::  kmxm, kmxp, kmzm, kmzp
79
80       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
81
82!
83!--    First calculate horizontal momentum flux v'u' at vertical walls,
84!--    if neccessary
85       IF ( topography /= 'flat' )  THEN
86          CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
87                            nzb_v_outer, wall_v )
88       ENDIF
89
90       DO  i = nxl, nxr
91          DO  j = nysv, nyn
92!
93!--          Compute horizontal diffusion
94             DO  k = nzb_v_outer(j,i)+1, nzt
95!
96!--             Interpolate eddy diffusivities on staggered gridpoints
97                kmxp = 0.25 * &
98                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
99                kmxm = 0.25 * &
100                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
101
102                tend(k,j,i) = tend(k,j,i)                                    &
103                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
104                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
105                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
106                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
107                      &   ) * ddx                                            &
108                      & + 2.0 * (                                            &
109                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
110                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
111                      &         ) * ddy2
112             ENDDO
113
114!
115!--          Wall functions at the left and right walls, respectively
116             IF ( wall_v(j,i) /= 0.0 )  THEN
117
118                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
119                   kmxp = 0.25 * &
120                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
121                   kmxm = 0.25 * &
122                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
123                   
124                   tend(k,j,i) = tend(k,j,i)                                   &
125                                 + 2.0 * (                                     &
126                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
127                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
128                                         ) * ddy2                              &
129                                 + (   fxp(j,i) * (                            &
130                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
131                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
132                                                  )                            &
133                                     - fxm(j,i) * (                            &
134                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
135                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
136                                                  )                            &
137                                     + wall_v(j,i) * vsus(k,j,i)               &
138                                   ) * ddx
139                ENDDO
140             ENDIF
141
142!
143!--          Compute vertical diffusion. In case of simulating a Prandtl
144!--          layer, index k starts at nzb_v_inner+2.
145             DO  k = nzb_diff_v(j,i), nzt_diff
146!
147!--             Interpolate eddy diffusivities on staggered gridpoints
148                kmzp = 0.25 * &
149                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
150                kmzm = 0.25 * &
151                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
152
153                tend(k,j,i) = tend(k,j,i)                                    &
154                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
155                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
156                      &            )                                         &
157                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
158                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
159                      &            )                                         &
160                      &   ) * ddzw(k)
161             ENDDO
162
163!
164!--          Vertical diffusion at the first grid point above the surface,
165!--          if the momentum flux at the bottom is given by the Prandtl law
166!--          or if it is prescribed by the user.
167!--          Difference quotient of the momentum flux is not formed over
168!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
169!--          comparison with other (LES) modell showed that the values of
170!--          the momentum flux becomes too large in this case.
171!--          The term containing w(k-1,..) (see above equation) is removed here
172!--          because the vertical velocity is assumed to be zero at the surface.
173             IF ( use_surface_fluxes )  THEN
174                k = nzb_v_inner(j,i)+1
175!
176!--             Interpolate eddy diffusivities on staggered gridpoints
177                kmzp = 0.25 * &
178                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
179                kmzm = 0.25 * &
180                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
181
182                tend(k,j,i) = tend(k,j,i)                                    &
183                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
184                      &   ) * ddzw(k)                                        &
185                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
186                      &   + vsws(j,i)                                        &
187                      &   ) * ddzw(k)
188             ENDIF
189
190!
191!--          Vertical diffusion at the first gridpoint below the top boundary,
192!--          if the momentum flux at the top is prescribed by the user
193             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
194                k = nzt
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                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
204                      &   ) * ddzw(k)                                        &
205                      & + ( -vswst(j,i)                                      &
206                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
207                      &   ) * ddzw(k)
208             ENDIF
209
210          ENDDO
211       ENDDO
212
213    END SUBROUTINE diffusion_v
214
215
216!------------------------------------------------------------------------------!
217! Call for grid point i,j
218!------------------------------------------------------------------------------!
219    SUBROUTINE diffusion_v_ij( i, j )
220
221       USE arrays_3d
222       USE control_parameters
223       USE grid_variables
224       USE indices
225
226       IMPLICIT NONE
227
228       INTEGER ::  i, j, k
229       REAL    ::  kmxm, kmxp, kmzm, kmzp
230
231       REAL, DIMENSION(nzb:nzt+1) ::  vsus
232
233!
234!--    Compute horizontal diffusion
235       DO  k = nzb_v_outer(j,i)+1, nzt
236!
237!--       Interpolate eddy diffusivities on staggered gridpoints
238          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
239          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
240
241          tend(k,j,i) = tend(k,j,i)                                          &
242                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
243                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
244                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
245                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
246                      &   ) * ddx                                            &
247                      & + 2.0 * (                                            &
248                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )    &
249                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )    &
250                      &         ) * ddy2
251       ENDDO
252
253!
254!--    Wall functions at the left and right walls, respectively
255       IF ( wall_v(j,i) /= 0.0 )  THEN
256
257!
258!--       Calculate the horizontal momentum flux v'u'
259          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
260                            vsus, 0.0, 1.0, 0.0, 0.0 )
261
262          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
263             kmxp = 0.25 * &
264                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
265             kmxm = 0.25 * &
266                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
267
268             tend(k,j,i) = tend(k,j,i)                                         &
269                                 + 2.0 * (                                     &
270                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
271                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
272                                         ) * ddy2                              &
273                                 + (   fxp(j,i) * (                            &
274                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
275                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
276                                                  )                            &
277                                     - fxm(j,i) * (                            &
278                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
279                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
280                                                  )                            &
281                                     + wall_v(j,i) * vsus(k)                   &
282                                   ) * ddx
283          ENDDO
284       ENDIF
285
286!
287!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
288!--    index k starts at nzb_v_inner+2.
289       DO  k = nzb_diff_v(j,i), nzt_diff
290!
291!--       Interpolate eddy diffusivities on staggered gridpoints
292          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
293          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
294
295          tend(k,j,i) = tend(k,j,i)                                          &
296                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)   &
297                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy         &
298                      &            )                                         &
299                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k) &
300                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy     &
301                      &            )                                         &
302                      &   ) * ddzw(k)
303       ENDDO
304
305!
306!--    Vertical diffusion at the first grid point above the surface, if the
307!--    momentum flux at the bottom is given by the Prandtl law or if it is
308!--    prescribed by the user.
309!--    Difference quotient of the momentum flux is not formed over half of
310!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
311!--    other (LES) modell showed that the values of the momentum flux becomes
312!--    too large in this case.
313!--    The term containing w(k-1,..) (see above equation) is removed here
314!--    because the vertical velocity is assumed to be zero at the surface.
315       IF ( use_surface_fluxes )  THEN
316          k = nzb_v_inner(j,i)+1
317!
318!--       Interpolate eddy diffusivities on staggered gridpoints
319          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
320          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
321
322          tend(k,j,i) = tend(k,j,i)                                          &
323                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
324                      &   ) * ddzw(k)                                        &
325                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
326                      &   + vsws(j,i)                                        &
327                      &   ) * ddzw(k)
328       ENDIF
329
330!
331!--    Vertical diffusion at the first gridpoint below the top boundary,
332!--    if the momentum flux at the top is prescribed by the user
333       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
334          k = nzt
335!
336!--       Interpolate eddy diffusivities on staggered gridpoints
337          kmzp = 0.25 * &
338                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
339          kmzm = 0.25 * &
340                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
341
342          tend(k,j,i) = tend(k,j,i)                                          &
343                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
344                      &   ) * ddzw(k)                                        &
345                      & + ( -vswst(j,i)                                      &
346                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
347                      &   ) * ddzw(k)
348       ENDIF
349
350    END SUBROUTINE diffusion_v_ij
351
352 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.