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

Last change on this file since 1811 was 1683, checked in by knoop, 9 years ago

last commit documented

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