Home

Context Navigation

init_1d_model.f90 @ 100

Last change on this file since 100 was 83, checked in by raasch, 18 years ago

New:
---

Changed:

PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:

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

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

Context Navigation

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

Download in other formats: