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

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

merge with branch rans

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