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

Last change on this file since 3637 was 3634, checked in by knoop, 6 years ago

OpenACC port for SPEC

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