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

Last change on this file since 978 was 978, checked in by fricke, 9 years ago

merge fricke branch back into trunk

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