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

Last change on this file since 3819 was 3658, checked in by knoop, 6 years ago

OpenACC: flow_statistics partly ported to GPU

  • Property svn:keywords set to Id
File size: 17.2 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: timestep.f90 3658 2019-01-07 20:28:54Z hellstea $
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    !$ACC UPDATE &
215    !$ACC HOST(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
216    !$ACC HOST(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
217    !$ACC HOST(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
218    !$ACC HOST(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
219    !$ACC HOST(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
220
221!
222!--    Save former time step as reference
223       old_dt = dt_3d
224
225!
226!-- In case of Galilei-transform not using the geostrophic wind as translation
227!-- velocity, compute the volume-averaged horizontal velocity components, which
228!-- will then be subtracted from the horizontal wind for the time step and
229!-- horizontal advection routines.
230    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
231       IF ( flow_statistics_called )  THEN
232!
233!--       Horizontal averages already existent, just need to average them
234!--       vertically.
235          u_gtrans = 0.0_wp
236          v_gtrans = 0.0_wp
237          DO  k = nzb+1, nzt
238             u_gtrans = u_gtrans + hom(k,1,1,0)
239             v_gtrans = v_gtrans + hom(k,1,2,0)
240          ENDDO
241          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
242          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
243       ELSE
244!
245!--       Averaging over the entire model domain.
246          u_gtrans_l = 0.0_wp
247          v_gtrans_l = 0.0_wp
248          DO  i = nxl, nxr
249             DO  j = nys, nyn
250                DO  k = nzb+1, nzt
251                   u_gtrans_l = u_gtrans_l + u(k,j,i)
252                   v_gtrans_l = v_gtrans_l + v(k,j,i)
253                ENDDO
254             ENDDO
255          ENDDO
256          uv_gtrans_l(1) = u_gtrans_l /                                        &
257                           REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
258          uv_gtrans_l(2) = v_gtrans_l /                                        &
259                           REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
260#if defined( __parallel )
261          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
262          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM,    &
263                              comm2d, ierr )
264          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
265          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
266#else
267          u_gtrans = uv_gtrans_l(1)
268          v_gtrans = uv_gtrans_l(2)
269#endif
270       ENDIF
271    ENDIF
272
273!
274!-- Determine the maxima of the velocity components, including their
275!-- grid index positions.
276    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
277                         u_max, u_max_ijk )
278    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
279                         v_max, v_max_ijk )
280    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
281                         w_max, w_max_ijk )
282
283    IF ( .NOT. dt_fixed )  THEN
284!
285!--    Variable time step:
286!--    Calculate the maximum time step according to the CFL-criterion,
287!--    individually for each velocity component
288       dt_u_l = 999999.9_wp
289       dt_v_l = 999999.9_wp
290       dt_w_l = 999999.9_wp
291       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
292       !$ACC COPY(dt_u_l, dt_v_l, dt_w_l, u_stokes_zu, v_stokes_zu) &
293       !$ACC REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l) &
294       !$ACC PRESENT(u, v, w, dzu)
295       DO  i = nxl, nxr
296          DO  j = nys, nyn
297             DO  k = nzb+1, nzt
298                dt_u_l = MIN( dt_u_l, ( dx     /                               &
299                                 ( ABS( u(k,j,i) - u_gtrans + u_stokes_zu(k) ) &
300                                   + 1.0E-10_wp ) ) )
301                dt_v_l = MIN( dt_v_l, ( dy     /                               &
302                                 ( ABS( v(k,j,i) - v_gtrans + v_stokes_zu(k) ) &
303                                   + 1.0E-10_wp ) ) )
304                dt_w_l = MIN( dt_w_l, ( dzu(k) /                               &
305                                 ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
306             ENDDO
307          ENDDO
308       ENDDO
309
310#if defined( __parallel )
311       reduce_l(1) = dt_u_l
312       reduce_l(2) = dt_v_l
313       reduce_l(3) = dt_w_l
314       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
315       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
316       dt_u = reduce(1)
317       dt_v = reduce(2)
318       dt_w = reduce(3)
319#else
320       dt_u = dt_u_l
321       dt_v = dt_v_l
322       dt_w = dt_w_l
323#endif
324
325!
326!--    Compute time step according to the diffusion criterion.
327!--    First calculate minimum grid spacing which only depends on index k.
328!--    When using the dynamic subgrid model, negative km are possible.
329       dt_diff_l = 999999.0_wp
330
331       !$ACC PARALLEL LOOP PRESENT(dxyz2_min, dzw)
332       DO  k = nzb+1, nzt
333           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
334       ENDDO
335
336       !$OMP PARALLEL private(i,j,k) reduction(MIN: dt_diff_l)
337       !$OMP DO
338       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
339       !$ACC COPY(dt_diff_l) REDUCTION(MIN: dt_diff_l) &
340       !$ACC PRESENT(dxyz2_min, kh, km)
341       DO  i = nxl, nxr
342          DO  j = nys, nyn
343             DO  k = nzb+1, nzt
344                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
345                                    ( MAX( kh(k,j,i), ABS( km(k,j,i) ) )       &
346                                      + 1E-20_wp ) )
347             ENDDO
348          ENDDO
349       ENDDO
350       !$OMP END PARALLEL
351#if defined( __parallel )
352       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
353       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d,   &
354                           ierr )
355#else
356       dt_diff = dt_diff_l
357#endif
358
359!
360!--    The time step is the minimum of the 3-4 components and the diffusion time
361!--    step minus a reduction (cfl_factor) to be on the safe side.
362!--    The time step must not exceed the maximum allowed value.
363       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
364       dt_3d = MIN( dt_3d, dt_max )
365!
366!--    In RANS mode, the time step must not increase by more than a factor of 2
367       IF ( rans_mode )  dt_3d = MIN( dt_3d, dt_max, 2.0_wp * old_dt )
368
369!
370!--    Remember the restricting time step criterion for later output.
371       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
372          timestep_reason = 'A'
373       ELSE
374          timestep_reason = 'D'
375       ENDIF
376
377!
378!--    Set flag if the time step becomes too small.
379       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
380          stop_dt = .TRUE.
381
382!
383!--       Determine the maxima of the diffusion coefficients, including their
384!--       grid index positions.
385          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, km, 'abs',  &
386                               0.0_wp, km_max, km_max_ijk )
387          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, kh, 'abs',  &
388                               0.0_wp, kh_max, kh_max_ijk )
389
390          WRITE( message_string, * ) 'Time step has reached minimum limit.',   &
391               '&dt              = ', dt_3d, ' s  Simulation is terminated.',  &
392               '&old_dt          = ', old_dt, ' s',                            &
393               '&dt_u            = ', dt_u, ' s',                              &
394               '&dt_v            = ', dt_v, ' s',                              &
395               '&dt_w            = ', dt_w, ' s',                              &
396               '&dt_diff         = ', dt_diff, ' s',                           &
397               '&u_max           = ', u_max, ' m/s    k=', u_max_ijk(1),       &
398               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                     &
399               '&v_max           = ', v_max, ' m/s    k=', v_max_ijk(1),       &
400               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                     &
401               '&w_max           = ', w_max, ' m/s    k=', w_max_ijk(1),       &
402               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3),                     &
403               '&km_max          = ', km_max, ' m2/s2  k=', km_max_ijk(1),     &
404               '  j=', km_max_ijk(2), '  i=', km_max_ijk(3),                   &
405               '&kh_max          = ', kh_max, ' m2/s2  k=', kh_max_ijk(1),     &
406                '  j=', kh_max_ijk(2), '  i=', kh_max_ijk(3)
407          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
408!
409!--       In case of coupled runs inform the remote model of the termination
410!--       and its reason, provided the remote model has not already been
411!--       informed of another termination reason (terminate_coupled > 0) before.
412#if defined( __parallel )
413          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
414             terminate_coupled = 2
415             IF ( myid == 0 )  THEN
416                CALL MPI_SENDRECV( &
417                     terminate_coupled,        1, MPI_INTEGER, target_id,  0,  &
418                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0,  &
419                     comm_inter, status, ierr )
420             ENDIF
421             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0,      &
422                             comm2d, ierr)
423          ENDIF
424#endif
425       ENDIF
426
427!
428!--    In case of nested runs all parent/child processes have to terminate if
429!--    one process has set the stop flag, i.e. they need to set the stop flag
430!--    too.
431       IF ( nested_run )  THEN
432          stop_dt_local = stop_dt
433#if defined( __parallel )
434          CALL MPI_ALLREDUCE( stop_dt_local, stop_dt, 1, MPI_LOGICAL, MPI_LOR, &
435                              MPI_COMM_WORLD, ierr )
436#endif
437       ENDIF
438
439!
440!--    Ensure a smooth value (two significant digits) of the timestep.
441       div = 1000.0_wp
442       DO  WHILE ( dt_3d < div )
443          div = div / 10.0_wp
444       ENDDO
445       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
446
447    ENDIF
448
449!
450!-- Vertical nesting: coarse and fine grid timestep has to be identical   
451    IF ( vnested )  CALL vnest_timestep_sync
452
453    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
454
455 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.