Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |