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

Last change on this file since 2769 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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