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

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

all OpenACC directives and related parts removed from the code

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