Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 35.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |