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

Last change on this file since 3048 was 3046, checked in by Giersch, 6 years ago

Remaining error messages revised, comments extended

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