Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |