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

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

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 9.5 KB
Line 
1 SUBROUTINE timestep
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: timestep.f90 4 2007-02-13 11:33:16Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.21  2006/02/23 12:59:44  raasch
14! nt_anz renamed current_timestep_number
15!
16! Revision 1.1  1997/08/11 06:26:19  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Compute the time step under consideration of the FCL and diffusion criterion.
23!------------------------------------------------------------------------------!
24
25    USE arrays_3d
26    USE control_parameters
27    USE cpulog
28    USE grid_variables
29    USE indices
30    USE interfaces
31    USE pegrid
32    USE statistics
33
34    IMPLICIT NONE
35
36    INTEGER ::  i, j, k
37
38    REAL ::  div, dt_diff, dt_diff_l, dt_u, dt_v, dt_w, percent_change, &
39             u_gtrans_l, value, v_gtrans_l
40
41    REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
42    REAL, DIMENSION(nzb+1:nzt) ::  dxyz2_min
43
44
45    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
46
47!
48!-- Determine the maxima of the velocity components.
49    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, u, 'abs', &
50                         u_max, u_max_ijk )
51    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, v, 'abs', &
52                         v_max, v_max_ijk )
53    CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, w, 'abs', &
54                         w_max, w_max_ijk )
55
56!
57!-- In case maxima of the horizontal velocity components have been found at the
58!-- bottom boundary (k=nzb), the corresponding maximum at level k=1 is chosen
59!-- if the Dirichlet-boundary condition ('mirror') has been set. This is
60!-- necessary, because otherwise in case of Galilei-transform a far too large
61!-- velocity (having the respective opposite sign) would be used for the time
62!-- step determination (almost double the mean flow velocity).
63    IF ( ibc_uv_b == 0 )  THEN
64       IF ( u_max_ijk(1) == nzb )  THEN
65          u_max        = -u_max
66          u_max_ijk(1) = nzb + 1
67       ENDIF
68       IF ( v_max_ijk(1) == nzb )  THEN
69          v_max        = -v_max
70          v_max_ijk(1) = nzb + 1
71       ENDIF
72    ENDIF
73
74!
75!-- In case of Galilei-transform not using the geostrophic wind as translation
76!-- velocity, compute the volume-averaged horizontal velocity components, which
77!-- will then be subtracted from the horizontal wind for the time step and
78!-- horizontal advection routines.
79    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
80       IF ( flow_statistics_called )  THEN
81!
82!--       Horizontal averages already existent, just need to average them
83!--       vertically.
84          u_gtrans = 0.0
85          v_gtrans = 0.0
86          DO  k = nzb+1, nzt
87             u_gtrans = u_gtrans + hom(k,1,1,0)
88             v_gtrans = v_gtrans + hom(k,1,2,0)
89          ENDDO
90          u_gtrans = u_gtrans / REAL( nzt - nzb )
91          v_gtrans = v_gtrans / REAL( nzt - nzb )
92       ELSE
93!
94!--       Averaging over the entire model domain.
95          uv_gtrans_l = 0.0
96          DO  i = nxl, nxr
97             DO  j = nys, nyn
98                DO  k = nzb+1, nzt
99                   uv_gtrans_l(1) = uv_gtrans_l(1) + u(k,j,i)
100                   uv_gtrans_l(2) = uv_gtrans_l(2) + v(k,j,i)
101                ENDDO
102             ENDDO
103          ENDDO
104          uv_gtrans_l = uv_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
105#if defined( __parallel )
106          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
107                              comm2d, ierr )
108          u_gtrans = uv_gtrans(1) / REAL( numprocs )
109          v_gtrans = uv_gtrans(2) / REAL( numprocs )
110#else
111          u_gtrans = uv_gtrans_l(1)
112          v_gtrans = uv_gtrans_l(2)
113#endif
114       ENDIF
115    ENDIF
116
117    IF ( .NOT. dt_fixed )  THEN
118!
119!--    Variable time step:
120!
121!--    For each component, compute the maximum time step according to the
122!--    FCL-criterion.
123       dt_u = dx / ( ABS( u_max - u_gtrans ) + 1.0E-10 )
124       dt_v = dy / ( ABS( v_max - v_gtrans ) + 1.0E-10 )
125       dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 )
126
127!
128!--    Compute time step according to the diffusion criterion.
129!--    First calculate minimum grid spacing which only depends on index k
130!--    Note: also at k=nzb+1 a full grid length is being assumed, although
131!--          in the Prandtl-layer friction term only dz/2 is used.
132!--          Experience from the old model seems to justify this.
133       dt_diff_l = 999999.0
134
135       DO  k = nzb+1, nzt
136           dxyz2_min(k) = MIN( dx2, dy2, dzu(k)*dzu(k) ) * 0.125
137       ENDDO
138
139!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
140!$OMP DO
141       DO  i = nxl, nxr
142          DO  j = nys, nyn
143             DO  k = nzb+1, nzt
144                value = dxyz2_min(k) / ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 )
145
146                dt_diff_l = MIN( value, dt_diff_l )
147             ENDDO
148          ENDDO
149       ENDDO
150!$OMP END PARALLEL
151#if defined( __parallel )
152       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, &
153                           ierr )
154#else
155       dt_diff = dt_diff_l
156#endif
157
158!
159!--    In case of non-cyclic lateral boundaries, the diffusion time step
160!--    may be further restricted by the lateral damping layer (damping only
161!--    along x and y)
162       IF ( bc_lr /= 'cyclic' )  THEN
163          dt_diff = MIN( dt_diff, 0.125 * dx2 / ( km_damp_max + 1E-20 ) )
164       ELSEIF ( bc_ns /= 'cyclic' )  THEN
165          dt_diff = MIN( dt_diff, 0.125 * dy2 / ( km_damp_max + 1E-20 ) )
166       ENDIF
167
168!
169!--    The time step is the minimum of the 3 components and the diffusion time
170!--    step minus a reduction to be on the safe side. Factor 0.5 is necessary
171!--    since the leap-frog scheme always progresses by 2 * delta t.
172!--    The user has to set the cfl_factor small enough to ensure that the
173!--    divergences do not become too large.
174!--    The time step must not exceed the maximum allowed value.
175       IF ( timestep_scheme(1:5) == 'runge' )  THEN
176          dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w )
177       ELSE
178          dt_3d = cfl_factor * 0.5 * MIN( dt_diff, dt_u, dt_v, dt_w )
179       ENDIF
180       dt_3d = MIN( dt_3d, dt_max )
181
182!
183!--    Remember the restricting time step criterion for later output.
184       IF ( dt_diff > MIN( dt_u, dt_v, dt_w ) )  THEN
185          timestep_reason = 'A'
186       ELSE
187          timestep_reason = 'D'
188       ENDIF
189
190!
191!--    Set flag if the time step becomes too small.
192       IF ( dt_3d < ( 0.00001 * dt_max ) )  THEN
193          stop_dt = .TRUE.
194          IF ( myid == 0 )  THEN
195             PRINT*,'+++ time_step: Time step has reached minimum limit.'
196             PRINT*,'    dt      = ', dt_3d, ' s   Simulation is terminated.'
197             PRINT*,'    old_dt  = ', old_dt, ' s'
198             PRINT*,'    dt_u    = ', dt_u, ' s'
199             PRINT*,'    dt_v    = ', dt_v, ' s'
200             PRINT*,'    dt_w    = ', dt_w, ' s'
201             PRINT*,'    dt_diff = ', dt_diff, ' s'
202             PRINT*,'    u_max   = ', u_max, ' m/s   k=', u_max_ijk(1), &
203                    '  j=', u_max_ijk(2), '  i=', u_max_ijk(3)
204             PRINT*,'    v_max   = ', v_max, ' m/s   k=', v_max_ijk(1), &
205                    '  j=', v_max_ijk(2), '  i=', v_max_ijk(3)
206             PRINT*,'    w_max   = ', w_max, ' m/s   k=', w_max_ijk(1), &
207                    '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
208          ENDIF
209       ENDIF
210
211!
212!--    Ensure a smooth value (two significant digits) of the timestep. For
213!--    other schemes than Runge-Kutta, the following restrictions appear:
214!--    The current timestep is only then changed, if the change relative to
215!--    its previous value exceeds +5 % or -2 %. In case of a timestep
216!--    reduction, at least 30 iterations have to be performed before a timestep
217!--    enlargement is permitted again.
218       percent_change = dt_3d / old_dt - 1.0
219       IF ( percent_change > 0.05  .OR.  percent_change < -0.02  .OR. &
220            timestep_scheme(1:5) == 'runge' )  THEN
221
222!
223!--       Time step enlargement by no more than 2 %.
224          IF ( percent_change > 0.0  .AND.  simulated_time /= 0.0  .AND. &
225               timestep_scheme(1:5) /= 'runge' )  THEN
226             dt_3d = 1.02 * old_dt
227          ENDIF
228
229!
230!--       A relatively smooth value of the time step is ensured by taking
231!--       only the first two significant digits.
232          div = 1000.0
233          DO  WHILE ( dt_3d < div )
234             div = div / 10.0
235          ENDDO
236          dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0
237
238!
239!--       Now the time step can be adjusted.
240          IF ( percent_change < 0.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
241          THEN
242!
243!--          Time step reduction.
244             old_dt = dt_3d
245             dt_changed = .TRUE.
246          ELSE
247!
248!--          For other timestep schemes , the time step is only enlarged
249!--          after at least 30 iterations since the previous time step
250!--          change or, of course, after model initialization.
251             IF ( current_timestep_number >= last_dt_change + 30  .OR. &
252                  simulated_time == 0.0 )  THEN
253                old_dt = dt_3d
254                dt_changed = .TRUE.
255             ELSE
256                dt_3d = old_dt
257                dt_changed = .FALSE.
258             ENDIF
259
260          ENDIF
261       ELSE
262!
263!--       No time step change since the difference is too small.
264          dt_3d = old_dt
265          dt_changed = .FALSE.
266       ENDIF
267
268       IF ( dt_changed )  last_dt_change = current_timestep_number
269
270    ENDIF
271
272    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
273
274
275 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.