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

Last change on this file since 3058 was 3049, checked in by Giersch, 7 years ago

Revision history corrected

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