source: palm/trunk/SOURCE/timestep.f90 @ 519

Last change on this file since 519 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

  • Property svn:keywords set to Id
File size: 15.2 KB
RevLine 
[1]1 SUBROUTINE timestep
2
3!------------------------------------------------------------------------------!
[258]4! Current revisions:
[1]5! -----------------
[392]6!
[316]7!
[1]8! Former revisions:
9! -----------------
[3]10! $Id: timestep.f90 392 2009-09-24 10:39:14Z raasch $
[110]11!
[392]12! 343 2009-06-24 12:59:09Z maronga
13! Additional timestep criterion in case of simulations with plant canopy
14! Output of messages replaced by message handling routine.
15!
[226]16! 222 2009-01-12 16:04:16Z letzel
17! Implementation of a MPI-1 Coupling: replaced myid with target_id
18! Bugfix for nonparallel execution
19!
[110]20! 108 2007-08-24 15:10:38Z letzel
21! modifications to terminate coupled runs
22!
[3]23! RCS Log replace by Id keyword, revision history cleaned up
24!
[1]25! Revision 1.21  2006/02/23 12:59:44  raasch
26! nt_anz renamed current_timestep_number
27!
28! Revision 1.1  1997/08/11 06:26:19  raasch
29! Initial revision
30!
31!
32! Description:
33! ------------
34! Compute the time step under consideration of the FCL and diffusion criterion.
35!------------------------------------------------------------------------------!
36
37    USE arrays_3d
38    USE control_parameters
39    USE cpulog
40    USE grid_variables
41    USE indices
42    USE interfaces
43    USE pegrid
44    USE statistics
45
46    IMPLICIT NONE
47
48    INTEGER ::  i, j, k
49
[318]50    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,                 &
51             dt_plant_canopy_l,                                        &
52             dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,  & 
53             dt_u, dt_v, dt_w, lad_max, percent_change,                &
[316]54             u_gtrans_l, vabs_max, value, v_gtrans_l
[1]55
56    REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
57    REAL, DIMENSION(nzb+1:nzt) ::  dxyz2_min
58
59    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
60
61!
62!-- Determine the maxima of the velocity components.
63    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, u, 'abs', &
64                         u_max, u_max_ijk )
65    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, v, 'abs', &
66                         v_max, v_max_ijk )
67    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, w, 'abs', &
68                         w_max, w_max_ijk )
69
70!
71!-- In case maxima of the horizontal velocity components have been found at the
72!-- bottom boundary (k=nzb), the corresponding maximum at level k=1 is chosen
73!-- if the Dirichlet-boundary condition ('mirror') has been set. This is
74!-- necessary, because otherwise in case of Galilei-transform a far too large
75!-- velocity (having the respective opposite sign) would be used for the time
76!-- step determination (almost double the mean flow velocity).
77    IF ( ibc_uv_b == 0 )  THEN
78       IF ( u_max_ijk(1) == nzb )  THEN
79          u_max        = -u_max
80          u_max_ijk(1) = nzb + 1
81       ENDIF
82       IF ( v_max_ijk(1) == nzb )  THEN
83          v_max        = -v_max
84          v_max_ijk(1) = nzb + 1
85       ENDIF
86    ENDIF
87
88!
89!-- In case of Galilei-transform not using the geostrophic wind as translation
90!-- velocity, compute the volume-averaged horizontal velocity components, which
91!-- will then be subtracted from the horizontal wind for the time step and
92!-- horizontal advection routines.
93    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
94       IF ( flow_statistics_called )  THEN
95!
96!--       Horizontal averages already existent, just need to average them
97!--       vertically.
98          u_gtrans = 0.0
99          v_gtrans = 0.0
100          DO  k = nzb+1, nzt
101             u_gtrans = u_gtrans + hom(k,1,1,0)
102             v_gtrans = v_gtrans + hom(k,1,2,0)
103          ENDDO
104          u_gtrans = u_gtrans / REAL( nzt - nzb )
105          v_gtrans = v_gtrans / REAL( nzt - nzb )
106       ELSE
107!
108!--       Averaging over the entire model domain.
109          uv_gtrans_l = 0.0
110          DO  i = nxl, nxr
111             DO  j = nys, nyn
112                DO  k = nzb+1, nzt
113                   uv_gtrans_l(1) = uv_gtrans_l(1) + u(k,j,i)
114                   uv_gtrans_l(2) = uv_gtrans_l(2) + v(k,j,i)
115                ENDDO
116             ENDDO
117          ENDDO
118          uv_gtrans_l = uv_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
119#if defined( __parallel )
120          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
121                              comm2d, ierr )
122          u_gtrans = uv_gtrans(1) / REAL( numprocs )
123          v_gtrans = uv_gtrans(2) / REAL( numprocs )
124#else
125          u_gtrans = uv_gtrans_l(1)
126          v_gtrans = uv_gtrans_l(2)
127#endif
128       ENDIF
129    ENDIF
130
131    IF ( .NOT. dt_fixed )  THEN
132!
133!--    Variable time step:
134!
135!--    For each component, compute the maximum time step according to the
136!--    FCL-criterion.
137       dt_u = dx / ( ABS( u_max - u_gtrans ) + 1.0E-10 )
138       dt_v = dy / ( ABS( v_max - v_gtrans ) + 1.0E-10 )
139       dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 )
140
141!
142!--    Compute time step according to the diffusion criterion.
143!--    First calculate minimum grid spacing which only depends on index k
144!--    Note: also at k=nzb+1 a full grid length is being assumed, although
145!--          in the Prandtl-layer friction term only dz/2 is used.
146!--          Experience from the old model seems to justify this.
147       dt_diff_l = 999999.0
148
149       DO  k = nzb+1, nzt
150           dxyz2_min(k) = MIN( dx2, dy2, dzu(k)*dzu(k) ) * 0.125
151       ENDDO
152
153!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
154!$OMP DO
155       DO  i = nxl, nxr
156          DO  j = nys, nyn
157             DO  k = nzb+1, nzt
158                value = dxyz2_min(k) / ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 )
159
160                dt_diff_l = MIN( value, dt_diff_l )
161             ENDDO
162          ENDDO
163       ENDDO
164!$OMP END PARALLEL
165#if defined( __parallel )
166       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, &
167                           ierr )
168#else
169       dt_diff = dt_diff_l
170#endif
171
172!
173!--    In case of non-cyclic lateral boundaries, the diffusion time step
174!--    may be further restricted by the lateral damping layer (damping only
175!--    along x and y)
176       IF ( bc_lr /= 'cyclic' )  THEN
177          dt_diff = MIN( dt_diff, 0.125 * dx2 / ( km_damp_max + 1E-20 ) )
178       ELSEIF ( bc_ns /= 'cyclic' )  THEN
179          dt_diff = MIN( dt_diff, 0.125 * dy2 / ( km_damp_max + 1E-20 ) )
180       ENDIF
181
182!
[316]183!--    Additional timestep criterion with plant canopies:
184!--    it is not allowed to extract more than the available momentum
185       IF ( plant_canopy ) THEN
[318]186
187          dt_plant_canopy_l = 0.0
188          DO  i = nxl, nxr
189             DO  j = nys, nyn
190                DO k = nzb+1, nzt
191                   dt_plant_canopy_u = cdc(k,j,i) * lad_u(k,j,i) *  &
192                                       SQRT(     u(k,j,i)**2     +  &
193                                             ( ( v(k,j,i-1)      +  &
194                                                 v(k,j,i)        +  &
195                                                 v(k,j+1,i)      +  &
196                                                 v(k,j+1,i-1) )     &
197                                               / 4.0 )**2        +  &
198                                             ( ( w(k-1,j,i-1)    +  &
199                                                 w(k-1,j,i)      +  &
200                                                 w(k,j,i-1)      +  &
201                                                 w(k,j,i) )         &
202                                                 / 4.0 )**2 ) 
203                   IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
204                      dt_plant_canopy_l = dt_plant_canopy_u 
205                   ENDIF
206                   dt_plant_canopy_v = cdc(k,j,i) * lad_v(k,j,i) *  &
207                                       SQRT( ( ( u(k,j-1,i)      +  &
208                                                 u(k,j-1,i+1)    +  &
209                                                 u(k,j,i)        +  &
210                                                 u(k,j,i+1) )       &
211                                               / 4.0 )**2        +  &
212                                                 v(k,j,i)**2     +  &
213                                             ( ( w(k-1,j-1,i)    +  &
214                                                 w(k-1,j,i)      +  &
215                                                 w(k,j-1,i)      +  &
216                                                 w(k,j,i) )         &
217                                                 / 4.0 )**2 ) 
218                   IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
219                      dt_plant_canopy_l = dt_plant_canopy_v
220                   ENDIF                   
221                   dt_plant_canopy_w = cdc(k,j,i) * lad_w(k,j,i) *  &
222                                       SQRT( ( ( u(k,j,i)        +  &
223                                                 u(k,j,i+1)      +  &
224                                                 u(k+1,j,i)      +  &
225                                                 u(k+1,j,i+1) )     &
226                                               / 4.0 )**2        +  &
227                                             ( ( v(k,j,i)        +  &
228                                                 v(k,j+1,i)      +  &
229                                                 v(k+1,j,i)      +  &
230                                                 v(k+1,j+1,i) )     &
231                                               / 4.0 )**2        +  &
232                                                 w(k,j,i)**2 )     
233                   IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
234                      dt_plant_canopy_l = dt_plant_canopy_w
235                   ENDIF
236                ENDDO
237             ENDDO
238          ENDDO 
239
240          IF ( dt_plant_canopy_l > 0.0 ) THEN
[320]241!
242!--          Invert dt_plant_canopy_l and apply a security timestep factor 0.1
[318]243             dt_plant_canopy_l = 0.1 / dt_plant_canopy_l
[320]244          ELSE
245!
246!--          In case of inhomogeneous plant canopy, some processors may have no
247!--          canopy at all. Then use dt_max as dummy instead.
248             dt_plant_canopy_l = dt_max
[318]249          ENDIF
[320]250
[316]251!
[318]252!--       Determine the global minumum
253#if defined( __parallel )
254          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL,  &
255                              MPI_MIN, comm2d, ierr )
256#else
257          dt_plant_canopy = dt_plant_canopy_l
258#endif
[316]259
260       ELSE
261!
262!--       Use dt_diff as dummy value to avoid extra IF branches further below
263          dt_plant_canopy = dt_diff
264
265       ENDIF
266
267!
268!--    The time step is the minimum of the 3-4 components and the diffusion time
[1]269!--    step minus a reduction to be on the safe side. Factor 0.5 is necessary
270!--    since the leap-frog scheme always progresses by 2 * delta t.
271!--    The user has to set the cfl_factor small enough to ensure that the
272!--    divergences do not become too large.
273!--    The time step must not exceed the maximum allowed value.
274       IF ( timestep_scheme(1:5) == 'runge' )  THEN
[316]275          dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
[1]276       ELSE
[316]277          dt_3d = cfl_factor * 0.5 *  &
278                  MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
[1]279       ENDIF
280       dt_3d = MIN( dt_3d, dt_max )
281
282!
283!--    Remember the restricting time step criterion for later output.
[316]284       IF ( MIN( dt_u, dt_v, dt_w ) < MIN( dt_diff, dt_plant_canopy ) )  THEN
[1]285          timestep_reason = 'A'
[316]286       ELSEIF ( dt_plant_canopy < dt_diff )  THEN
287          timestep_reason = 'C'
[1]288       ELSE
289          timestep_reason = 'D'
290       ENDIF
291
292!
293!--    Set flag if the time step becomes too small.
294       IF ( dt_3d < ( 0.00001 * dt_max ) )  THEN
295          stop_dt = .TRUE.
[108]296
[320]297          WRITE( message_string, * ) 'Time step has reached minimum limit.',  &
298               '&dt              = ', dt_3d, ' s  Simulation is terminated.', &
299               '&old_dt          = ', old_dt, ' s',                           &
300               '&dt_u            = ', dt_u, ' s',                             &
301               '&dt_v            = ', dt_v, ' s',                             &
302               '&dt_w            = ', dt_w, ' s',                             &
303               '&dt_diff         = ', dt_diff, ' s',                          &
304               '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
305               '&u_max   = ', u_max, ' m/s   k=', u_max_ijk(1),               &
306               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
307               '&v_max   = ', v_max, ' m/s   k=', v_max_ijk(1),               &
308               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
309               '&w_max   = ', w_max, ' m/s   k=', w_max_ijk(1),               &
310               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
[258]311          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
[108]312!
313!--       In case of coupled runs inform the remote model of the termination
314!--       and its reason, provided the remote model has not already been
315!--       informed of another termination reason (terminate_coupled > 0) before.
[222]316#if defined( __parallel )
[108]317          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
318             terminate_coupled = 2
319             CALL MPI_SENDRECV( &
[206]320                  terminate_coupled,        1, MPI_INTEGER, target_id,  0, &
321                  terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0, &
[108]322                  comm_inter, status, ierr )
323          ENDIF
[222]324#endif
[1]325       ENDIF
326
327!
328!--    Ensure a smooth value (two significant digits) of the timestep. For
329!--    other schemes than Runge-Kutta, the following restrictions appear:
330!--    The current timestep is only then changed, if the change relative to
331!--    its previous value exceeds +5 % or -2 %. In case of a timestep
332!--    reduction, at least 30 iterations have to be performed before a timestep
333!--    enlargement is permitted again.
334       percent_change = dt_3d / old_dt - 1.0
335       IF ( percent_change > 0.05  .OR.  percent_change < -0.02  .OR. &
336            timestep_scheme(1:5) == 'runge' )  THEN
337
338!
339!--       Time step enlargement by no more than 2 %.
340          IF ( percent_change > 0.0  .AND.  simulated_time /= 0.0  .AND. &
341               timestep_scheme(1:5) /= 'runge' )  THEN
342             dt_3d = 1.02 * old_dt
343          ENDIF
344
345!
346!--       A relatively smooth value of the time step is ensured by taking
347!--       only the first two significant digits.
348          div = 1000.0
349          DO  WHILE ( dt_3d < div )
350             div = div / 10.0
351          ENDDO
352          dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0
353
354!
355!--       Now the time step can be adjusted.
356          IF ( percent_change < 0.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
357          THEN
358!
359!--          Time step reduction.
360             old_dt = dt_3d
361             dt_changed = .TRUE.
362          ELSE
363!
364!--          For other timestep schemes , the time step is only enlarged
365!--          after at least 30 iterations since the previous time step
366!--          change or, of course, after model initialization.
367             IF ( current_timestep_number >= last_dt_change + 30  .OR. &
368                  simulated_time == 0.0 )  THEN
369                old_dt = dt_3d
370                dt_changed = .TRUE.
371             ELSE
372                dt_3d = old_dt
373                dt_changed = .FALSE.
374             ENDIF
375
376          ENDIF
377       ELSE
378!
379!--       No time step change since the difference is too small.
380          dt_3d = old_dt
381          dt_changed = .FALSE.
382       ENDIF
383
384       IF ( dt_changed )  last_dt_change = current_timestep_number
385
386    ENDIF
387    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
388
389 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.