Home

Context Navigation

init_1d_model.f90 @ 2337

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

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

Context Navigation

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

Download in other formats: