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

Last change on this file since 3294 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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