Home

Context Navigation

source: palm/trunk/SOURCE/init_1d_model.f90 @ 243

Last change on this file since 243 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 34.0 KB

Rev	Line
[1]	1	SUBROUTINE init_1d_model
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[198]	6	!
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: init_1d_model.f90 198 2008-09-17 08:55:28Z letzel $
[77]	11	!
[198]	12	! 184 2008-08-04 15:53:39Z letzel
	13	! provisional solution for run_control_1d output: add 'CALL check_open( 15 )'
	14	!
[139]	15	! 135 2007-11-22 12:24:23Z raasch
	16	! Bugfix: absolute value of f must be used when calculating the Blackadar
	17	! mixing length
	18	!
[83]	19	! 82 2007-04-16 15:40:52Z raasch
	20	! Preprocessor strings for different linux clusters changed to "lc",
	21	! routine local_flush is used for buffer flushing
	22	!
[77]	23	! 75 2007-03-22 09:54:05Z raasch
	24	! Bugfix: preset of tendencies te_em, te_um, te_vm,
	25	! moisture renamed humidity
	26	!
[3]	27	! RCS Log replace by Id keyword, revision history cleaned up
	28	!
[1]	29	! Revision 1.21 2006/06/02 15:19:57 raasch
	30	! cpp-directives extended for lctit
	31	!
	32	! Revision 1.1 1998/03/09 16:22:10 raasch
	33	! Initial revision
	34	!
	35	!
	36	! Description:
	37	! ------------
	38	! 1D-model to initialize the 3D-arrays.
	39	! The temperature profile is set as steady and a corresponding steady solution
	40	! of the wind profile is being computed.
	41	! All subroutines required can be found within this file.
	42	!------------------------------------------------------------------------------!
	43
	44	USE arrays_3d
	45	USE indices
	46	USE model_1d
	47	USE control_parameters
	48
	49	IMPLICIT NONE
	50
	51	CHARACTER (LEN=9) :: time_to_string
	52	INTEGER :: k
	53	REAL :: lambda
	54
	55	!
	56	!-- Allocate required 1D-arrays
	57	ALLOCATE( e1d(nzb:nzt+1), e1d_m(nzb:nzt+1), e1d_p(nzb:nzt+1), &
	58	kh1d(nzb:nzt+1), kh1d_m(nzb:nzt+1), km1d(nzb:nzt+1), &
	59	km1d_m(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1), &
	60	l1d_m(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1), &
	61	te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &
	62	te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &
	63	u1d_m(nzb:nzt+1), u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), &
	64	v1d_m(nzb:nzt+1), v1d_p(nzb:nzt+1) )
	65
	66	!
	67	!-- Initialize arrays
	68	IF ( constant_diffusion ) THEN
	69	km1d = km_constant
	70	km1d_m = km_constant
	71	kh1d = km_constant / prandtl_number
	72	kh1d_m = km_constant / prandtl_number
	73	ELSE
	74	e1d = 0.0; e1d_m = 0.0; e1d_p = 0.0
	75	kh1d = 0.0; kh1d_m = 0.0; km1d = 0.0; km1d_m = 0.0
	76	rif1d = 0.0
	77	!
	78	!-- Compute the mixing length
	79	l_black(nzb) = 0.0
	80
	81	IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN
	82	!
	83	!-- Blackadar mixing length
	84	IF ( f /= 0.0 ) THEN
[135]	85	lambda = 2.7E-4 * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
	86	ABS( f ) + 1E-10
[1]	87	ELSE
	88	lambda = 30.0
	89	ENDIF
	90
	91	DO k = nzb+1, nzt+1
	92	l_black(k) = kappa * zu(k) / ( 1.0 + kappa * zu(k) / lambda )
	93	ENDDO
	94
	95	ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN
	96	!
	97	!-- Use the same mixing length as in 3D model
	98	l_black(1:nzt) = l_grid
	99	l_black(nzt+1) = l_black(nzt)
	100
	101	ENDIF
	102
	103	!
	104	!-- Adjust mixing length to the prandtl mixing length (within the prandtl
	105	!-- layer)
	106	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	107	k = nzb+1
	108	l_black(k) = MIN( l_black(k), kappa * zu(k) )
	109	ENDIF
	110	ENDIF
	111	l1d = l_black
	112	l1d_m = l_black
	113	u1d = u_init
	114	u1d_m = u_init
	115	u1d_p = u_init
	116	v1d = v_init
	117	v1d_m = v_init
	118	v1d_p = v_init
	119
	120	!
	121	!-- Set initial horizontal velocities at the lowest grid levels to a very small
	122	!-- value in order to avoid too small time steps caused by the diffusion limit
	123	!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
	124	u1d(0:1) = 0.1
	125	u1d_m(0:1) = 0.1
	126	u1d_p(0:1) = 0.1
	127	v1d(0:1) = 0.1
	128	v1d_m(0:1) = 0.1
	129	v1d_p(0:1) = 0.1
	130
	131	!
	132	!-- For u, theta and the momentum fluxes plausible values are set
	133	IF ( prandtl_layer ) THEN
	134	us1d = 0.1 ! without initial friction the flow would not change
	135	ELSE
	136	e1d(nzb+1) = 1.0
	137	km1d(nzb+1) = 1.0
	138	us1d = 0.0
	139	ENDIF
	140	ts1d = 0.0
	141	usws1d = 0.0; usws1d_m = 0.0
	142	vsws1d = 0.0; vsws1d_m = 0.0
	143	z01d = roughness_length
[75]	144	IF ( humidity .OR. passive_scalar ) qs1d = 0.0
[1]	145
	146	!
[46]	147	!-- Tendencies must be preset in order to avoid runtime errors within the
	148	!-- first Runge-Kutta step
	149	te_em = 0.0
	150	te_um = 0.0
	151	te_vm = 0.0
	152
	153	!
[1]	154	!-- Set start time in hh:mm:ss - format
	155	simulated_time_chr = time_to_string( simulated_time_1d )
	156
	157	!
	158	!-- Integrate the 1D-model equations using the leap-frog scheme
	159	CALL time_integration_1d
	160
	161
	162	END SUBROUTINE init_1d_model
	163
	164
	165
	166	SUBROUTINE time_integration_1d
	167
	168	!------------------------------------------------------------------------------!
	169	! Description:
	170	! ------------
	171	! Leap-frog time differencing scheme for the 1D-model.
	172	!------------------------------------------------------------------------------!
	173
	174	USE arrays_3d
	175	USE control_parameters
	176	USE indices
	177	USE model_1d
	178	USE pegrid
	179
	180	IMPLICIT NONE
	181
	182	CHARACTER (LEN=9) :: time_to_string
	183	INTEGER :: k
	184	REAL :: a, b, dissipation, dpt_dz, flux, kmzm, kmzp, l_stable, pt_0, &
	185	uv_total
	186
	187	!
	188	!-- Determine the time step at the start of a 1D-simulation and
	189	!-- determine and printout quantities used for run control
	190	CALL timestep_1d
	191	CALL run_control_1d
	192
	193	!
	194	!-- Start of time loop
	195	DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d )
	196
	197	!
	198	!-- Depending on the timestep scheme, carry out one or more intermediate
	199	!-- timesteps
	200
	201	intermediate_timestep_count = 0
	202	DO WHILE ( intermediate_timestep_count < &
	203	intermediate_timestep_count_max )
	204
	205	intermediate_timestep_count = intermediate_timestep_count + 1
	206
	207	CALL timestep_scheme_steering
	208
	209	!
	210	!-- Compute all tendency terms. If a Prandtl-layer is simulated, k starts
	211	!-- at nzb+2.
	212	DO k = nzb_diff, nzt
	213
	214	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
	215	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
	216	!
	217	!-- u-component
	218	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
	219	kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) &
	220	- kmzm * ( u1d_m(k) - u1d_m(k-1) ) * ddzu(k) &
	221	) * ddzw(k)
	222	!
	223	!-- v-component
	224	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	225	kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) &
	226	- kmzm * ( v1d_m(k) - v1d_m(k-1) ) * ddzu(k) &
	227	) * ddzw(k)
	228	ENDDO
	229	IF ( .NOT. constant_diffusion ) THEN
	230	DO k = nzb_diff, nzt
	231	!
	232	!-- TKE
	233	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
	234	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
[75]	235	IF ( .NOT. humidity ) THEN
[1]	236	pt_0 = pt_init(k)
	237	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	238	ELSE
	239	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
	240	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
	241	0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) &
	242	) * dd2zu(k)
	243	ENDIF
	244
	245	IF ( dissipation_1d == 'detering' ) THEN
	246	!
	247	!-- According to Detering, c_e=0.064
	248	dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
	249	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
	250	dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) &
	251	* e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
	252	ENDIF
	253
	254	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
	255	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
	256	) &
	257	- g / pt_0 * kh1d(k) * flux &
	258	+ ( &
	259	kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &
	260	- kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &
	261	) * ddzw(k) &
	262	- dissipation
	263	ENDDO
	264	ENDIF
	265
	266	!
	267	!-- Tendency terms at the top of the Prandtl-layer.
	268	!-- Finite differences of the momentum fluxes are computed using half the
	269	!-- normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
	270	IF ( prandtl_layer ) THEN
	271
	272	k = nzb+1
	273	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
	274	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
[75]	275	IF ( .NOT. humidity ) THEN
[1]	276	pt_0 = pt_init(k)
	277	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	278	ELSE
	279	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
	280	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
	281	0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) &
	282	) * dd2zu(k)
	283	ENDIF
	284
	285	IF ( dissipation_1d == 'detering' ) THEN
	286	!
	287	!-- According to Detering, c_e=0.064
	288	dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
	289	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
	290	dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) &
	291	* e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
	292	ENDIF
	293
	294	!
	295	!-- u-component
	296	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
	297	kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) + usws1d_m &
	298	) * 2.0 * ddzw(k)
	299	!
	300	!-- v-component
	301	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	302	kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) + vsws1d_m &
	303	) * 2.0 * ddzw(k)
	304	!
	305	!-- TKE
	306	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	307	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	308	) &
	309	- g / pt_0 * kh1d(k) * flux &
	310	+ ( &
	311	kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &
	312	- kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &
	313	) * ddzw(k) &
	314	- dissipation
	315	ENDIF
	316
	317	!
	318	!-- Prognostic equations for all 1D variables
	319	DO k = nzb+1, nzt
	320
	321	u1d_p(k) = ( 1. - tsc(1) ) * u1d_m(k) + &
	322	tsc(1) * u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
	323	tsc(3) * te_um(k) )
	324	v1d_p(k) = ( 1. - tsc(1) ) * v1d_m(k) + &
	325	tsc(1) * v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
	326	tsc(3) * te_vm(k) )
	327
	328	ENDDO
	329	IF ( .NOT. constant_diffusion ) THEN
	330	DO k = nzb+1, nzt
	331
	332	e1d_p(k) = ( 1. - tsc(1) ) * e1d_m(k) + &
	333	tsc(1) * e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
	334	tsc(3) * te_em(k) )
	335
	336	ENDDO
	337	!
	338	!-- Eliminate negative TKE values, which can result from the
	339	!-- integration due to numerical inaccuracies. In such cases the TKE
	340	!-- value is reduced to 10 percent of its old value.
	341	WHERE ( e1d_p < 0.0 ) e1d_p = 0.1 * e1d
	342	ENDIF
	343
	344	!
	345	!-- Calculate tendencies for the next Runge-Kutta step
	346	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	347	IF ( intermediate_timestep_count == 1 ) THEN
	348
	349	DO k = nzb+1, nzt
	350	te_um(k) = te_u(k)
	351	te_vm(k) = te_v(k)
	352	ENDDO
	353
	354	IF ( .NOT. constant_diffusion ) THEN
	355	DO k = nzb+1, nzt
	356	te_em(k) = te_e(k)
	357	ENDDO
	358	ENDIF
	359
	360	ELSEIF ( intermediate_timestep_count < &
	361	intermediate_timestep_count_max ) THEN
	362
	363	DO k = nzb+1, nzt
	364	te_um(k) = -9.5625 * te_u(k) + 5.3125 * te_um(k)
	365	te_vm(k) = -9.5625 * te_v(k) + 5.3125 * te_vm(k)
	366	ENDDO
	367
	368	IF ( .NOT. constant_diffusion ) THEN
	369	DO k = nzb+1, nzt
	370	te_em(k) = -9.5625 * te_e(k) + 5.3125 * te_em(k)
	371	ENDDO
	372	ENDIF
	373
	374	ENDIF
	375	ENDIF
	376
	377
	378	!
	379	!-- Boundary conditions for the prognostic variables.
	380	!-- At the top boundary (nzt+1) u,v and e keep their initial values
	381	!-- (ug(nzt+1), vg(nzt+1), 0), at the bottom boundary the mirror
	382	!-- boundary condition applies to u and v.
	383	!-- The boundary condition for e is set further below ( (u/cm)*2 ).
	384	u1d_p(nzb) = -u1d_p(nzb+1)
	385	v1d_p(nzb) = -v1d_p(nzb+1)
	386
	387	!
	388	!-- If necessary, apply the time filter
	389	IF ( asselin_filter_factor /= 0.0 .AND. &
	390	timestep_scheme(1:5) /= 'runge' ) THEN
	391
	392	u1d = u1d + asselin_filter_factor * ( u1d_p - 2.0 * u1d + u1d_m )
	393	v1d = v1d + asselin_filter_factor * ( v1d_p - 2.0 * v1d + v1d_m )
	394
	395	IF ( .NOT. constant_diffusion ) THEN
	396	e1d = e1d + asselin_filter_factor * &
	397	( e1d_p - 2.0 * e1d + e1d_m )
	398	ENDIF
	399
	400	ENDIF
	401
	402	!
	403	!-- Swap the time levels in preparation for the next time step.
	404	IF ( timestep_scheme(1:4) == 'leap' ) THEN
	405	u1d_m = u1d
	406	v1d_m = v1d
	407	IF ( .NOT. constant_diffusion ) THEN
	408	e1d_m = e1d
	409	kh1d_m = kh1d ! The old diffusion quantities are required for
	410	km1d_m = km1d ! explicit diffusion in the leap-frog scheme.
	411	l1d_m = l1d
	412	IF ( prandtl_layer ) THEN
	413	usws1d_m = usws1d
	414	vsws1d_m = vsws1d
	415	ENDIF
	416	ENDIF
	417	ENDIF
	418	u1d = u1d_p
	419	v1d = v1d_p
	420	IF ( .NOT. constant_diffusion ) THEN
	421	e1d = e1d_p
	422	ENDIF
	423
	424	!
	425	!-- Compute diffusion quantities
	426	IF ( .NOT. constant_diffusion ) THEN
	427
	428	!
	429	!-- First compute the vertical fluxes in the Prandtl-layer
	430	IF ( prandtl_layer ) THEN
	431	!
	432	!-- Compute theta* using Rif numbers of the previous time step
	433	IF ( rif1d(1) >= 0.0 ) THEN
	434	!
	435	!-- Stable stratification
	436	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	437	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
	438	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	439	)
	440	ELSE
	441	!
	442	!-- Unstable stratification
	443	a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
	444	b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
	445	!
	446	!-- In the borderline case the formula for stable stratification
	447	!-- must be applied, because otherwise a zero division would
	448	!-- occur in the argument of the logarithm.
	449	IF ( a == 0.0 .OR. b == 0.0 ) THEN
	450	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	451	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
	452	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	453	)
	454	ELSE
	455	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	456	LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) )
	457	ENDIF
	458	ENDIF
	459
	460	ENDIF ! prandtl_layer
	461
	462	!
	463	!-- Compute the Richardson-flux numbers,
	464	!-- first at the top of the Prandtl-layer using u* of the previous
	465	!-- time step (+1E-30, if u* = 0), then in the remaining area. There
	466	!-- the rif-numbers of the previous time step are used.
	467
	468	IF ( prandtl_layer ) THEN
[75]	469	IF ( .NOT. humidity ) THEN
[1]	470	pt_0 = pt_init(nzb+1)
	471	flux = ts1d
	472	ELSE
	473	pt_0 = pt_init(nzb+1) * ( 1.0 + 0.61 * q_init(nzb+1) )
	474	flux = ts1d + 0.61 * pt_init(k) * qs1d
	475	ENDIF
	476	rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
	477	( pt_0 * ( us1d**2 + 1E-30 ) )
	478	ENDIF
	479
	480	DO k = nzb_diff, nzt
[75]	481	IF ( .NOT. humidity ) THEN
[1]	482	pt_0 = pt_init(k)
	483	flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
	484	ELSE
	485	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
	486	flux = ( ( pt_init(k+1) - pt_init(k-1) ) &
	487	+ 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )&
	488	) * dd2zu(k)
	489	ENDIF
	490	IF ( rif1d(k) >= 0.0 ) THEN
	491	rif1d(k) = g / pt_0 * flux / &
	492	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	493	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	494	+ 1E-30 &
	495	)
	496	ELSE
	497	rif1d(k) = g / pt_0 * flux / &
	498	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	499	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	500	+ 1E-30 &
	501	) * ( 1.0 - 16.0 * rif1d(k) )**0.25
	502	ENDIF
	503	ENDDO
	504	!
	505	!-- Richardson-numbers must remain restricted to a realistic value
	506	!-- range. It is exceeded excessively for very small velocities
	507	!-- (u,v --> 0).
	508	WHERE ( rif1d < rif_min ) rif1d = rif_min
	509	WHERE ( rif1d > rif_max ) rif1d = rif_max
	510
	511	!
	512	!-- Compute u* from the absolute velocity value
	513	IF ( prandtl_layer ) THEN
	514	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
	515
	516	IF ( rif1d(nzb+1) >= 0.0 ) THEN
	517	!
	518	!-- Stable stratification
	519	us1d = kappa * uv_total / ( &
	520	LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
	521	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	522	)
	523	ELSE
	524	!
	525	!-- Unstable stratification
	526	a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) )
	527	b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) &
	528	* z01d ) )
	529	!
	530	!-- In the borderline case the formula for stable stratification
	531	!-- must be applied, because otherwise a zero division would
	532	!-- occur in the argument of the logarithm.
	533	IF ( a == 1.0 .OR. b == 1.0 ) THEN
	534	us1d = kappa * uv_total / ( &
	535	LOG( zu(nzb+1) / z01d ) + &
	536	5.0 * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) / &
	537	zu(nzb+1) )
	538	ELSE
	539	us1d = kappa * uv_total / ( &
	540	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) +&
	541	2.0 * ( ATAN( b ) - ATAN( a ) ) &
	542	)
	543	ENDIF
	544	ENDIF
	545
	546	!
	547	!-- Compute the momentum fluxes for the diffusion terms
	548	usws1d = - u1d(nzb+1) / uv_total * us1d**2
	549	vsws1d = - v1d(nzb+1) / uv_total * us1d**2
	550
	551	!
	552	!-- Boundary condition for the turbulent kinetic energy at the top
	553	!-- of the Prandtl-layer. c_m = 0.4 according to Detering.
	554	!-- Additional Neumann condition de/dz = 0 at nzb is set to ensure
	555	!-- compatibility with the 3D model.
	556	IF ( ibc_e_b == 2 ) THEN
	557	e1d(nzb+1) = ( us1d / 0.1 )**2
	558	! e1d(nzb+1) = ( us1d / 0.4 )**2 !not used so far, see also
	559	!prandtl_fluxes
	560	ENDIF
	561	e1d(nzb) = e1d(nzb+1)
	562
[75]	563	IF ( humidity .OR. passive_scalar ) THEN
[1]	564	!
	565	!-- Compute q*
	566	IF ( rif1d(1) >= 0.0 ) THEN
	567	!
	568	!-- Stable stratification
	569	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
	570	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
	571	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	572	)
	573	ELSE
	574	!
	575	!-- Unstable stratification
	576	a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
	577	b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
	578	!
	579	!-- In the borderline case the formula for stable stratification
	580	!-- must be applied, because otherwise a zero division would
	581	!-- occur in the argument of the logarithm.
	582	IF ( a == 1.0 .OR. b == 1.0 ) THEN
	583	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
	584	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
	585	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	586	)
	587	ELSE
	588	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
	589	LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) )
	590	ENDIF
	591	ENDIF
	592	ELSE
	593	qs1d = 0.0
	594	ENDIF
	595
	596	ENDIF ! prandtl_layer
	597
	598	!
	599	!-- Compute the diabatic mixing length
	600	IF ( mixing_length_1d == 'blackadar' ) THEN
	601	DO k = nzb+1, nzt
	602	IF ( rif1d(k) >= 0.0 ) THEN
	603	l1d(k) = l_black(k) / ( 1.0 + 5.0 * rif1d(k) )
	604	ELSE
	605	l1d(k) = l_black(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25
	606	ENDIF
	607	l1d(k) = l_black(k)
	608	ENDDO
	609
	610	ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN
	611	DO k = nzb+1, nzt
	612	dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
	613	IF ( dpt_dz > 0.0 ) THEN
	614	l_stable = 0.76 * SQRT( e1d(k) ) / &
	615	SQRT( g / pt_init(k) * dpt_dz ) + 1E-5
	616	ELSE
	617	l_stable = l_grid(k)
	618	ENDIF
	619	l1d(k) = MIN( l_grid(k), l_stable )
	620	ENDDO
	621	ENDIF
	622
	623	!
	624	!-- Adjust mixing length to the prandtl mixing length
	625	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	626	k = nzb+1
	627	IF ( rif1d(k) >= 0.0 ) THEN
	628	l1d(k) = MIN( l1d(k), kappa * zu(k) / ( 1.0 + 5.0 * &
	629	rif1d(k) ) )
	630	ELSE
	631	l1d(k) = MIN( l1d(k), kappa * zu(k) * &
	632	SQRT( SQRT( 1.0 - 16.0 * rif1d(k) ) ) )
	633	ENDIF
	634	ENDIF
	635
	636	!
	637	!-- Compute the diffusion coefficients for momentum via the
	638	!-- corresponding Prandtl-layer relationship and according to
	639	!-- Prandtl-Kolmogorov, respectively. The unstable stratification is
	640	!-- computed via the adiabatic mixing length, for the unstability has
	641	!-- already been taken account of via the TKE (cf. also Diss.).
	642	IF ( prandtl_layer ) THEN
	643	IF ( rif1d(nzb+1) >= 0.0 ) THEN
	644	km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &
	645	( 1.0 + 5.0 * rif1d(nzb+1) )
	646	ELSE
	647	km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &
	648	( 1.0 - 16.0 * rif1d(nzb+1) )**0.25
	649	ENDIF
	650	ENDIF
	651	DO k = nzb_diff, nzt
	652	! km1d(k) = 0.4 * SQRT( e1d(k) ) !changed: adjustment to 3D-model
	653	km1d(k) = 0.1 * SQRT( e1d(k) )
	654	IF ( rif1d(k) >= 0.0 ) THEN
	655	km1d(k) = km1d(k) * l1d(k)
	656	ELSE
	657	km1d(k) = km1d(k) * l_black(k)
	658	ENDIF
	659	ENDDO
	660
	661	!
	662	!-- Add damping layer
	663	DO k = damp_level_ind_1d+1, nzt+1
	664	km1d(k) = 1.1 * km1d(k-1)
	665	km1d(k) = MIN( km1d(k), 10.0 )
	666	ENDDO
	667
	668	!
	669	!-- Compute the diffusion coefficient for heat via the relationship
	670	!-- kh = phim / phih * km
	671	DO k = nzb+1, nzt
	672	IF ( rif1d(k) >= 0.0 ) THEN
	673	kh1d(k) = km1d(k)
	674	ELSE
	675	kh1d(k) = km1d(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25
	676	ENDIF
	677	ENDDO
	678
	679	ENDIF ! .NOT. constant_diffusion
	680
	681	!
	682	!-- The Runge-Kutta scheme needs the recent diffusion quantities
	683	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	684	u1d_m = u1d
	685	v1d_m = v1d
	686	IF ( .NOT. constant_diffusion ) THEN
	687	e1d_m = e1d
	688	kh1d_m = kh1d
	689	km1d_m = km1d
	690	l1d_m = l1d
	691	IF ( prandtl_layer ) THEN
	692	usws1d_m = usws1d
	693	vsws1d_m = vsws1d
	694	ENDIF
	695	ENDIF
	696	ENDIF
	697
	698
	699	ENDDO ! intermediate step loop
	700
	701	!
	702	!-- Increment simulated time and output times
	703	current_timestep_number_1d = current_timestep_number_1d + 1
	704	simulated_time_1d = simulated_time_1d + dt_1d
	705	simulated_time_chr = time_to_string( simulated_time_1d )
	706	time_pr_1d = time_pr_1d + dt_1d
	707	time_run_control_1d = time_run_control_1d + dt_1d
	708
	709	!
	710	!-- Determine and print out quantities for run control
	711	IF ( time_run_control_1d >= dt_run_control_1d ) THEN
	712	CALL run_control_1d
	713	time_run_control_1d = time_run_control_1d - dt_run_control_1d
	714	ENDIF
	715
	716	!
	717	!-- Profile output on file
	718	IF ( time_pr_1d >= dt_pr_1d ) THEN
	719	CALL print_1d_model
	720	time_pr_1d = time_pr_1d - dt_pr_1d
	721	ENDIF
	722
	723	!
	724	!-- Determine size of next time step
	725	CALL timestep_1d
	726
	727	ENDDO ! time loop
	728
	729
	730	END SUBROUTINE time_integration_1d
	731
	732
	733	SUBROUTINE run_control_1d
	734
	735	!------------------------------------------------------------------------------!
	736	! Description:
	737	! ------------
	738	! Compute and print out quantities for run control of the 1D model.
	739	!------------------------------------------------------------------------------!
	740
	741	USE constants
	742	USE indices
	743	USE model_1d
	744	USE pegrid
	745	USE control_parameters
	746
	747	IMPLICIT NONE
	748
	749	INTEGER :: k
	750	REAL :: alpha, energy, umax, uv_total, vmax
	751
	752	!
	753	!-- Output
	754	IF ( myid == 0 ) THEN
	755	!
	756	!-- If necessary, write header
	757	IF ( .NOT. run_control_header_1d ) THEN
[184]	758	CALL check_open( 15 )
[1]	759	WRITE ( 15, 100 )
	760	run_control_header_1d = .TRUE.
	761	ENDIF
	762
	763	!
	764	!-- Compute control quantities
	765	!-- grid level nzp is excluded due to mirror boundary condition
	766	umax = 0.0; vmax = 0.0; energy = 0.0
	767	DO k = nzb+1, nzt+1
	768	umax = MAX( ABS( umax ), ABS( u1d(k) ) )
	769	vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
	770	energy = energy + 0.5 * ( u1d(k)2 + v1d(k)2 )
	771	ENDDO
	772	energy = energy / REAL( nzt - nzb + 1 )
	773
	774	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
	775	IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5 ) THEN
	776	alpha = ACOS( SIGN( 1.0 , u1d(nzb+1) ) )
	777	ELSE
	778	alpha = ACOS( u1d(nzb+1) / uv_total )
	779	IF ( v1d(nzb+1) <= 0.0 ) alpha = 2.0 * pi - alpha
	780	ENDIF
	781	alpha = alpha / ( 2.0 * pi ) * 360.0
	782
	783	WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &
	784	dt_1d, umax, vmax, us1d, alpha, energy
	785	!
	786	!-- Write buffer contents to disc immediately
[82]	787	CALL local_flush( 15 )
[1]	788
	789	ENDIF
	790
	791	!
	792	!-- formats
	793	100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ &
	794	&'------------------------------'// &
	795	&'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ &
	796	&'-------------------------------------------------------------')
	797	101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,2X,F5.3,2X,F5.1,2X,F7.2)
	798
	799
	800	END SUBROUTINE run_control_1d
	801
	802
	803
	804	SUBROUTINE timestep_1d
	805
	806	!------------------------------------------------------------------------------!
	807	! Description:
	808	! ------------
	809	! Compute the time step w.r.t. the diffusion criterion
	810	!------------------------------------------------------------------------------!
	811
	812	USE arrays_3d
	813	USE indices
	814	USE model_1d
	815	USE pegrid
	816	USE control_parameters
	817
	818	IMPLICIT NONE
	819
	820	INTEGER :: k
	821	REAL :: div, dt_diff, fac, percent_change, value
	822
	823
	824	!
	825	!-- Compute the currently feasible time step according to the diffusion
	826	!-- criterion. At nzb+1 the half grid length is used.
	827	IF ( timestep_scheme(1:4) == 'leap' ) THEN
	828	fac = 0.25
	829	ELSE
	830	fac = 0.35
	831	ENDIF
	832	dt_diff = dt_max_1d
	833	DO k = nzb+2, nzt
	834	value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20 )
	835	dt_diff = MIN( value, dt_diff )
	836	ENDDO
	837	value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20 )
	838	dt_1d = MIN( value, dt_diff )
	839
	840	!
	841	!-- Set flag when the time step becomes too small
	842	IF ( dt_1d < ( 0.00001 * dt_max_1d ) ) THEN
	843	stop_dt_1d = .TRUE.
	844	IF ( myid == 0 ) THEN
	845	PRINT*,'+++ timestep_1d: timestep has exceeded the lower limit'
	846	PRINT*,' dt_1d = ',dt_1d,' s simulation stopped!'
	847	ENDIF
	848	CALL local_stop
	849	ENDIF
	850
	851	IF ( timestep_scheme(1:4) == 'leap' ) THEN
	852
	853	!
	854	!-- The current time step will only be changed if the new time step exceeds
	855	!-- its previous value by 5 % or falls 2 % below. After a time step
	856	!-- reduction at least 30 iterations must be done with this value before a
	857	!-- new reduction will be allowed again.
	858	!-- The control parameters for application of Euler- or leap-frog schemes are
	859	!-- set accordingly.
	860	percent_change = dt_1d / old_dt_1d - 1.0
	861	IF ( percent_change > 0.05 .OR. percent_change < -0.02 ) THEN
	862
	863	!
	864	!-- Each time step increase is by at most 2 %
	865	IF ( percent_change > 0.0 .AND. simulated_time_1d /= 0.0 ) THEN
	866	dt_1d = 1.02 * old_dt_1d
	867	ENDIF
	868
	869	!
	870	!-- A more or less simple new time step value is obtained taking only the
	871	!-- first two significant digits
	872	div = 1000.0
	873	DO WHILE ( dt_1d < div )
	874	div = div / 10.0
	875	ENDDO
	876	dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0
	877
	878	!
	879	!-- Now the time step can be changed.
	880	IF ( percent_change < 0.0 ) THEN
	881	!
	882	!-- Time step reduction
	883	old_dt_1d = dt_1d
	884	last_dt_change_1d = current_timestep_number_1d
	885	ELSE
	886	!
	887	!-- Time step will only be increased if at least 30 iterations have
	888	!-- been done since the previous time step change, and of course at
	889	!-- simulation start, respectively.
	890	IF ( current_timestep_number_1d >= last_dt_change_1d + 30 .OR. &
	891	simulated_time_1d == 0.0 ) THEN
	892	old_dt_1d = dt_1d
	893	last_dt_change_1d = current_timestep_number_1d
	894	ELSE
	895	dt_1d = old_dt_1d
	896	ENDIF
	897	ENDIF
	898	ELSE
	899	!
	900	!-- No time step change since the difference is too small
	901	dt_1d = old_dt_1d
	902	ENDIF
	903
	904	ELSE ! Runge-Kutta
	905
	906	!-- A more or less simple new time step value is obtained taking only the
	907	!-- first two significant digits
	908	div = 1000.0
	909	DO WHILE ( dt_1d < div )
	910	div = div / 10.0
	911	ENDDO
	912	dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0
	913
	914	old_dt_1d = dt_1d
	915	last_dt_change_1d = current_timestep_number_1d
	916
	917	ENDIF
	918
	919	END SUBROUTINE timestep_1d
	920
	921
	922
	923	SUBROUTINE print_1d_model
	924
	925	!------------------------------------------------------------------------------!
	926	! Description:
	927	! ------------
	928	! List output of profiles from the 1D-model
	929	!------------------------------------------------------------------------------!
	930
	931	USE arrays_3d
	932	USE indices
	933	USE model_1d
	934	USE pegrid
	935	USE control_parameters
	936
	937	IMPLICIT NONE
	938
	939
	940	INTEGER :: k
	941
	942
	943	IF ( myid == 0 ) THEN
	944	!
	945	!-- Open list output file for profiles from the 1D-model
	946	CALL check_open( 17 )
	947
	948	!
	949	!-- Write Header
	950	WRITE ( 17, 100 ) TRIM( run_description_header ), &
	951	TRIM( simulated_time_chr )
	952	WRITE ( 17, 101 )
	953
	954	!
	955	!-- Write the values
	956	WRITE ( 17, 102 )
	957	WRITE ( 17, 101 )
	958	DO k = nzt+1, nzb, -1
	959	WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
	960	rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
	961	ENDDO
	962	WRITE ( 17, 101 )
	963	WRITE ( 17, 102 )
	964	WRITE ( 17, 101 )
	965
	966	!
	967	!-- Write buffer contents to disc immediately
[82]	968	CALL local_flush( 17 )
[1]	969
	970	ENDIF
	971
	972	!
	973	!-- Formats
	974	100 FORMAT (//1X,A/1X,10('-')/' 1d-model profiles'/ &
	975	' Time: ',A)
	976	101 FORMAT (1X,79('-'))
	977	102 FORMAT (' k zu u v pt e rif Km Kh ', &
	978	'l zu k')
	979	103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
	980	1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
	981
	982
	983	END SUBROUTINE print_1d_model

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |