Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |