Home

Context Navigation

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

Last change on this file since 331 was 254, checked in by heinze, 16 years ago
Output of messages replaced by message handling routine.
Property svn:keywords set to `Id`
File size: 34.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |