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

Last change on this file since 2347 was 2258, checked in by suehring, 8 years ago

Bugfix, add pre-preprocessor directives to enable non-parrallel mode

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