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

Last change on this file since 2152 was 2131, checked in by raasch, 7 years ago

last commit documented

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