Home


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

source: palm/trunk/SOURCE/diffusion_u.f90 @ 250

Last change on this file since 250 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.7 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_u.f90 110 2007-10-05 05:13:14Z letzel $
11!
12! 106 2007-08-16 14:30:26Z raasch
13! Momentumflux at top (uswst) included as boundary condition,
14! i loop is starting from nxlu (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, uxrp 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:35:35  raasch
26! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
27! or nzb_diff_u, 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:23:51  raasch
34! Initial revision
35!
36!
37! Description:
38! ------------
39! Diffusion term of the u-component
40! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
41!        and slows down the speed on NEC about 5-10%
42!------------------------------------------------------------------------------!
43
44    USE wall_fluxes_mod
45
46    PRIVATE
47    PUBLIC diffusion_u
48
49    INTERFACE diffusion_u
50       MODULE PROCEDURE diffusion_u
51       MODULE PROCEDURE diffusion_u_ij
52    END INTERFACE diffusion_u
53
54 CONTAINS
55
56
57!------------------------------------------------------------------------------!
58! Call for all grid points
59!------------------------------------------------------------------------------!
60    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
61                            v, w )
62
63       USE control_parameters
64       USE grid_variables
65       USE indices
66
67       IMPLICIT NONE
68
69       INTEGER ::  i, j, k
70       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
71       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
72       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
73       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
74       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
75       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
76
77!
78!--    First calculate horizontal momentum flux u'v' at vertical walls,
79!--    if neccessary
80       IF ( topography /= 'flat' )  THEN
81          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
82                            nzb_u_outer, wall_u )
83       ENDIF
84
85       DO  i = nxlu, nxr
86          DO  j = nys,nyn
87!
88!--          Compute horizontal diffusion
89             DO  k = nzb_u_outer(j,i)+1, nzt
90!
91!--             Interpolate eddy diffusivities on staggered gridpoints
92                kmyp_x = 0.25 * &
93                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
94                kmym_x = 0.25 * &
95                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
96                kmyp_y = kmyp_x
97                kmym_y = kmym_x
98!
99!--             Increase diffusion at the outflow boundary in case of
100!--             non-cyclic lateral boundaries. Damping is only needed for
101!--             velocity components parallel to the outflow boundary in
102!--             the direction normal to the outflow boundary.
103                IF ( bc_ns /= 'cyclic' )  THEN
104                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
105                   kmym_y = MAX( kmym_y, km_damp_y(j) )
106                ENDIF
107
108                tend(k,j,i) = tend(k,j,i)                                    &
109                      & + 2.0 * (                                            &
110                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
111                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
112                      &         ) * ddx2                                     &
113                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
114                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
115                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
116                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
117                      &   ) * ddy
118             ENDDO
119
120!
121!--          Wall functions at the north and south walls, respectively
122             IF ( wall_u(j,i) /= 0.0 )  THEN
123
124                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
125                   kmyp_x = 0.25 * &
126                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
127                   kmym_x = 0.25 * &
128                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
129                   kmyp_y = kmyp_x
130                   kmym_y = kmym_x
131!
132!--                Increase diffusion at the outflow boundary in case of
133!--                non-cyclic lateral boundaries. Damping is only needed for
134!--                velocity components parallel to the outflow boundary in
135!--                the direction normal to the outflow boundary.
136                   IF ( bc_ns /= 'cyclic' )  THEN
137                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
138                      kmym_y = MAX( kmym_y, km_damp_y(j) )
139                   ENDIF
140
141                   tend(k,j,i) = tend(k,j,i)                                   &
142                                 + 2.0 * (                                     &
143                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
144                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
145                                         ) * ddx2                              &
146                                 + (   fyp(j,i) * (                            &
147                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
148                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
149                                                  )                            &
150                                     - fym(j,i) * (                            &
151                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
152                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
153                                                  )                            &
154                                     + wall_u(j,i) * usvs(k,j,i)               &
155                                   ) * ddy
156                ENDDO
157             ENDIF
158
159!
160!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
161!--          index k starts at nzb_u_inner+2.
162             DO  k = nzb_diff_u(j,i), nzt_diff
163!
164!--             Interpolate eddy diffusivities on staggered gridpoints
165                kmzp = 0.25 * &
166                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
167                kmzm = 0.25 * &
168                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
169
170                tend(k,j,i) = tend(k,j,i)                                    &
171                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
172                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
173                      &            )                                         &
174                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
175                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
176                      &            )                                         &
177                      &   ) * ddzw(k)
178             ENDDO
179
180!
181!--          Vertical diffusion at the first grid point above the surface,
182!--          if the momentum flux at the bottom is given by the Prandtl law or
183!--          if it is prescribed by the user.
184!--          Difference quotient of the momentum flux is not formed over half
185!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
186!--          with other (LES) modell showed that the values of the momentum
187!--          flux becomes too large in this case.
188!--          The term containing w(k-1,..) (see above equation) is removed here
189!--          because the vertical velocity is assumed to be zero at the surface.
190             IF ( use_surface_fluxes )  THEN
191                k = nzb_u_inner(j,i)+1
192!
193!--             Interpolate eddy diffusivities on staggered gridpoints
194                kmzp = 0.25 * &
195                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
196                kmzm = 0.25 * &
197                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
198
199                tend(k,j,i) = tend(k,j,i)                                    &
200                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
201                      &   ) * ddzw(k)                                        &
202                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
203                      &   + usws(j,i)                                        &
204                      &   ) * ddzw(k)
205             ENDIF
206
207!
208!--          Vertical diffusion at the first gridpoint below the top boundary,
209!--          if the momentum flux at the top is prescribed by the user
210             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
211                k = nzt
212!
213!--             Interpolate eddy diffusivities on staggered gridpoints
214                kmzp = 0.25 * &
215                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
216                kmzm = 0.25 * &
217                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
218
219                tend(k,j,i) = tend(k,j,i)                                    &
220                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
221                      &   ) * ddzw(k)                                        &
222                      & + ( -uswst(j,i)                                      &
223                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
224                      &   ) * ddzw(k)
225             ENDIF
226
227          ENDDO
228       ENDDO
229
230    END SUBROUTINE diffusion_u
231
232
233!------------------------------------------------------------------------------!
234! Call for grid point i,j
235!------------------------------------------------------------------------------!
236    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
237                               uswst, v, w )
238
239       USE control_parameters
240       USE grid_variables
241       USE indices
242
243       IMPLICIT NONE
244
245       INTEGER ::  i, j, k
246       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
247       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
248       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
249       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
250       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
251       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
252
253!
254!--    Compute horizontal diffusion
255       DO  k = nzb_u_outer(j,i)+1, nzt
256!
257!--       Interpolate eddy diffusivities on staggered gridpoints
258          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
259          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
260          kmyp_y = kmyp_x
261          kmym_y = kmym_x
262
263!
264!--       Increase diffusion at the outflow boundary in case of non-cyclic
265!--       lateral boundaries. Damping is only needed for velocity components
266!--       parallel to the outflow boundary in the direction normal to the
267!--       outflow boundary.
268          IF ( bc_ns /= 'cyclic' )  THEN
269             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
270             kmym_y = MAX( kmym_y, km_damp_y(j) )
271          ENDIF
272
273          tend(k,j,i) = tend(k,j,i)                                          &
274                      & + 2.0 * (                                            &
275                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
276                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
277                      &         ) * ddx2                                     &
278                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
279                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
280                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
281                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
282                      &   ) * ddy
283       ENDDO
284
285!
286!--    Wall functions at the north and south walls, respectively
287       IF ( wall_u(j,i) .NE. 0.0 )  THEN
288
289!
290!--       Calculate the horizontal momentum flux u'v'
291          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
292                            usvs, 1.0, 0.0, 0.0, 0.0 )
293
294          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
295             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
296             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
297             kmyp_y = kmyp_x
298             kmym_y = kmym_x
299!
300!--          Increase diffusion at the outflow boundary in case of
301!--          non-cyclic lateral boundaries. Damping is only needed for
302!--          velocity components parallel to the outflow boundary in
303!--          the direction normal to the outflow boundary.
304             IF ( bc_ns /= 'cyclic' )  THEN
305                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
306                kmym_y = MAX( kmym_y, km_damp_y(j) )
307             ENDIF
308
309             tend(k,j,i) = tend(k,j,i)                                         &
310                                 + 2.0 * (                                     &
311                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
312                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
313                                         ) * ddx2                              &
314                                 + (   fyp(j,i) * (                            &
315                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
316                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
317                                                  )                            &
318                                     - fym(j,i) * (                            &
319                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
320                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
321                                                  )                            &
322                                     + wall_u(j,i) * usvs(k)                   &
323                                   ) * ddy
324          ENDDO
325       ENDIF
326
327!
328!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
329!--    index k starts at nzb_u_inner+2.
330       DO  k = nzb_diff_u(j,i), nzt_diff
331!
332!--       Interpolate eddy diffusivities on staggered gridpoints
333          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
334          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
335
336          tend(k,j,i) = tend(k,j,i)                                          &
337                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
338                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
339                      &            )                                         &
340                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
341                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
342                      &            )                                         &
343                      &   ) * ddzw(k)
344       ENDDO
345
346!
347!--    Vertical diffusion at the first grid point above the surface, if the
348!--    momentum flux at the bottom is given by the Prandtl law or if it is
349!--    prescribed by the user.
350!--    Difference quotient of the momentum flux is not formed over half of
351!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
352!--    other (LES) modell showed that the values of the momentum flux becomes
353!--    too large in this case.
354!--    The term containing w(k-1,..) (see above equation) is removed here
355!--    because the vertical velocity is assumed to be zero at the surface.
356       IF ( use_surface_fluxes )  THEN
357          k = nzb_u_inner(j,i)+1
358!
359!--       Interpolate eddy diffusivities on staggered gridpoints
360          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
361          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
362
363          tend(k,j,i) = tend(k,j,i)                                          &
364                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
365                      &   ) * ddzw(k)                                        &
366                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
367                      &   + usws(j,i)                                        &
368                      &   ) * ddzw(k)
369       ENDIF
370
371!
372!--    Vertical diffusion at the first gridpoint below the top boundary,
373!--    if the momentum flux at the top is prescribed by the user
374       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
375          k = nzt
376!
377!--       Interpolate eddy diffusivities on staggered gridpoints
378          kmzp = 0.25 * &
379                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
380          kmzm = 0.25 * &
381                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
382
383          tend(k,j,i) = tend(k,j,i)                                          &
384                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
385                      &   ) * ddzw(k)                                        &
386                      & + ( -uswst(j,i)                                      &
387                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
388                      &   ) * ddzw(k)
389       ENDIF
390
391    END SUBROUTINE diffusion_u_ij
392
393 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |