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

Last change on this file since 2119 was 2119, checked in by raasch, 4 years ago

last commit documented

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