Home

Context Navigation

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

Last change on this file since 1433 was 1343, checked in by kanani, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 21.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |