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

Last change on this file since 2089 was 2001, checked in by knoop, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 17.6 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!
23!
24! Former revisions:
25! -----------------
26! $Id: timestep.f90 2001 2016-08-20 18:41:22Z suehring $
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          !$acc parallel present( u, v )
185          DO  i = nxl, nxr
186             DO  j = nys, nyn
187                DO  k = nzb+1, nzt
188                   u_gtrans_l = u_gtrans_l + u(k,j,i)
189                   v_gtrans_l = v_gtrans_l + v(k,j,i)
190                ENDDO
191             ENDDO
192          ENDDO
193          !$acc end parallel
194          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
195          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
196#if defined( __parallel )
197          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
198          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
199                              comm2d, ierr )
200          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
201          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
202#else
203          u_gtrans = uv_gtrans_l(1)
204          v_gtrans = uv_gtrans_l(2)
205#endif
206       ENDIF
207    ENDIF
208
209!
210!-- Determine the maxima of the velocity components, including their
211!-- grid index positions.
212#if defined( __openacc )
213    IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
214       u_max_l = -999999.9_wp
215       u_min_l =  999999.9_wp
216       v_max_l = -999999.9_wp
217       v_min_l =  999999.9_wp
218       w_max_l = -999999.9_wp
219       w_min_l =  999999.9_wp
220       !$acc parallel present( u, v, w )
221       DO  i = nxl, nxr
222          DO  j = nys, nyn
223             DO  k = nzb+1, nzt
224                u_max_l = MAX( u_max_l, u(k,j,i) )
225                u_min_l = MIN( u_min_l, u(k,j,i) )
226                v_max_l = MAX( v_max_l, v(k,j,i) )
227                v_min_l = MIN( v_min_l, v(k,j,i) )
228                w_max_l = MAX( w_max_l, w(k,j,i) )
229                w_min_l = MIN( w_min_l, w(k,j,i) )
230             ENDDO
231          ENDDO
232       ENDDO
233       !$acc end parallel
234#if defined( __parallel )
235       reduce_l(1) = u_max_l
236       reduce_l(2) = v_max_l
237       reduce_l(3) = w_max_l
238       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
239       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
240       u_max = reduce(1)
241       v_max = reduce(2)
242       w_max = reduce(3)
243       reduce_l(1) = u_min_l
244       reduce_l(2) = v_min_l
245       reduce_l(3) = w_min_l
246       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
247       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
248       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
249       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
250       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
251#else
252       IF ( ABS( u_min_l ) > u_max_l )  THEN
253          u_max = u_min_l
254       ELSE
255          u_max = u_max_l
256       ENDIF
257       IF ( ABS( v_min_l ) > v_max_l )  THEN
258          v_max = v_min_l
259       ELSE
260          v_max = v_max_l
261       ENDIF
262       IF ( ABS( w_min_l ) > w_max_l )  THEN
263          w_max = w_min_l
264       ELSE
265          w_max = w_max_l
266       ENDIF
267#endif
268    ENDIF
269#else
270    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
271                         u_max, u_max_ijk )
272    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
273                         v_max, v_max_ijk )
274    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
275                         w_max, w_max_ijk )
276#endif
277
278    IF ( .NOT. dt_fixed )  THEN
279#if defined( __openacc )
280!
281!--    Variable time step:
282!--    Calculate the maximum time step according to the CFL-criterion,
283!--    individually for each velocity component
284       dt_u_l  =  999999.9_wp
285       dt_v_l  =  999999.9_wp
286       dt_w_l  =  999999.9_wp
287       u_max_l = -999999.9_wp
288       u_min_l =  999999.9_wp
289       v_max_l = -999999.9_wp
290       v_min_l =  999999.9_wp
291       w_max_l = -999999.9_wp
292       w_min_l =  999999.9_wp
293       !$acc parallel loop collapse(3) present( u, v, w )
294       DO  i = nxl, nxr
295          DO  j = nys, nyn
296             DO  k = nzb+1, nzt
297                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
298                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
299                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
300                u_max_l = MAX( u_max_l, u(k,j,i) )
301                u_min_l = MIN( u_min_l, u(k,j,i) )
302                v_max_l = MAX( v_max_l, v(k,j,i) )
303                v_min_l = MIN( v_min_l, v(k,j,i) )
304                w_max_l = MAX( w_max_l, w(k,j,i) )
305                w_min_l = MIN( w_min_l, w(k,j,i) )
306             ENDDO
307          ENDDO
308       ENDDO
309       !$acc end parallel
310
311#if defined( __parallel )
312       reduce_l(1) = dt_u_l
313       reduce_l(2) = dt_v_l
314       reduce_l(3) = dt_w_l
315       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
316       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
317       dt_u = reduce(1)
318       dt_v = reduce(2)
319       dt_w = reduce(3)
320
321       reduce_l(1) = u_max_l
322       reduce_l(2) = v_max_l
323       reduce_l(3) = w_max_l
324       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
325       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
326       u_max = reduce(1)
327       v_max = reduce(2)
328       w_max = reduce(3)
329       reduce_l(1) = u_min_l
330       reduce_l(2) = v_min_l
331       reduce_l(3) = w_min_l
332       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
333       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
334       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
335       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
336       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
337#else
338       dt_u = dt_u_l
339       dt_v = dt_v_l
340       dt_w = dt_w_l
341
342       IF ( ABS( u_min_l ) > u_max_l )  THEN
343          u_max = u_min_l
344       ELSE
345          u_max = u_max_l
346       ENDIF
347       IF ( ABS( v_min_l ) > v_max_l )  THEN
348          v_max = v_min_l
349       ELSE
350          v_max = v_max_l
351       ENDIF
352       IF ( ABS( w_min_l ) > w_max_l )  THEN
353          w_max = w_min_l
354       ELSE
355          w_max = w_max_l
356       ENDIF
357#endif
358
359#else
360!
361!--    Variable time step:
362!--    Calculate the maximum time step according to the CFL-criterion,
363!--    individually for each velocity component
364       dt_u_l = 999999.9_wp
365       dt_v_l = 999999.9_wp
366       dt_w_l = 999999.9_wp
367       DO  i = nxl, nxr
368          DO  j = nys, nyn
369             DO  k = nzb+1, nzt
370                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
371                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
372                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
373             ENDDO
374          ENDDO
375       ENDDO
376
377#if defined( __parallel )
378       reduce_l(1) = dt_u_l
379       reduce_l(2) = dt_v_l
380       reduce_l(3) = dt_w_l
381       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
382       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
383       dt_u = reduce(1)
384       dt_v = reduce(2)
385       dt_w = reduce(3)
386#else
387       dt_u = dt_u_l
388       dt_v = dt_v_l
389       dt_w = dt_w_l
390#endif
391
392#endif
393
394!
395!--    Compute time step according to the diffusion criterion.
396!--    First calculate minimum grid spacing which only depends on index k
397!--    Note: also at k=nzb+1 a full grid length is being assumed, although
398!--          in the Prandtl-layer friction term only dz/2 is used.
399!--          Experience from the old model seems to justify this.
400       dt_diff_l = 999999.0_wp
401
402       DO  k = nzb+1, nzt
403           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
404       ENDDO
405
406!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
407!$OMP DO
408       !$acc parallel loop collapse(3) present( kh, km )
409       DO  i = nxl, nxr
410          DO  j = nys, nyn
411             DO  k = nzb+1, nzt
412                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
413                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
414             ENDDO
415          ENDDO
416       ENDDO
417       !$acc end parallel
418!$OMP END PARALLEL
419#if defined( __parallel )
420       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
421       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, &
422                           ierr )
423#else
424       dt_diff = dt_diff_l
425#endif
426
427!
428!--    The time step is the minimum of the 3-4 components and the diffusion time
429!--    step minus a reduction (cfl_factor) to be on the safe side.
430!--    The time step must not exceed the maximum allowed value.
431       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w,   &
432                                 dt_precipitation )
433       dt_3d = MIN( dt_3d, dt_max )
434
435!
436!--    Remember the restricting time step criterion for later output.
437       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
438          timestep_reason = 'A'
439       ELSE
440          timestep_reason = 'D'
441       ENDIF
442
443!
444!--    Set flag if the time step becomes too small.
445       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
446          stop_dt = .TRUE.
447
448          WRITE( message_string, * ) 'Time step has reached minimum limit.',  &
449               '&dt              = ', dt_3d, ' s  Simulation is terminated.', &
450               '&old_dt          = ', old_dt, ' s',                           &
451               '&dt_u            = ', dt_u, ' s',                             &
452               '&dt_v            = ', dt_v, ' s',                             &
453               '&dt_w            = ', dt_w, ' s',                             &
454               '&dt_diff         = ', dt_diff, ' s',                          &
455               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
456               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
457               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),       &
458               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
459               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),       &
460               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
461          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
462!
463!--       In case of coupled runs inform the remote model of the termination
464!--       and its reason, provided the remote model has not already been
465!--       informed of another termination reason (terminate_coupled > 0) before.
466#if defined( __parallel )
467          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
468             terminate_coupled = 2
469             IF ( myid == 0 ) THEN
470                CALL MPI_SENDRECV( &
471                     terminate_coupled,        1, MPI_INTEGER, target_id,  0, &
472                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0, &
473                     comm_inter, status, ierr )
474             ENDIF
475             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
476          ENDIF
477#endif
478       ENDIF
479
480!
481!--    Ensure a smooth value (two significant digits) of the timestep.
482       div = 1000.0_wp
483       DO  WHILE ( dt_3d < div )
484          div = div / 10.0_wp
485       ENDDO
486       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
487
488!
489!--    Adjust the time step
490       old_dt = dt_3d
491
492    ENDIF
493
494    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
495
496 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.