Home

Context Navigation

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

Last change on this file since 690 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 34.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |