Home


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

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

Last change on this file since 84 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

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