source: palm/trunk/SOURCE/time_integration.f90 @ 77

Last change on this file since 77 was 77, checked in by raasch, 14 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: 16.3 KB
Line 
1 SUBROUTINE time_integration
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: time_integration.f90 77 2007-03-29 04:26:56Z raasch $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! Move call of user_actions( 'after_integration' ) below increment of times
14! and counters,
15! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
16! .._vector, these calls are now controlled by switch loop_optimization,
17! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz,
18! moisture renamed humidity
19!
20! RCS Log replace by Id keyword, revision history cleaned up
21!
22! Revision 1.8  2006/08/22 14:16:05  raasch
23! Disturbances are imposed only for the last Runge-Kutta-substep
24!
25! Revision 1.2  2004/04/30 13:03:40  raasch
26! decalpha-specific warning removed, routine name changed to time_integration,
27! particle advection is carried out only once during the intermediate steps,
28! impulse_advec renamed momentum_advec
29!
30! Revision 1.1  1997/08/11 06:19:04  raasch
31! Initial revision
32!
33!
34! Description:
35! ------------
36! Integration in time of the model equations, statistical analysis and graphic
37! output
38!------------------------------------------------------------------------------!
39
40    USE arrays_3d
41    USE averaging
42    USE control_parameters
43    USE cpulog
44#if defined( __dvrp_graphics )
45    USE DVRP
46#endif
47    USE grid_variables
48    USE indices
49    USE interaction_droplets_ptq_mod
50    USE interfaces
51    USE particle_attributes
52    USE pegrid
53    USE prognostic_equations_mod
54    USE statistics
55    USE user_actions_mod
56
57    IMPLICIT NONE
58
59    CHARACTER (LEN=9) ::  time_to_string
60    INTEGER ::  i, j, k
61
62!
63!-- At the beginning of a simulation determine the time step as well as
64!-- determine and print out the run control parameters
65    IF ( simulated_time == 0.0 )  CALL timestep
66    CALL run_control
67
68#if defined( __dvrp_graphics )
69!
70!-- Time measurement with dvrp software 
71    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
72#endif
73
74!
75!-- Start of the time loop
76    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
77                .NOT. terminate_run )
78
79       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
80
81!
82!--    Determine size of next time step
83       IF ( simulated_time /= 0.0 )  CALL timestep
84
85!
86!--    Execute the user-defined actions
87       CALL user_actions( 'before_timestep' )
88
89!
90!--    Start of intermediate step loop
91       intermediate_timestep_count = 0
92       DO  WHILE ( intermediate_timestep_count < &
93                   intermediate_timestep_count_max )
94
95          intermediate_timestep_count = intermediate_timestep_count + 1
96
97!
98!--       Set the steering factors for the prognostic equations which depend
99!--       on the timestep scheme
100          CALL timestep_scheme_steering
101
102!
103!--       Solve the prognostic equations. A fast cache optimized version with
104!--       only one single loop is used in case of Piascek-Williams advection
105!--       scheme. NEC vector machines use a different version, because
106!--       in the other versions a good vectorization is prohibited due to
107!--       inlining problems.
108          IF ( loop_optimization == 'vector' )  THEN
109             CALL prognostic_equations_vector
110          ELSE
111             IF ( momentum_advec == 'ups-scheme'  .OR.  &
112                  scalar_advec == 'ups-scheme'   .OR.  &
113                  scalar_advec == 'bc-scheme' )        &
114             THEN
115                CALL prognostic_equations_noopt
116             ELSE
117                CALL prognostic_equations_cache
118             ENDIF
119          ENDIF
120
121!
122!--       Particle advection (only once during intermediate steps, because
123!--       it uses an Euler-step)
124          IF ( particle_advection  .AND.                         &
125               simulated_time >= particle_advection_start  .AND. &
126               intermediate_timestep_count == 1 )  THEN
127             CALL advec_particles
128             first_call_advec_particles = .FALSE.
129          ENDIF
130
131!
132!--       Interaction of droplets with temperature and specific humidity.
133!--       Droplet condensation and evaporation is calculated within
134!--       advec_particles.
135          IF ( cloud_droplets  .AND.  &
136               intermediate_timestep_count == intermediate_timestep_count_max )&
137          THEN
138             CALL interaction_droplets_ptq
139          ENDIF
140
141!
142!--       Exchange of ghost points (lateral boundary conditions)
143          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
144          CALL exchange_horiz( u_p )
145          CALL exchange_horiz( v_p )
146          CALL exchange_horiz( w_p )
147          CALL exchange_horiz( pt_p )
148          IF ( .NOT. constant_diffusion       )  CALL exchange_horiz( e_p )
149          IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz( q_p )
150          IF ( cloud_droplets )  THEN
151             CALL exchange_horiz( ql )
152             CALL exchange_horiz( ql_c )
153             CALL exchange_horiz( ql_v )
154             CALL exchange_horiz( ql_vp )
155          ENDIF
156
157          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
158
159!
160!--       Apply time filter in case of leap-frog timestep
161          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
162             CALL asselin_filter
163          ENDIF
164
165!
166!--       Boundary conditions for the prognostic quantities (except of the
167!--       velocities at the outflow in case of a non-cyclic lateral wall)
168          CALL boundary_conds( 'main' )
169
170!
171!--       Swap the time levels in preparation for the next time step.
172          CALL swap_timelevel
173
174!
175!--       Temperature offset must be imposed at cyclic boundaries in x-direction
176!--       when a sloping surface is used
177          IF ( sloping_surface )  THEN
178             IF ( nxl ==  0 )  pt(:,:,nxl-1) = pt(:,:,nxl-1) - pt_slope_offset
179             IF ( nxr == nx )  pt(:,:,nxr+1) = pt(:,:,nxr+1) + pt_slope_offset
180          ENDIF
181
182!
183!--       Impose a random perturbation on the horizontal velocity field
184          IF ( create_disturbances  .AND.  &
185               intermediate_timestep_count == intermediate_timestep_count_max )&
186          THEN
187             time_disturb = time_disturb + dt_3d
188             IF ( time_disturb >= dt_disturb )  THEN
189                IF ( hom(nzb+5,1,var_hom,0) < disturbance_energy_limit )  THEN
190                   CALL disturb_field( nzb_u_inner, tend, u )
191                   CALL disturb_field( nzb_v_inner, tend, v )
192                ELSEIF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
193!
194!--                Runs with a non-cyclic lateral wall need perturbations
195!--                near the inflow throughout the whole simulation
196                   dist_range = 1
197                   CALL disturb_field( nzb_u_inner, tend, u )
198                   CALL disturb_field( nzb_v_inner, tend, v )
199                   dist_range = 0
200                ENDIF
201                time_disturb = time_disturb - dt_disturb
202             ENDIF
203          ENDIF
204
205!
206!--       Reduce the velocity divergence via the equation for perturbation
207!--       pressure.
208          IF ( intermediate_timestep_count == intermediate_timestep_count_max &
209               .OR.  call_psolver_at_all_substeps )  THEN
210             CALL pres
211          ENDIF
212
213!
214!--       In case of a non-cyclic lateral wall, set the boundary conditions for
215!--       the velocities at the outflow
216!          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
217!             CALL boundary_conds( 'outflow_uvw' )
218!          ENDIF
219
220!
221!--       If required, compute virtuell potential temperature
222          IF ( humidity ) CALL compute_vpt
223
224!
225!--       If required, compute liquid water content
226          IF ( cloud_physics ) CALL calc_liquid_water_content
227
228!
229!--       Compute the diffusion quantities
230          IF ( .NOT. constant_diffusion )  THEN
231
232!
233!--          First the vertical fluxes in the Prandtl layer are being computed
234             IF ( prandtl_layer )  THEN
235                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
236                CALL prandtl_fluxes
237                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
238             ENDIF
239
240!
241!--          Compute the diffusion coefficients
242             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
243             IF ( .NOT. humidity ) THEN
244                CALL diffusivities( pt )
245             ELSE
246                CALL diffusivities( vpt )
247             ENDIF
248             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
249
250          ENDIF
251
252       ENDDO   ! Intermediate step loop
253
254!
255!--    Increase simulation time and output times
256       current_timestep_number = current_timestep_number + 1
257       simulated_time     = simulated_time   + dt_3d
258       simulated_time_chr = time_to_string( simulated_time )
259       IF ( simulated_time >= skip_time_data_output_av )  THEN
260          time_do_av         = time_do_av       + dt_3d
261       ENDIF
262       IF ( simulated_time >= skip_time_do2d_xy )  THEN
263          time_do2d_xy       = time_do2d_xy     + dt_3d
264       ENDIF
265       IF ( simulated_time >= skip_time_do2d_xz )  THEN
266          time_do2d_xz       = time_do2d_xz     + dt_3d
267       ENDIF
268       IF ( simulated_time >= skip_time_do2d_yz )  THEN
269          time_do2d_yz       = time_do2d_yz     + dt_3d
270       ENDIF
271       IF ( simulated_time >= skip_time_do3d    )  THEN
272          time_do3d          = time_do3d        + dt_3d
273       ENDIF
274       time_dvrp          = time_dvrp        + dt_3d
275       IF ( simulated_time >= skip_time_dosp )  THEN
276          time_dosp       = time_dosp        + dt_3d
277       ENDIF
278       time_dots          = time_dots        + dt_3d
279       IF ( .NOT. first_call_advec_particles )  THEN
280          time_dopts      = time_dopts       + dt_3d
281       ENDIF
282       IF ( simulated_time >= skip_time_dopr )  THEN
283          time_dopr       = time_dopr        + dt_3d
284       ENDIF
285       time_dopr_listing          = time_dopr_listing        + dt_3d
286       time_run_control   = time_run_control + dt_3d
287
288!
289!--    Execute user-defined actions
290       CALL user_actions( 'after_integration' )
291
292!
293!--    If Galilei transformation is used, determine the distance that the
294!--    model has moved so far
295       IF ( galilei_transformation )  THEN
296          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
297          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
298       ENDIF
299
300!
301!--    Check, if restart is necessary (because cpu-time is expiring or
302!--    because it is forced by user) and set stop flag
303       CALL check_for_restart
304
305!
306!--    Carry out statistical analysis and output at the requested output times.
307!--    The MOD function is used for calculating the output time counters (like
308!--    time_dopr) in order to regard a possible decrease of the output time
309!--    interval in case of restart runs
310
311!
312!--    Set a flag indicating that so far no statistics have been created
313!--    for this time step
314       flow_statistics_called = .FALSE.
315
316!
317!--    If required, call flow_statistics for averaging in time
318       IF ( averaging_interval_pr /= 0.0  .AND.  &
319            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
320            simulated_time >= skip_time_dopr )  THEN
321          time_dopr_av = time_dopr_av + dt_3d
322          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
323             do_sum = .TRUE.
324             time_dopr_av = MOD( time_dopr_av, &
325                                    MAX( dt_averaging_input_pr, dt_3d ) )
326          ENDIF
327       ENDIF
328       IF ( do_sum )  CALL flow_statistics
329
330!
331!--    Sum-up 3d-arrays for later output of time-averaged data
332       IF ( averaging_interval /= 0.0  .AND.                                &
333            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
334            simulated_time >= skip_time_data_output_av )                    &
335       THEN
336          time_do_sla = time_do_sla + dt_3d
337          IF ( time_do_sla >= dt_averaging_input )  THEN
338             CALL sum_up_3d_data
339             average_count_3d = average_count_3d + 1
340             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
341          ENDIF
342       ENDIF
343
344!
345!--    Calculate spectra for time averaging
346       IF ( averaging_interval_sp /= 0.0  .AND.  &
347            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
348            simulated_time >= skip_time_dosp )  THEN
349          time_dosp_av = time_dosp_av + dt_3d
350          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
351             CALL calc_spectra
352             time_dosp_av = MOD( time_dosp_av, &
353                                 MAX( dt_averaging_input_pr, dt_3d ) )
354          ENDIF
355       ENDIF
356
357!
358!--    Computation and output of run control parameters.
359!--    This is also done whenever the time step has changed or perturbations
360!--    have been imposed
361       IF ( time_run_control >= dt_run_control  .OR.                     &
362            ( dt_changed  .AND.  timestep_scheme(1:5) /= 'runge' )  .OR. &
363            disturbance_created )                                        &
364       THEN
365          CALL run_control
366          IF ( time_run_control >= dt_run_control )  THEN
367             time_run_control = MOD( time_run_control, &
368                                     MAX( dt_run_control, dt_3d ) )
369          ENDIF
370       ENDIF
371
372!
373!--    Profile output (ASCII) on file
374       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
375          CALL print_1d
376          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
377                                                           dt_3d ) )
378       ENDIF
379
380!
381!--    Graphic output for PROFIL
382       IF ( time_dopr >= dt_dopr )  THEN
383          IF ( dopr_n /= 0 )  CALL data_output_profiles
384          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
385          time_dopr_av = 0.0    ! due to averaging (see above)
386       ENDIF
387
388!
389!--    Graphic output for time series
390       IF ( time_dots >= dt_dots )  THEN
391          CALL data_output_tseries
392          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
393       ENDIF
394
395!
396!--    Output of spectra (formatted for use with PROFIL), in case of no
397!--    time averaging, spectra has to be calculated before
398       IF ( time_dosp >= dt_dosp )  THEN
399          IF ( average_count_sp == 0 )  CALL calc_spectra
400          CALL data_output_spectra
401          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
402       ENDIF
403
404!
405!--    2d-data output (cross-sections)
406       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
407          CALL data_output_2d( 'xy', 0 )
408          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
409       ENDIF
410       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
411          CALL data_output_2d( 'xz', 0 )
412          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
413       ENDIF
414       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
415          CALL data_output_2d( 'yz', 0 )
416          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
417       ENDIF
418
419!
420!--    3d-data output (volume data)
421       IF ( time_do3d >= dt_do3d )  THEN
422          CALL data_output_3d( 0 )
423          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
424       ENDIF
425
426!
427!--    Output of time-averaged 2d/3d-data
428       IF ( time_do_av >= dt_data_output_av )  THEN
429          CALL average_3d_data
430          CALL data_output_2d( 'xy', 1 )
431          CALL data_output_2d( 'xz', 1 )
432          CALL data_output_2d( 'yz', 1 )
433          CALL data_output_3d( 1 )
434          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
435       ENDIF
436
437!
438!--    Output of particle time series
439       IF ( time_dopts >= dt_dopts  .OR. &
440            ( simulated_time >= particle_advection_start  .AND. &
441              first_call_advec_particles ) )  THEN
442          CALL data_output_ptseries
443          time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
444       ENDIF
445
446!
447!--    Output of dvrp-graphics (isosurface, particles, slicer)
448#if defined( __dvrp_graphics )
449       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
450#endif
451       IF ( time_dvrp >= dt_dvrp )  THEN
452          CALL data_output_dvrp
453          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
454       ENDIF
455#if defined( __dvrp_graphics )
456       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
457#endif
458
459!
460!--    If required, set the heat flux for the next time step at a random value
461       IF ( constant_heatflux  .AND.  random_heatflux )  CALL disturb_heatflux
462
463!
464!--    Execute user-defined actions
465       CALL user_actions( 'after_timestep' )
466
467       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
468
469    ENDDO   ! time loop
470
471#if defined( __dvrp_graphics )
472    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
473#endif
474
475 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.