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

Last change on this file since 4124 was 4101, checked in by gronemeier, 5 years ago

timestep.f90:

  • consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
  • in RANS mode, instead of considering each wind component individually use the wind speed of 3d wind vector in CFL criterion
  • do not limit the increase of dt based on its previous value in RANS mode

other:

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