Home

Context Navigation

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

Last change on this file since 2304 was 2299, checked in by maronga, 7 years ago
improvements for spinup mechanism
Property svn:keywords set to `Id`
File size: 36.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |