Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 354 was 110, checked in by raasch, 17 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:

Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:

Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

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