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

Last change on this file since 3516 was 3311, checked in by raasch, 5 years ago

Stokes drift is regarded in timestep calculation, check if ocean mode is used for invalid combinations

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