Home

Context Navigation

timestep.f90 @ 4116

Last change on this file since 4116 was 4101, checked in by gronemeier, 5 years ago

timestep.f90:

consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
in RANS mode, instead of considering each wind component individually use the wind speed of 3d wind vector in CFL criterion
do not limit the increase of dt based on its previous value in RANS mode

other:

remove dt_old

Property svn:keywords set to Id

File size: 18.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: