Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |