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

Last change on this file since 707 was 707, checked in by raasch, 10 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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