Home

Context Navigation

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

Last change on this file since 70 was 46, checked in by raasch, 18 years ago
updating parts of Marcus changes
Property svn:keywords set to `Id`
File size: 33.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |