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

Last change on this file since 2698 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

  • Property svn:keywords set to Id
File size: 14.4 KB
RevLine 
[1682]1!> @file timestep.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[2101]17! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[258]20! Current revisions:
[866]21! ------------------
[1485]22!
[2258]23!
[1485]24! Former revisions:
25! -----------------
26! $Id: timestep.f90 2696 2017-12-14 17:12:51Z suehring $
[2365]27! Vertical grid nesting: Sync fine and coarse grid timestep (SadiqHuq)
28!
29! 2258 2017-06-08 07:55:13Z suehring
[2258]30! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
[1485]31!
[2258]32! 2168 2017-03-06 13:08:38Z suehring
33!
[2131]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!
[2119]39! 2118 2017-01-17 16:38:49Z raasch
40! OpenACC directives and related part of code removed
41!
[2001]42! 2000 2016-08-20 18:09:15Z knoop
43! Forced header and separation lines into 80 columns
44!
[1851]45! 1849 2016-04-08 11:33:18Z hoffmann
46! Adapted for modularization of microphysics
[1852]47!
[1683]48! 1682 2015-10-07 23:56:08Z knoop
49! Code annotations made doxygen readable
50!
[1485]51! 1484 2014-10-21 10:53:05Z kanani
[1484]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)
[1343]57!
58! 1342 2014-03-26 17:04:47Z kanani
59! REAL constants defined as wp-kind
60!
[1323]61! 1322 2014-03-20 16:38:49Z raasch
62! REAL functions provided with KIND-attribute
63!
[1321]64! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1321]72!
[1258]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!
[1093]78! 1092 2013-02-02 11:24:22Z raasch
79! unused variables removed
80!
[1054]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!
[1037]85! 1036 2012-10-22 13:43:42Z raasch
86! code put under GPL (PALM 3.9)
87!
[1002]88! 1001 2012-09-13 14:08:46Z raasch
89! all actions concerning leapfrog scheme removed
90!
[979]91! 978 2012-08-09 08:28:32Z fricke
92! restriction of the outflow damping layer in the diffusion criterion removed
93!
[867]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!
[1]98! Revision 1.1  1997/08/11 06:26:19  raasch
99! Initial revision
100!
101!
102! Description:
103! ------------
[1682]104!> Compute the time step under consideration of the FCL and diffusion criterion.
[1]105!------------------------------------------------------------------------------!
[1682]106 SUBROUTINE timestep
107 
[1]108
[1320]109    USE arrays_3d,                                                             &
[1484]110        ONLY:  dzu, dzw, kh, km, u, v, w
[1320]111
112    USE control_parameters,                                                    &
113        ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
[1484]114               galilei_transformation, old_dt, message_string,                 &
[1320]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
[1]127    USE interfaces
[1320]128
129    USE kinds
130
[1849]131    USE microphysics_mod,                                                      &
132        ONLY:  dt_precipitation
133
[1]134    USE pegrid
135
[2130]136    USE pmc_interface,                                                         &
137        ONLY:  nested_run
138
[1320]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
[2365]143    USE vertical_nesting_mod,                                                  &
144        ONLY:  vnested, vnest_timestep_sync
145
[1]146    IMPLICIT NONE
147
[1682]148    INTEGER(iwp) ::  i !<
149    INTEGER(iwp) ::  j !<
150    INTEGER(iwp) ::  k !<
[1]151
[2130]152    LOGICAL ::  stop_dt_local !< local switch for controlling the time stepping
153
[1682]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           !<
[1320]172 
[1682]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   !< 
[1]178
179
[667]180
[1]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.
[1342]193          u_gtrans = 0.0_wp
194          v_gtrans = 0.0_wp
[1]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
[1322]199          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
200          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
[1]201       ELSE
202!
203!--       Averaging over the entire model domain.
[1342]204          u_gtrans_l = 0.0_wp
205          v_gtrans_l = 0.0_wp
[1]206          DO  i = nxl, nxr
207             DO  j = nys, nyn
208                DO  k = nzb+1, nzt
[1257]209                   u_gtrans_l = u_gtrans_l + u(k,j,i)
210                   v_gtrans_l = v_gtrans_l + v(k,j,i)
[1]211                ENDDO
212             ENDDO
213          ENDDO
[2130]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 )
[1]218#if defined( __parallel )
[622]219          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[2130]220          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM,    &
[1]221                              comm2d, ierr )
[1322]222          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
223          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
[1]224#else
225          u_gtrans = uv_gtrans_l(1)
226          v_gtrans = uv_gtrans_l(2)
227#endif
228       ENDIF
229    ENDIF
230
[866]231!
[1257]232!-- Determine the maxima of the velocity components, including their
233!-- grid index positions.
[1320]234    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
[866]235                         u_max, u_max_ijk )
[1320]236    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
[866]237                         v_max, v_max_ijk )
[1320]238    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
[866]239                         w_max, w_max_ijk )
240
[1257]241    IF ( .NOT. dt_fixed )  THEN
[866]242!
[1257]243!--    Variable time step:
244!--    Calculate the maximum time step according to the CFL-criterion,
245!--    individually for each velocity component
[1342]246       dt_u_l = 999999.9_wp
247       dt_v_l = 999999.9_wp
248       dt_w_l = 999999.9_wp
[1257]249       DO  i = nxl, nxr
250          DO  j = nys, nyn
251             DO  k = nzb+1, nzt
[2130]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 ) ) )
[1257]258             ENDDO
259          ENDDO
260       ENDDO
[1]261
[1257]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
[1]277!
278!--    Compute time step according to the diffusion criterion.
279!--    First calculate minimum grid spacing which only depends on index k
[1342]280       dt_diff_l = 999999.0_wp
[1]281
282       DO  k = nzb+1, nzt
[1342]283           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
[1]284       ENDDO
285
[2118]286       !$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
287       !$OMP DO
[1]288       DO  i = nxl, nxr
289          DO  j = nys, nyn
290             DO  k = nzb+1, nzt
[2130]291                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
292                                    ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
[1]293             ENDDO
294          ENDDO
295       ENDDO
[2118]296       !$OMP END PARALLEL
[1]297#if defined( __parallel )
[622]298       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[2130]299       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d,   &
[1]300                           ierr )
301#else
302       dt_diff = dt_diff_l
303#endif
304
305!
[316]306!--    The time step is the minimum of the 3-4 components and the diffusion time
[1001]307!--    step minus a reduction (cfl_factor) to be on the safe side.
[1]308!--    The time step must not exceed the maximum allowed value.
[2130]309       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
[1]310       dt_3d = MIN( dt_3d, dt_max )
311
312!
313!--    Remember the restricting time step criterion for later output.
[1484]314       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
[1]315          timestep_reason = 'A'
316       ELSE
317          timestep_reason = 'D'
318       ENDIF
319
320!
321!--    Set flag if the time step becomes too small.
[1342]322       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
[1]323          stop_dt = .TRUE.
[108]324
[2130]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),        &
[320]337               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
[258]338          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
[108]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.
[222]343#if defined( __parallel )
[108]344          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
345             terminate_coupled = 2
[2130]346             IF ( myid == 0 )  THEN
[667]347                CALL MPI_SENDRECV( &
[2130]348                     terminate_coupled,        1, MPI_INTEGER, target_id,  0,  &
349                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0,  &
[667]350                     comm_inter, status, ierr )
351             ENDIF
[2130]352             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0,      &
353                             comm2d, ierr)
[108]354          ENDIF
[222]355#endif
[1]356       ENDIF
357
358!
[2130]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
[2258]364#if defined( __parallel )
[2130]365          CALL MPI_ALLREDUCE( stop_dt_local, stop_dt, 1, MPI_LOGICAL, MPI_LOR, &
366                              MPI_COMM_WORLD, ierr )
[2258]367#endif
[2130]368       ENDIF
369
370!
[1001]371!--    Ensure a smooth value (two significant digits) of the timestep.
[1342]372       div = 1000.0_wp
[1001]373       DO  WHILE ( dt_3d < div )
[1342]374          div = div / 10.0_wp
[1001]375       ENDDO
[1342]376       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
[1]377
378!
[1001]379!--    Adjust the time step
380       old_dt = dt_3d
[1]381
[1001]382    ENDIF
[1]383
[2365]384!
385!-- Vertical nesting: coarse and fine grid timestep has to be identical   
386    IF ( vnested )  CALL vnest_timestep_sync
387
[1]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.