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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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