Home

Context Navigation

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

Last change on this file since 4556 was 4540, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 19.0 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	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4540]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4540]	19	!
[258]	20	! Current revisions:
[4540]	21	! -----------------
[1485]	22	!
[3049]	23	!
[1485]	24	! Former revisions:
	25	! -----------------
	26	! $Id: timestep.f90 4540 2020-05-18 15:23:29Z raasch $
[4540]	27	! File re-formatted to follow the PALM coding standard
	28	!
	29	!
	30	! 4444 2020-03-05 15:59:50Z raasch
	31	! Bugfix: cpp-directives for serial mode added
	32	!
[4444]	33	! 4360 2020-01-07 11:25:50Z suehring
[4237]	34	! Added missing OpenMP directives
[4540]	35	!
[4237]	36	! 4233 2019-09-20 09:55:54Z knoop
[4233]	37	! OpenACC data update host removed
[4540]	38	!
[4233]	39	! 4182 2019-08-22 15:20:23Z scharf
[4182]	40	! Corrected "Former revisions" section
[4540]	41	!
[4182]	42	! 4101 2019-07-17 15:14:26Z gronemeier
[4540]	43	! - Consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
	44	! - in RANS mode, instead of considering each wind component individually use the wind speed of 3d
	45	! wind vector in CFL criterion
	46	! - Do not limit the increase of dt based on its previous value in RANS mode
[4101]	47	!
	48	! 3658 2019-01-07 20:28:54Z knoop
[3634]	49	! OpenACC port for SPEC
[2716]	50	!
[4182]	51	! Revision 1.1 1997/08/11 06:26:19 raasch
	52	! Initial revision
	53	!
	54	!
[1]	55	! Description:
	56	! ------------
[1682]	57	!> Compute the time step under consideration of the FCL and diffusion criterion.
[4540]	58	!--------------------------------------------------------------------------------------------------!
[1682]	59	SUBROUTINE timestep
[1]	60
[1320]	61
[4540]	62	USE arrays_3d, &
	63	ONLY: dzu, &
	64	dzw, &
	65	kh, &
	66	km, &
	67	u, &
	68	u_stokes_zu, &
	69	v, &
	70	v_stokes_zu, &
	71	w
[1320]	72
[4540]	73	USE control_parameters, &
	74	ONLY: cfl_factor, &
	75	dt_3d, &
	76	dt_fixed, &
	77	dt_max, &
	78	galilei_transformation, &
	79	message_string, &
	80	rans_mode, &
	81	stop_dt, &
	82	timestep_reason, &
	83	u_gtrans, &
	84	use_ug_for_galilei_tr, &
	85	v_gtrans
	86
[4444]	87	#if defined( __parallel )
[4540]	88	USE control_parameters, &
	89	ONLY: coupling_mode, &
	90	terminate_coupled, &
	91	terminate_coupled_remote
[4444]	92	#endif
	93
[4540]	94	USE cpulog, &
	95	ONLY: cpu_log, &
	96	log_point
[1320]	97
[4540]	98	USE grid_variables, &
	99	ONLY: dx, &
	100	dx2, &
	101	dy, &
	102	dy2
[1320]	103
[4540]	104	USE indices, &
	105	ONLY: nxl, &
	106	nxlg, &
	107	nxr, &
	108	nxrg, &
	109	nyn, &
	110	nyng, &
	111	nys, &
	112	nysg, &
	113	nzb, &
	114	nzt
[1320]	115
[1]	116	USE interfaces
[1320]	117
	118	USE kinds
	119
[4540]	120	USE bulk_cloud_model_mod, &
[1849]	121	ONLY: dt_precipitation
	122
[1]	123	USE pegrid
	124
[4540]	125	USE pmc_interface, &
[2130]	126	ONLY: nested_run
	127
[4540]	128	USE statistics, &
	129	ONLY: flow_statistics_called, &
	130	hom, &
	131	u_max, &
	132	u_max_ijk, &
	133	v_max, &
	134	v_max_ijk, &
	135	w_max, &
	136	w_max_ijk
[1320]	137
[4444]	138	#if defined( __parallel )
[4540]	139	USE vertical_nesting_mod, &
	140	ONLY: vnested, &
	141	vnest_timestep_sync
[4444]	142	#endif
[2365]	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.
	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 )
	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
[3241]	313	!$OMP PARALLEL private(i,j,k) reduction(MIN: dt_diff_l)
[2118]	314	!$OMP DO
[3634]	315	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	316	!$ACC COPY(dt_diff_l) REDUCTION(MIN: dt_diff_l) &
	317	!$ACC PRESENT(dxyz2_min, kh, km)
[1]	318	DO i = nxl, nxr
	319	DO j = nys, nyn
	320	DO k = nzb+1, nzt
[4540]	321	dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / ( MAX( kh(k,j,i), 2.0_wp * &
	322	ABS( km(k,j,i) ) ) + 1E-20_wp ) )
[1]	323	ENDDO
	324	ENDDO
	325	ENDDO
[2118]	326	!$OMP END PARALLEL
[1]	327	#if defined( __parallel )
[622]	328	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4540]	329	CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
[1]	330	#else
	331	dt_diff = dt_diff_l
	332	#endif
	333
	334	!
[4540]	335	!-- The time step is the minimum of the 3-4 components and the diffusion time step minus a
	336	!-- reduction (cfl_factor) to be on the safe side.
[3084]	337	!-- The time step must not exceed the maximum allowed value.
[2130]	338	dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
[3084]	339	dt_3d = MIN( dt_3d, dt_max )
[1]	340
	341	!
	342	!-- Remember the restricting time step criterion for later output.
[1484]	343	IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff ) THEN
[1]	344	timestep_reason = 'A'
	345	ELSE
	346	timestep_reason = 'D'
	347	ENDIF
	348
	349	!
	350	!-- Set flag if the time step becomes too small.
[1342]	351	IF ( dt_3d < ( 0.00001_wp * dt_max ) ) THEN
[1]	352	stop_dt = .TRUE.
[108]	353
[3083]	354	!
[4540]	355	!-- Determine the maxima of the diffusion coefficients, including their grid index positions.
	356	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, km, 'abs', 0.0_wp, km_max, &
	357	km_max_ijk )
	358	CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, kh, 'abs', 0.0_wp, kh_max, &
	359	kh_max_ijk )
[3083]	360
[4540]	361	WRITE( message_string, * ) 'Time step has reached minimum limit.', &
	362	'&dt = ', dt_3d, ' s Simulation is terminated.', &
	363	'&dt_u = ', dt_u, ' s', &
	364	'&dt_v = ', dt_v, ' s', &
	365	'&dt_w = ', dt_w, ' s', &
	366	'&dt_diff = ', dt_diff, ' s', &
	367	'&u_max = ', u_max, ' m/s k=', u_max_ijk(1), &
	368	' j=', u_max_ijk(2), ' i=', u_max_ijk(3), &
	369	'&v_max = ', v_max, ' m/s k=', v_max_ijk(1), &
	370	' j=', v_max_ijk(2), ' i=', v_max_ijk(3), &
	371	'&w_max = ', w_max, ' m/s k=', w_max_ijk(1), &
	372	' j=', w_max_ijk(2), ' i=', w_max_ijk(3), &
	373	'&km_max = ', km_max, ' m2/s2 k=', km_max_ijk(1), &
	374	' j=', km_max_ijk(2), ' i=', km_max_ijk(3), &
	375	'&kh_max = ', kh_max, ' m2/s2 k=', kh_max_ijk(1), &
[3083]	376	' j=', kh_max_ijk(2), ' i=', kh_max_ijk(3)
[258]	377	CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
[108]	378	!
[4540]	379	!-- In case of coupled runs inform the remote model of the termination and its reason,
	380	!-- provided the remote model has not already been informed of another termination reason
	381	!-- (terminate_coupled > 0).
[222]	382	#if defined( __parallel )
[108]	383	IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 ) THEN
	384	terminate_coupled = 2
[2130]	385	IF ( myid == 0 ) THEN
[4540]	386	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, target_id, 0, &
	387	terminate_coupled_remote, 1, MPI_INTEGER, target_id, 0, &
	388	comm_inter, status, ierr )
[667]	389	ENDIF
[4540]	390	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
[108]	391	ENDIF
[222]	392	#endif
[1]	393	ENDIF
	394
	395	!
[4540]	396	!-- In case of nested runs all parent/child processes have to terminate if one process has set
	397	!-- the stop flag, i.e. they need to set the stop flag too.
[2130]	398	IF ( nested_run ) THEN
	399	stop_dt_local = stop_dt
[2258]	400	#if defined( __parallel )
[4540]	401	CALL MPI_ALLREDUCE( stop_dt_local, stop_dt, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
[2258]	402	#endif
[2130]	403	ENDIF
	404
	405	!
[1001]	406	!-- Ensure a smooth value (two significant digits) of the timestep.
[1342]	407	div = 1000.0_wp
[1001]	408	DO WHILE ( dt_3d < div )
[1342]	409	div = div / 10.0_wp
[1001]	410	ENDDO
[1342]	411	dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
[1]	412
[1001]	413	ENDIF
[1]	414
[4444]	415	#if defined( __parallel )
[2365]	416	!
[4540]	417	!-- Vertical nesting: coarse and fine grid timestep has to be identical
[2365]	418	IF ( vnested ) CALL vnest_timestep_sync
[4444]	419	#endif
[2365]	420
[1]	421	CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
	422
	423	END SUBROUTINE timestep

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |