Home

Context Navigation

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

Last change on this file since 1332 was 1323, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 34.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |