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

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

implementation of dynamic sgs model

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