Home

Context Navigation

timestep.f90 @ 4875

Last change on this file since 4875 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 18.8 KB

Rev	Line
[1682]	1	!> @file timestep.f90
[4540]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4540]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4540]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4540]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4540]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4540]	19	!
[258]	20	! Current revisions:
[4540]	21	! -----------------
[4717]	22	!
	23	!
[1485]	24	! Former revisions:
	25	! -----------------
	26	! $Id: timestep.f90 4828 2021-01-05 11:21:41Z scharf $
[4717]	27	! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
	28	! directives (J. Resler)
	29	!
	30	! 4564 2020-06-12 14:03:36Z raasch
[4564]	31	! Vertical nesting method of Huq et al. (2019) removed
	32	!
	33	! 4540 2020-05-18 15:23:29Z raasch
[4540]	34	! File re-formatted to follow the PALM coding standard
	35	!
	36	! 4444 2020-03-05 15:59:50Z raasch
	37	! Bugfix: cpp-directives for serial mode added
	38	!
[4444]	39	! 4360 2020-01-07 11:25:50Z suehring
[4237]	40	! Added missing OpenMP directives
[4540]	41	!
[4237]	42	! 4233 2019-09-20 09:55:54Z knoop
[4233]	43	! OpenACC data update host removed
[4540]	44	!
[4233]	45	! 4182 2019-08-22 15:20:23Z scharf
[4182]	46	! Corrected "Former revisions" section
[4540]	47	!
[4182]	48	! 4101 2019-07-17 15:14:26Z gronemeier
[4540]	49	! - Consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
	50	! - in RANS mode, instead of considering each wind component individually use the wind speed of 3d
	51	! wind vector in CFL criterion
	52	! - Do not limit the increase of dt based on its previous value in RANS mode
[4101]	53	!
	54	! 3658 2019-01-07 20:28:54Z knoop
[3634]	55	! OpenACC port for SPEC
[2716]	56	!
[4182]	57	! Revision 1.1 1997/08/11 06:26:19 raasch
	58	! Initial revision
	59	!
	60	!
[1]	61	! Description:
	62	! ------------
[1682]	63	!> Compute the time step under consideration of the FCL and diffusion criterion.
[4540]	64	!--------------------------------------------------------------------------------------------------!
[1682]	65	SUBROUTINE timestep
[1]	66
[1320]	67
[4540]	68	USE arrays_3d, &
	69	ONLY: dzu, &
	70	dzw, &
	71	kh, &
	72	km, &
	73	u, &
	74	u_stokes_zu, &
	75	v, &
	76	v_stokes_zu, &
	77	w
[1320]	78
[4540]	79	USE control_parameters, &
	80	ONLY: cfl_factor, &
	81	dt_3d, &
	82	dt_fixed, &
	83	dt_max, &
	84	galilei_transformation, &
	85	message_string, &
	86	rans_mode, &
	87	stop_dt, &
	88	timestep_reason, &
	89	u_gtrans, &
	90	use_ug_for_galilei_tr, &
	91	v_gtrans
	92
[4444]	93	#if defined( __parallel )
[4540]	94	USE control_parameters, &
	95	ONLY: coupling_mode, &
	96	terminate_coupled, &
	97	terminate_coupled_remote
[4444]	98	#endif
	99
[4540]	100	USE cpulog, &
	101	ONLY: cpu_log, &
	102	log_point
[1320]	103
[4540]	104	USE grid_variables, &
	105	ONLY: dx, &
	106	dx2, &
	107	dy, &
	108	dy2
[1320]	109
[4540]	110	USE indices, &
	111	ONLY: nxl, &
	112	nxlg, &
	113	nxr, &
	114	nxrg, &
	115	nyn, &
	116	nyng, &
	117	nys, &
	118	nysg, &
	119	nzb, &
	120	nzt
[1320]	121
[1]	122	USE interfaces
[1320]	123
	124	USE kinds
	125
[4540]	126	USE bulk_cloud_model_mod, &
[1849]	127	ONLY: dt_precipitation
	128
[1]	129	USE pegrid
	130
[4540]	131	USE pmc_interface, &
[2130]	132	ONLY: nested_run
	133
[4540]	134	USE statistics, &
	135	ONLY: flow_statistics_called, &
	136	hom, &
	137	u_max, &
	138	u_max_ijk, &
	139	v_max, &
	140	v_max_ijk, &
	141	w_max, &
	142	w_max_ijk
[1320]	143
[1]	144	IMPLICIT NONE
	145
[4540]	146	INTEGER(iwp) :: i !<
	147	INTEGER(iwp) :: j !<
	148	INTEGER(iwp) :: k !<
[3083]	149	INTEGER(iwp) :: km_max_ijk(3) = -1 !< index values (i,j,k) of location where km_max occurs
	150	INTEGER(iwp) :: kh_max_ijk(3) = -1 !< index values (i,j,k) of location where kh_max occurs
[1]	151
[4540]	152	LOGICAL :: stop_dt_local !< local switch for controlling the time stepping
[2130]	153
[4540]	154	REAL(wp) :: div !<
	155	REAL(wp) :: dt_diff !<
	156	REAL(wp) :: dt_diff_l !<
	157	REAL(wp) :: dt_u !<
	158	REAL(wp) :: dt_u_l !<
	159	REAL(wp) :: dt_v !<
	160	REAL(wp) :: dt_v_l !<
	161	REAL(wp) :: dt_w !<
	162	REAL(wp) :: dt_w_l !<
	163	REAL(wp) :: km_max !< maximum of Km in entire domain
	164	REAL(wp) :: kh_max !< maximum of Kh in entire domain
	165	REAL(wp) :: u_gtrans_l !<
	166	REAL(wp) :: v_gtrans_l !<
	167
	168	REAL(wp), DIMENSION(2) :: uv_gtrans_l !<
[4444]	169	#if defined( __parallel )
[4540]	170	REAL(wp), DIMENSION(2) :: uv_gtrans !<
	171	REAL(wp), DIMENSION(3) :: reduce !<
	172	REAL(wp), DIMENSION(3) :: reduce_l !<
[4444]	173	#endif
[4540]	174	REAL(wp), DIMENSION(nzb+1:nzt) :: dxyz2_min !<
[3634]	175	!$ACC DECLARE CREATE(dxyz2_min)
[1]	176
	177
	178	CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
[3658]	179
[3083]	180	!
[4540]	181	!-- In case of Galilei-transform not using the geostrophic wind as translation velocity, compute the
	182	!-- volume-averaged horizontal velocity components, which will then be subtracted from the
	183	!-- horizontal wind for the time step and horizontal advection routines.
[1]	184	IF ( galilei_transformation .AND. .NOT. use_ug_for_galilei_tr ) THEN
	185	IF ( flow_statistics_called ) THEN
	186	!
[4540]	187	!-- Horizontal averages already existent, just need to average them vertically.
[1342]	188	u_gtrans = 0.0_wp
	189	v_gtrans = 0.0_wp
[1]	190	DO k = nzb+1, nzt
	191	u_gtrans = u_gtrans + hom(k,1,1,0)
	192	v_gtrans = v_gtrans + hom(k,1,2,0)
	193	ENDDO
[4540]	194	u_gtrans = u_gtrans / REAL( nzt - nzb, KIND = wp )
	195	v_gtrans = v_gtrans / REAL( nzt - nzb, KIND = wp )
[1]	196	ELSE
	197	!
	198	!-- Averaging over the entire model domain.
[1342]	199	u_gtrans_l = 0.0_wp
	200	v_gtrans_l = 0.0_wp
[1]	201	DO i = nxl, nxr
	202	DO j = nys, nyn
	203	DO k = nzb+1, nzt
[1257]	204	u_gtrans_l = u_gtrans_l + u(k,j,i)
	205	v_gtrans_l = v_gtrans_l + v(k,j,i)
[1]	206	ENDDO
	207	ENDDO
	208	ENDDO
[4540]	209	uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1) * (nyn-nys+1) * (nzt-nzb), KIND = wp )
	210	uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1) * (nyn-nys+1) * (nzt-nzb), KIND = wp )
[1]	211	#if defined( __parallel )
[622]	212	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4540]	213	CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, comm2d, ierr )
	214	u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND = wp )
	215	v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND = wp )
[1]	216	#else
	217	u_gtrans = uv_gtrans_l(1)
	218	v_gtrans = uv_gtrans_l(2)
	219	#endif
	220	ENDIF
	221	ENDIF
	222
[866]	223	!
[4540]	224	!-- Determine the maxima of the velocity components, including their grid index positions.
[4823]	225	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, u_max, u_max_ijk )
	226	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, v_max, v_max_ijk )
[4540]	227	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, w_max, w_max_ijk )
[866]	228
[1257]	229	IF ( .NOT. dt_fixed ) THEN
[866]	230	!
[1257]	231	!-- Variable time step:
[4101]	232	!-- Calculate the maximum time step according to the CFL-criterion
[1342]	233	dt_u_l = 999999.9_wp
	234	dt_v_l = 999999.9_wp
	235	dt_w_l = 999999.9_wp
[4101]	236
	237	IF ( .NOT. rans_mode ) THEN
	238	!
	239	!-- Consider each velocity component individually
	240
	241	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	242	!$ACC COPY(dt_u_l, dt_v_l, dt_w_l, u_stokes_zu, v_stokes_zu) &
	243	!$ACC REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l) &
	244	!$ACC PRESENT(u, v, w, dzu)
[4237]	245	!$OMP PARALLEL DO PRIVATE(i,j,k) &
	246	!$OMP REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l)
[4101]	247	DO i = nxl, nxr
	248	DO j = nys, nyn
	249	DO k = nzb+1, nzt
[4540]	250	dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans + u_stokes_zu(k) ) &
	251	+ 1.0E-10_wp ) ) )
	252	dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans + v_stokes_zu(k) ) &
	253	+ 1.0E-10_wp ) ) )
	254	dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10_wp ) ) )
[4101]	255	ENDDO
[1257]	256	ENDDO
	257	ENDDO
[1]	258
[4101]	259	ELSE
	260	!
	261	!-- Consider the wind speed at the scalar-grid point
[4540]	262	!-- !> @note Considering the wind speed instead of each individual wind component is only a
	263	!-- !> workaround so far. This has to be changed in the future.
[4101]	264
	265	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	266	!$ACC COPY(dt_u_l, u_stokes_zu, v_stokes_zu) &
	267	!$ACC REDUCTION(MIN: dt_u_l) &
	268	!$ACC PRESENT(u, v, w, dzu)
[4237]	269	!$OMP PARALLEL DO PRIVATE(i,j,k) &
	270	!$OMP REDUCTION(MIN: dt_u_l)
[4101]	271	DO i = nxl, nxr
	272	DO j = nys, nyn
	273	DO k = nzb+1, nzt
[4540]	274	dt_u_l = MIN( dt_u_l, ( MIN( dx, dy, dzu(k) ) / ( SQRT( &
	275	( 0.5 * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans + u_stokes_zu(k) )**2 &
	276	+ ( 0.5 * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans + v_stokes_zu(k) )**2 &
	277	+ ( 0.5 * ( w(k,j,i) + w(k-1,j,i) ) )**2 ) + 1.0E-10_wp ) ) )
[4101]	278	ENDDO
	279	ENDDO
	280	ENDDO
[4540]	281
[4101]	282	dt_v_l = dt_u_l
	283	dt_w_l = dt_u_l
	284
	285	ENDIF
	286
[1257]	287	#if defined( __parallel )
	288	reduce_l(1) = dt_u_l
	289	reduce_l(2) = dt_v_l
	290	reduce_l(3) = dt_w_l
	291	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	292	CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
	293	dt_u = reduce(1)
	294	dt_v = reduce(2)
	295	dt_w = reduce(3)
	296	#else
	297	dt_u = dt_u_l
	298	dt_v = dt_v_l
	299	dt_w = dt_w_l
	300	#endif
	301
[1]	302	!
	303	!-- Compute time step according to the diffusion criterion.
[4540]	304	!-- First calculate minimum grid spacing which only depends on index k. When using the dynamic
	305	!-- subgrid model, negative km are possible.
[1342]	306	dt_diff_l = 999999.0_wp
[1]	307
[3634]	308	!$ACC PARALLEL LOOP PRESENT(dxyz2_min, dzw)
[1]	309	DO k = nzb+1, nzt
[4540]	310	dxyz2_min(k) = MIN( dx2, dy2, dzw(k) * dzw(k) ) * 0.125_wp
[1]	311	ENDDO
	312
[4717]	313	!$OMP PARALLEL DO private(i,j,k) reduction(MIN: dt_diff_l)
[3634]	314	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	315	!$ACC COPY(dt_diff_l) REDUCTION(MIN: dt_diff_l) &
	316	!$ACC PRESENT(dxyz2_min, kh, km)
[1]	317	DO i = nxl, nxr
	318	DO j = nys, nyn
	319	DO k = nzb+1, nzt
[4540]	320	dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / ( MAX( kh(k,j,i), 2.0_wp * &
	321	ABS( km(k,j,i) ) ) + 1E-20_wp ) )
[1]	322	ENDDO
	323	ENDDO
	324	ENDDO
	325	#if defined( __parallel )
[622]	326	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4540]	327	CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
[1]	328	#else
	329	dt_diff = dt_diff_l
	330	#endif
	331
	332	!
[4540]	333	!-- The time step is the minimum of the 3-4 components and the diffusion time step minus a
	334	!-- reduction (cfl_factor) to be on the safe side.
[3084]	335	!-- The time step must not exceed the maximum allowed value.
[2130]	336	dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
[3084]	337	dt_3d = MIN( dt_3d, dt_max )
[1]	338
	339	!
	340	!-- Remember the restricting time step criterion for later output.
[1484]	341	IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff ) THEN
[1]	342	timestep_reason = 'A'
	343	ELSE
	344	timestep_reason = 'D'
	345	ENDIF
	346
	347	!
	348	!-- Set flag if the time step becomes too small.
[1342]	349	IF ( dt_3d < ( 0.00001_wp * dt_max ) ) THEN
[1]	350	stop_dt = .TRUE.
[108]	351
[3083]	352	!
[4540]	353	!-- Determine the maxima of the diffusion coefficients, including their grid index positions.
	354	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, km, 'abs', 0.0_wp, km_max, &
	355	km_max_ijk )
	356	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, kh, 'abs', 0.0_wp, kh_max, &
	357	kh_max_ijk )
[3083]	358
[4540]	359	WRITE( message_string, * ) 'Time step has reached minimum limit.', &
	360	'&dt = ', dt_3d, ' s Simulation is terminated.', &
	361	'&dt_u = ', dt_u, ' s', &
	362	'&dt_v = ', dt_v, ' s', &
	363	'&dt_w = ', dt_w, ' s', &
	364	'&dt_diff = ', dt_diff, ' s', &
	365	'&u_max = ', u_max, ' m/s k=', u_max_ijk(1), &
	366	' j=', u_max_ijk(2), ' i=', u_max_ijk(3), &
	367	'&v_max = ', v_max, ' m/s k=', v_max_ijk(1), &
	368	' j=', v_max_ijk(2), ' i=', v_max_ijk(3), &
	369	'&w_max = ', w_max, ' m/s k=', w_max_ijk(1), &
	370	' j=', w_max_ijk(2), ' i=', w_max_ijk(3), &
	371	'&km_max = ', km_max, ' m2/s2 k=', km_max_ijk(1), &
	372	' j=', km_max_ijk(2), ' i=', km_max_ijk(3), &
	373	'&kh_max = ', kh_max, ' m2/s2 k=', kh_max_ijk(1), &
[3083]	374	' j=', kh_max_ijk(2), ' i=', kh_max_ijk(3)
[258]	375	CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
[108]	376	!
[4540]	377	!-- In case of coupled runs inform the remote model of the termination and its reason,
	378	!-- provided the remote model has not already been informed of another termination reason
	379	!-- (terminate_coupled > 0).
[222]	380	#if defined( __parallel )
[108]	381	IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 ) THEN
	382	terminate_coupled = 2
[2130]	383	IF ( myid == 0 ) THEN
[4540]	384	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, target_id, 0, &
	385	terminate_coupled_remote, 1, MPI_INTEGER, target_id, 0, &
	386	comm_inter, status, ierr )
[667]	387	ENDIF
[4540]	388	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
[108]	389	ENDIF
[222]	390	#endif
[1]	391	ENDIF
	392
	393	!
[4540]	394	!-- In case of nested runs all parent/child processes have to terminate if one process has set
	395	!-- the stop flag, i.e. they need to set the stop flag too.
[2130]	396	IF ( nested_run ) THEN
	397	stop_dt_local = stop_dt
[2258]	398	#if defined( __parallel )
[4540]	399	CALL MPI_ALLREDUCE( stop_dt_local, stop_dt, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
[2258]	400	#endif
[2130]	401	ENDIF
	402
	403	!
[1001]	404	!-- Ensure a smooth value (two significant digits) of the timestep.
[1342]	405	div = 1000.0_wp
[1001]	406	DO WHILE ( dt_3d < div )
[1342]	407	div = div / 10.0_wp
[1001]	408	ENDDO
[1342]	409	dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
[1]	410
[1001]	411	ENDIF
[1]	412
	413	CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
	414
	415	END SUBROUTINE timestep

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

Context Navigation

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

Download in other formats: