Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |