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

Last change on this file since 3084 was 3084, checked in by gronemeier, 6 years ago

bugfix: limit increase of dt_3d only in case of RANS mode

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