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

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

bugfix: in case of nested runs the stop condition in case of too small timesteps is communicated to all parent/child processes

  • Property svn:keywords set to Id
File size: 13.8 KB
Line 
1!> @file timestep.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! bugfix: in case of nested runs the stop condition in case of too small
23! timesteps is communicated to all parent/child processes,
24! some formatting done
25!
26! Former revisions:
27! -----------------
28! $Id: timestep.f90 2130 2017-01-24 16:25:39Z raasch $
29!
30! 2118 2017-01-17 16:38:49Z raasch
31! OpenACC directives and related part of code removed
32!
33! 2000 2016-08-20 18:09:15Z knoop
34! Forced header and separation lines into 80 columns
35!
36! 1849 2016-04-08 11:33:18Z hoffmann
37! Adapted for modularization of microphysics
38!
39! 1682 2015-10-07 23:56:08Z knoop
40! Code annotations made doxygen readable
41!
42! 1484 2014-10-21 10:53:05Z kanani
43! Changes due to new module structure of the plant canopy model:
44!   calculations and parameters related to the plant canopy model removed
45!   (the limitation of the canopy drag, i.e. that the canopy drag itself should
46!   not change the sign of the velocity components, is now assured for in the
47!   calculation of the canopy tendency terms in subroutine plant_canopy_model)
48!
49! 1342 2014-03-26 17:04:47Z kanani
50! REAL constants defined as wp-kind
51!
52! 1322 2014-03-20 16:38:49Z raasch
53! REAL functions provided with KIND-attribute
54!
55! 1320 2014-03-20 08:40:49Z raasch
56! ONLY-attribute added to USE-statements,
57! kind-parameters added to all INTEGER and REAL declaration statements,
58! kinds are defined in new module kinds,
59! old module precision_kind is removed,
60! revision history before 2012 removed,
61! comment fields (!:) to be used for variable explanations added to
62! all variable declaration statements
63!
64! 1257 2013-11-08 15:18:40Z raasch
65! openacc porting
66! bugfix for calculation of advective timestep in case of vertically stretched
67! grids
68!
69! 1092 2013-02-02 11:24:22Z raasch
70! unused variables removed
71!
72! 1053 2012-11-13 17:11:03Z hoffmann
73! timestep is reduced in two-moment cloud scheme according to the maximum
74! terminal velocity of rain drops
75!
76! 1036 2012-10-22 13:43:42Z raasch
77! code put under GPL (PALM 3.9)
78!
79! 1001 2012-09-13 14:08:46Z raasch
80! all actions concerning leapfrog scheme removed
81!
82! 978 2012-08-09 08:28:32Z fricke
83! restriction of the outflow damping layer in the diffusion criterion removed
84!
85! 866 2012-03-28 06:44:41Z raasch
86! bugfix for timestep calculation in case of Galilei transformation,
87! special treatment in case of mirror velocity boundary condition removed
88!
89! Revision 1.1  1997/08/11 06:26:19  raasch
90! Initial revision
91!
92!
93! Description:
94! ------------
95!> Compute the time step under consideration of the FCL and diffusion criterion.
96!------------------------------------------------------------------------------!
97 SUBROUTINE timestep
98 
99
100    USE arrays_3d,                                                             &
101        ONLY:  dzu, dzw, kh, km, u, v, w
102
103    USE control_parameters,                                                    &
104        ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
105               galilei_transformation, old_dt, message_string,                 &
106               stop_dt, terminate_coupled, terminate_coupled_remote,           &
107               timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
108
109    USE cpulog,                                                                &
110        ONLY:  cpu_log, log_point
111
112    USE grid_variables,                                                        &
113        ONLY:  dx, dx2, dy, dy2
114
115    USE indices,                                                               &
116        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
117
118    USE interfaces
119
120    USE kinds
121
122    USE microphysics_mod,                                                      &
123        ONLY:  dt_precipitation
124
125    USE pegrid
126
127    USE pmc_interface,                                                         &
128        ONLY:  nested_run
129
130    USE statistics,                                                            &
131        ONLY:  flow_statistics_called, hom, u_max, u_max_ijk, v_max, v_max_ijk,&
132               w_max, w_max_ijk
133
134    IMPLICIT NONE
135
136    INTEGER(iwp) ::  i !<
137    INTEGER(iwp) ::  j !<
138    INTEGER(iwp) ::  k !<
139
140    LOGICAL ::  stop_dt_local !< local switch for controlling the time stepping
141
142    REAL(wp) ::  div               !<
143    REAL(wp) ::  dt_diff           !<
144    REAL(wp) ::  dt_diff_l         !<
145    REAL(wp) ::  dt_u              !<
146    REAL(wp) ::  dt_u_l            !<
147    REAL(wp) ::  dt_v              !<
148    REAL(wp) ::  dt_v_l            !<
149    REAL(wp) ::  dt_w              !<
150    REAL(wp) ::  dt_w_l            !<
151    REAL(wp) ::  u_gtrans_l        !<
152    REAL(wp) ::  u_max_l           !<
153    REAL(wp) ::  u_min_l           !<
154    REAL(wp) ::  value             !<
155    REAL(wp) ::  v_gtrans_l        !<
156    REAL(wp) ::  v_max_l           !<
157    REAL(wp) ::  v_min_l           !<
158    REAL(wp) ::  w_max_l           !<
159    REAL(wp) ::  w_min_l           !<
160 
161    REAL(wp), DIMENSION(2)         ::  uv_gtrans   !<
162    REAL(wp), DIMENSION(2)         ::  uv_gtrans_l !<
163    REAL(wp), DIMENSION(3)         ::  reduce      !<
164    REAL(wp), DIMENSION(3)         ::  reduce_l    !<
165    REAL(wp), DIMENSION(nzb+1:nzt) ::  dxyz2_min   !< 
166
167
168
169    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
170
171!
172!-- In case of Galilei-transform not using the geostrophic wind as translation
173!-- velocity, compute the volume-averaged horizontal velocity components, which
174!-- will then be subtracted from the horizontal wind for the time step and
175!-- horizontal advection routines.
176    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
177       IF ( flow_statistics_called )  THEN
178!
179!--       Horizontal averages already existent, just need to average them
180!--       vertically.
181          u_gtrans = 0.0_wp
182          v_gtrans = 0.0_wp
183          DO  k = nzb+1, nzt
184             u_gtrans = u_gtrans + hom(k,1,1,0)
185             v_gtrans = v_gtrans + hom(k,1,2,0)
186          ENDDO
187          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
188          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
189       ELSE
190!
191!--       Averaging over the entire model domain.
192          u_gtrans_l = 0.0_wp
193          v_gtrans_l = 0.0_wp
194          DO  i = nxl, nxr
195             DO  j = nys, nyn
196                DO  k = nzb+1, nzt
197                   u_gtrans_l = u_gtrans_l + u(k,j,i)
198                   v_gtrans_l = v_gtrans_l + v(k,j,i)
199                ENDDO
200             ENDDO
201          ENDDO
202          uv_gtrans_l(1) = u_gtrans_l /                                        &
203                           REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
204          uv_gtrans_l(2) = v_gtrans_l /                                        &
205                           REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
206#if defined( __parallel )
207          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
208          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM,    &
209                              comm2d, ierr )
210          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
211          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
212#else
213          u_gtrans = uv_gtrans_l(1)
214          v_gtrans = uv_gtrans_l(2)
215#endif
216       ENDIF
217    ENDIF
218
219!
220!-- Determine the maxima of the velocity components, including their
221!-- grid index positions.
222    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
223                         u_max, u_max_ijk )
224    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
225                         v_max, v_max_ijk )
226    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
227                         w_max, w_max_ijk )
228
229    IF ( .NOT. dt_fixed )  THEN
230!
231!--    Variable time step:
232!--    Calculate the maximum time step according to the CFL-criterion,
233!--    individually for each velocity component
234       dt_u_l = 999999.9_wp
235       dt_v_l = 999999.9_wp
236       dt_w_l = 999999.9_wp
237       DO  i = nxl, nxr
238          DO  j = nys, nyn
239             DO  k = nzb+1, nzt
240                dt_u_l = MIN( dt_u_l, ( dx     /                               &
241                                 ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
242                dt_v_l = MIN( dt_v_l, ( dy     /                               &
243                                 ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
244                dt_w_l = MIN( dt_w_l, ( dzu(k) /                               &
245                                 ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
246             ENDDO
247          ENDDO
248       ENDDO
249
250#if defined( __parallel )
251       reduce_l(1) = dt_u_l
252       reduce_l(2) = dt_v_l
253       reduce_l(3) = dt_w_l
254       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
255       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
256       dt_u = reduce(1)
257       dt_v = reduce(2)
258       dt_w = reduce(3)
259#else
260       dt_u = dt_u_l
261       dt_v = dt_v_l
262       dt_w = dt_w_l
263#endif
264
265!
266!--    Compute time step according to the diffusion criterion.
267!--    First calculate minimum grid spacing which only depends on index k
268       dt_diff_l = 999999.0_wp
269
270       DO  k = nzb+1, nzt
271           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
272       ENDDO
273
274       !$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
275       !$OMP DO
276       DO  i = nxl, nxr
277          DO  j = nys, nyn
278             DO  k = nzb+1, nzt
279                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
280                                    ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
281             ENDDO
282          ENDDO
283       ENDDO
284       !$OMP END PARALLEL
285#if defined( __parallel )
286       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
287       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d,   &
288                           ierr )
289#else
290       dt_diff = dt_diff_l
291#endif
292
293!
294!--    The time step is the minimum of the 3-4 components and the diffusion time
295!--    step minus a reduction (cfl_factor) to be on the safe side.
296!--    The time step must not exceed the maximum allowed value.
297       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
298       dt_3d = MIN( dt_3d, dt_max )
299
300!
301!--    Remember the restricting time step criterion for later output.
302       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
303          timestep_reason = 'A'
304       ELSE
305          timestep_reason = 'D'
306       ENDIF
307
308!
309!--    Set flag if the time step becomes too small.
310       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
311          stop_dt = .TRUE.
312
313          WRITE( message_string, * ) 'Time step has reached minimum limit.',   &
314               '&dt              = ', dt_3d, ' s  Simulation is terminated.',  &
315               '&old_dt          = ', old_dt, ' s',                            &
316               '&dt_u            = ', dt_u, ' s',                              &
317               '&dt_v            = ', dt_v, ' s',                              &
318               '&dt_w            = ', dt_w, ' s',                              &
319               '&dt_diff         = ', dt_diff, ' s',                           &
320               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),        &
321               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                     &
322               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),        &
323               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                     &
324               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),        &
325               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
326          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
327!
328!--       In case of coupled runs inform the remote model of the termination
329!--       and its reason, provided the remote model has not already been
330!--       informed of another termination reason (terminate_coupled > 0) before.
331#if defined( __parallel )
332          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
333             terminate_coupled = 2
334             IF ( myid == 0 )  THEN
335                CALL MPI_SENDRECV( &
336                     terminate_coupled,        1, MPI_INTEGER, target_id,  0,  &
337                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0,  &
338                     comm_inter, status, ierr )
339             ENDIF
340             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0,      &
341                             comm2d, ierr)
342          ENDIF
343#endif
344       ENDIF
345
346!
347!--    In case of nested runs all parent/child processes have to terminate if
348!--    one process has set the stop flag, i.e. they need to set the stop flag
349!--    too.
350       IF ( nested_run )  THEN
351          stop_dt_local = stop_dt
352          CALL MPI_ALLREDUCE( stop_dt_local, stop_dt, 1, MPI_LOGICAL, MPI_LOR, &
353                              MPI_COMM_WORLD, ierr )
354       ENDIF
355
356!
357!--    Ensure a smooth value (two significant digits) of the timestep.
358       div = 1000.0_wp
359       DO  WHILE ( dt_3d < div )
360          div = div / 10.0_wp
361       ENDDO
362       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
363
364!
365!--    Adjust the time step
366       old_dt = dt_3d
367
368    ENDIF
369
370    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
371
372 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.