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

Last change on this file since 2481 was 2365, checked in by kanani, 7 years ago

Vertical nesting implemented (SadiqHuq?)

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