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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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