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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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