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

Last change on this file since 1510 was 1485, checked in by kanani, 10 years ago

last commit documented

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