Home

Context Navigation

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

Last change on this file since 17 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 33.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |