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

Last change on this file since 1969 was 1852, checked in by hoffmann, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 17.5 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 terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: timestep.f90 1852 2016-04-08 14:07:36Z suehring $
26!
27! 1849 2016-04-08 11:33:18Z hoffmann
28! Adapted for modularization of microphysics
29!
30! 1682 2015-10-07 23:56:08Z knoop
31! Code annotations made doxygen readable
32!
33! 1484 2014-10-21 10:53:05Z kanani
34! Changes due to new module structure of the plant canopy model:
35!   calculations and parameters related to the plant canopy model removed
36!   (the limitation of the canopy drag, i.e. that the canopy drag itself should
37!   not change the sign of the velocity components, is now assured for in the
38!   calculation of the canopy tendency terms in subroutine plant_canopy_model)
39!
40! 1342 2014-03-26 17:04:47Z kanani
41! REAL constants defined as wp-kind
42!
43! 1322 2014-03-20 16:38:49Z raasch
44! REAL functions provided with KIND-attribute
45!
46! 1320 2014-03-20 08:40:49Z raasch
47! ONLY-attribute added to USE-statements,
48! kind-parameters added to all INTEGER and REAL declaration statements,
49! kinds are defined in new module kinds,
50! old module precision_kind is removed,
51! revision history before 2012 removed,
52! comment fields (!:) to be used for variable explanations added to
53! all variable declaration statements
54!
55! 1257 2013-11-08 15:18:40Z raasch
56! openacc porting
57! bugfix for calculation of advective timestep in case of vertically stretched
58! grids
59!
60! 1092 2013-02-02 11:24:22Z raasch
61! unused variables removed
62!
63! 1053 2012-11-13 17:11:03Z hoffmann
64! timestep is reduced in two-moment cloud scheme according to the maximum
65! terminal velocity of rain drops
66!
67! 1036 2012-10-22 13:43:42Z raasch
68! code put under GPL (PALM 3.9)
69!
70! 1001 2012-09-13 14:08:46Z raasch
71! all actions concerning leapfrog scheme removed
72!
73! 978 2012-08-09 08:28:32Z fricke
74! restriction of the outflow damping layer in the diffusion criterion removed
75!
76! 866 2012-03-28 06:44:41Z raasch
77! bugfix for timestep calculation in case of Galilei transformation,
78! special treatment in case of mirror velocity boundary condition removed
79!
80! Revision 1.1  1997/08/11 06:26:19  raasch
81! Initial revision
82!
83!
84! Description:
85! ------------
86!> Compute the time step under consideration of the FCL and diffusion criterion.
87!------------------------------------------------------------------------------!
88 SUBROUTINE timestep
89 
90
91    USE arrays_3d,                                                             &
92        ONLY:  dzu, dzw, kh, km, u, v, w
93
94    USE control_parameters,                                                    &
95        ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
96               galilei_transformation, old_dt, message_string,                 &
97               stop_dt, terminate_coupled, terminate_coupled_remote,           &
98               timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
99
100    USE cpulog,                                                                &
101        ONLY:  cpu_log, log_point
102
103    USE grid_variables,                                                        &
104        ONLY:  dx, dx2, dy, dy2
105
106    USE indices,                                                               &
107        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
108
109    USE interfaces
110
111    USE kinds
112
113    USE microphysics_mod,                                                      &
114        ONLY:  dt_precipitation
115
116    USE pegrid
117
118    USE statistics,                                                            &
119        ONLY:  flow_statistics_called, hom, u_max, u_max_ijk, v_max, v_max_ijk,&
120               w_max, w_max_ijk
121
122    IMPLICIT NONE
123
124    INTEGER(iwp) ::  i !<
125    INTEGER(iwp) ::  j !<
126    INTEGER(iwp) ::  k !<
127
128    REAL(wp) ::  div               !<
129    REAL(wp) ::  dt_diff           !<
130    REAL(wp) ::  dt_diff_l         !<
131    REAL(wp) ::  dt_u              !<
132    REAL(wp) ::  dt_u_l            !<
133    REAL(wp) ::  dt_v              !<
134    REAL(wp) ::  dt_v_l            !<
135    REAL(wp) ::  dt_w              !<
136    REAL(wp) ::  dt_w_l            !<
137    REAL(wp) ::  u_gtrans_l        !<
138    REAL(wp) ::  u_max_l           !<
139    REAL(wp) ::  u_min_l           !<
140    REAL(wp) ::  value             !<
141    REAL(wp) ::  v_gtrans_l        !<
142    REAL(wp) ::  v_max_l           !<
143    REAL(wp) ::  v_min_l           !<
144    REAL(wp) ::  w_max_l           !<
145    REAL(wp) ::  w_min_l           !<
146 
147    REAL(wp), DIMENSION(2)         ::  uv_gtrans   !<
148    REAL(wp), DIMENSION(2)         ::  uv_gtrans_l !<
149    REAL(wp), DIMENSION(3)         ::  reduce      !<
150    REAL(wp), DIMENSION(3)         ::  reduce_l    !<
151    REAL(wp), DIMENSION(nzb+1:nzt) ::  dxyz2_min   !< 
152
153
154
155    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
156
157!
158!-- In case of Galilei-transform not using the geostrophic wind as translation
159!-- velocity, compute the volume-averaged horizontal velocity components, which
160!-- will then be subtracted from the horizontal wind for the time step and
161!-- horizontal advection routines.
162    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
163       IF ( flow_statistics_called )  THEN
164!
165!--       Horizontal averages already existent, just need to average them
166!--       vertically.
167          u_gtrans = 0.0_wp
168          v_gtrans = 0.0_wp
169          DO  k = nzb+1, nzt
170             u_gtrans = u_gtrans + hom(k,1,1,0)
171             v_gtrans = v_gtrans + hom(k,1,2,0)
172          ENDDO
173          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
174          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
175       ELSE
176!
177!--       Averaging over the entire model domain.
178          u_gtrans_l = 0.0_wp
179          v_gtrans_l = 0.0_wp
180          !$acc parallel present( u, v )
181          DO  i = nxl, nxr
182             DO  j = nys, nyn
183                DO  k = nzb+1, nzt
184                   u_gtrans_l = u_gtrans_l + u(k,j,i)
185                   v_gtrans_l = v_gtrans_l + v(k,j,i)
186                ENDDO
187             ENDDO
188          ENDDO
189          !$acc end parallel
190          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
191          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
192#if defined( __parallel )
193          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
194          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
195                              comm2d, ierr )
196          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
197          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
198#else
199          u_gtrans = uv_gtrans_l(1)
200          v_gtrans = uv_gtrans_l(2)
201#endif
202       ENDIF
203    ENDIF
204
205!
206!-- Determine the maxima of the velocity components, including their
207!-- grid index positions.
208#if defined( __openacc )
209    IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
210       u_max_l = -999999.9_wp
211       u_min_l =  999999.9_wp
212       v_max_l = -999999.9_wp
213       v_min_l =  999999.9_wp
214       w_max_l = -999999.9_wp
215       w_min_l =  999999.9_wp
216       !$acc parallel present( u, v, w )
217       DO  i = nxl, nxr
218          DO  j = nys, nyn
219             DO  k = nzb+1, nzt
220                u_max_l = MAX( u_max_l, u(k,j,i) )
221                u_min_l = MIN( u_min_l, u(k,j,i) )
222                v_max_l = MAX( v_max_l, v(k,j,i) )
223                v_min_l = MIN( v_min_l, v(k,j,i) )
224                w_max_l = MAX( w_max_l, w(k,j,i) )
225                w_min_l = MIN( w_min_l, w(k,j,i) )
226             ENDDO
227          ENDDO
228       ENDDO
229       !$acc end parallel
230#if defined( __parallel )
231       reduce_l(1) = u_max_l
232       reduce_l(2) = v_max_l
233       reduce_l(3) = w_max_l
234       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
235       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
236       u_max = reduce(1)
237       v_max = reduce(2)
238       w_max = reduce(3)
239       reduce_l(1) = u_min_l
240       reduce_l(2) = v_min_l
241       reduce_l(3) = w_min_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       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
245       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
246       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
247#else
248       IF ( ABS( u_min_l ) > u_max_l )  THEN
249          u_max = u_min_l
250       ELSE
251          u_max = u_max_l
252       ENDIF
253       IF ( ABS( v_min_l ) > v_max_l )  THEN
254          v_max = v_min_l
255       ELSE
256          v_max = v_max_l
257       ENDIF
258       IF ( ABS( w_min_l ) > w_max_l )  THEN
259          w_max = w_min_l
260       ELSE
261          w_max = w_max_l
262       ENDIF
263#endif
264    ENDIF
265#else
266    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
267                         u_max, u_max_ijk )
268    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
269                         v_max, v_max_ijk )
270    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
271                         w_max, w_max_ijk )
272#endif
273
274    IF ( .NOT. dt_fixed )  THEN
275#if defined( __openacc )
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       u_max_l = -999999.9_wp
284       u_min_l =  999999.9_wp
285       v_max_l = -999999.9_wp
286       v_min_l =  999999.9_wp
287       w_max_l = -999999.9_wp
288       w_min_l =  999999.9_wp
289       !$acc parallel loop collapse(3) present( u, v, w )
290       DO  i = nxl, nxr
291          DO  j = nys, nyn
292             DO  k = nzb+1, nzt
293                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
294                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
295                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
296                u_max_l = MAX( u_max_l, u(k,j,i) )
297                u_min_l = MIN( u_min_l, u(k,j,i) )
298                v_max_l = MAX( v_max_l, v(k,j,i) )
299                v_min_l = MIN( v_min_l, v(k,j,i) )
300                w_max_l = MAX( w_max_l, w(k,j,i) )
301                w_min_l = MIN( w_min_l, w(k,j,i) )
302             ENDDO
303          ENDDO
304       ENDDO
305       !$acc end parallel
306
307#if defined( __parallel )
308       reduce_l(1) = dt_u_l
309       reduce_l(2) = dt_v_l
310       reduce_l(3) = dt_w_l
311       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
312       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
313       dt_u = reduce(1)
314       dt_v = reduce(2)
315       dt_w = reduce(3)
316
317       reduce_l(1) = u_max_l
318       reduce_l(2) = v_max_l
319       reduce_l(3) = w_max_l
320       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
321       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
322       u_max = reduce(1)
323       v_max = reduce(2)
324       w_max = reduce(3)
325       reduce_l(1) = u_min_l
326       reduce_l(2) = v_min_l
327       reduce_l(3) = w_min_l
328       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
329       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
330       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
331       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
332       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
333#else
334       dt_u = dt_u_l
335       dt_v = dt_v_l
336       dt_w = dt_w_l
337
338       IF ( ABS( u_min_l ) > u_max_l )  THEN
339          u_max = u_min_l
340       ELSE
341          u_max = u_max_l
342       ENDIF
343       IF ( ABS( v_min_l ) > v_max_l )  THEN
344          v_max = v_min_l
345       ELSE
346          v_max = v_max_l
347       ENDIF
348       IF ( ABS( w_min_l ) > w_max_l )  THEN
349          w_max = w_min_l
350       ELSE
351          w_max = w_max_l
352       ENDIF
353#endif
354
355#else
356!
357!--    Variable time step:
358!--    Calculate the maximum time step according to the CFL-criterion,
359!--    individually for each velocity component
360       dt_u_l = 999999.9_wp
361       dt_v_l = 999999.9_wp
362       dt_w_l = 999999.9_wp
363       DO  i = nxl, nxr
364          DO  j = nys, nyn
365             DO  k = nzb+1, nzt
366                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
367                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
368                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
369             ENDDO
370          ENDDO
371       ENDDO
372
373#if defined( __parallel )
374       reduce_l(1) = dt_u_l
375       reduce_l(2) = dt_v_l
376       reduce_l(3) = dt_w_l
377       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
378       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
379       dt_u = reduce(1)
380       dt_v = reduce(2)
381       dt_w = reduce(3)
382#else
383       dt_u = dt_u_l
384       dt_v = dt_v_l
385       dt_w = dt_w_l
386#endif
387
388#endif
389
390!
391!--    Compute time step according to the diffusion criterion.
392!--    First calculate minimum grid spacing which only depends on index k
393!--    Note: also at k=nzb+1 a full grid length is being assumed, although
394!--          in the Prandtl-layer friction term only dz/2 is used.
395!--          Experience from the old model seems to justify this.
396       dt_diff_l = 999999.0_wp
397
398       DO  k = nzb+1, nzt
399           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
400       ENDDO
401
402!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
403!$OMP DO
404       !$acc parallel loop collapse(3) present( kh, km )
405       DO  i = nxl, nxr
406          DO  j = nys, nyn
407             DO  k = nzb+1, nzt
408                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
409                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
410             ENDDO
411          ENDDO
412       ENDDO
413       !$acc end parallel
414!$OMP END PARALLEL
415#if defined( __parallel )
416       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
417       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, &
418                           ierr )
419#else
420       dt_diff = dt_diff_l
421#endif
422
423!
424!--    The time step is the minimum of the 3-4 components and the diffusion time
425!--    step minus a reduction (cfl_factor) to be on the safe side.
426!--    The time step must not exceed the maximum allowed value.
427       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w,   &
428                                 dt_precipitation )
429       dt_3d = MIN( dt_3d, dt_max )
430
431!
432!--    Remember the restricting time step criterion for later output.
433       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
434          timestep_reason = 'A'
435       ELSE
436          timestep_reason = 'D'
437       ENDIF
438
439!
440!--    Set flag if the time step becomes too small.
441       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
442          stop_dt = .TRUE.
443
444          WRITE( message_string, * ) 'Time step has reached minimum limit.',  &
445               '&dt              = ', dt_3d, ' s  Simulation is terminated.', &
446               '&old_dt          = ', old_dt, ' s',                           &
447               '&dt_u            = ', dt_u, ' s',                             &
448               '&dt_v            = ', dt_v, ' s',                             &
449               '&dt_w            = ', dt_w, ' s',                             &
450               '&dt_diff         = ', dt_diff, ' s',                          &
451               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
452               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
453               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),       &
454               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
455               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),       &
456               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
457          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
458!
459!--       In case of coupled runs inform the remote model of the termination
460!--       and its reason, provided the remote model has not already been
461!--       informed of another termination reason (terminate_coupled > 0) before.
462#if defined( __parallel )
463          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
464             terminate_coupled = 2
465             IF ( myid == 0 ) THEN
466                CALL MPI_SENDRECV( &
467                     terminate_coupled,        1, MPI_INTEGER, target_id,  0, &
468                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0, &
469                     comm_inter, status, ierr )
470             ENDIF
471             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
472          ENDIF
473#endif
474       ENDIF
475
476!
477!--    Ensure a smooth value (two significant digits) of the timestep.
478       div = 1000.0_wp
479       DO  WHILE ( dt_3d < div )
480          div = div / 10.0_wp
481       ENDDO
482       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
483
484!
485!--    Adjust the time step
486       old_dt = dt_3d
487
488    ENDIF
489
490    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
491
492 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.