Home

Context Navigation

source: palm/trunk/SOURCE/model_1d_mod.f90 @ 4573

Last change on this file since 4573 was 4449, checked in by suehring, 5 years ago
Bugfix in 1D model: Set intermediate_timestep_count back to zero after 1D-model integration. This is required e.g. for initial calls of calc_mean_profile.
Property svn:keywords set to `Id`
File size: 41.2 KB

Rev	Line
[2338]	1	!> @file model_1d_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	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	!
[4360]	17	! Copyright 1997-2020 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[1961]	22	!
[3049]	23	!
[1961]	24	! Former revisions:
	25	! -----------------
	26	! $Id: model_1d_mod.f90 4449 2020-03-09 14:43:16Z oliver.maas $
[4449]	27	! Set intermediate_timestep_count back to zero after 1D-model integration.
	28	! This is required e.g. for initial calls of calc_mean_profile.
	29	!
	30	! 4360 2020-01-07 11:25:50Z suehring
[4182]	31	! Corrected "Former revisions" section
[4449]	32	!
[4182]	33	! 3655 2019-01-07 16:51:22Z knoop
[3274]	34	! Modularization of all bulk cloud physics code components
[1961]	35	!
[4182]	36	! Revision 1.1 1998/03/09 16:22:10 raasch
	37	! Initial revision
	38	!
	39	!
[1]	40	! Description:
	41	! ------------
[1682]	42	!> 1D-model to initialize the 3D-arrays.
[4449]	43	!> The temperature profile is set as steady and a corresponding steady solution
[1682]	44	!> of the wind profile is being computed.
	45	!> All subroutines required can be found within this file.
[1691]	46	!>
	47	!> @todo harmonize code with new surface_layer_fluxes module
[1709]	48	!> @bug 1D model crashes when using small grid spacings in the order of 1 m
[2965]	49	!> @fixme option "as_in_3d_model" seems to be an inappropriate option because
[2918]	50	!> the 1D model uses different turbulence closure approaches at least if
	51	!> the 3D model is set to LES-mode.
[1]	52	!------------------------------------------------------------------------------!
[2338]	53	MODULE model_1d_mod
[1]	54
[1320]	55	USE arrays_3d, &
[2918]	56	ONLY: dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, &
[2338]	57	vg, v_init, zu
[3274]	58
	59	USE basic_constants_and_equations_mod, &
	60	ONLY: g, kappa, pi
[4449]	61
[2338]	62	USE control_parameters, &
[3274]	63	ONLY: constant_diffusion, constant_flux_layer, dissipation_1d, f, &
[2338]	64	humidity, ibc_e_b, intermediate_timestep_count, &
[3274]	65	intermediate_timestep_count_max, km_constant, &
[2338]	66	message_string, mixing_length_1d, prandtl_number, &
[2696]	67	roughness_length, run_description_header, simulated_time_chr, &
	68	timestep_scheme, tsc, z0h_factor
[2338]	69
[1320]	70	USE indices, &
[2338]	71	ONLY: nzb, nzb_diff, nzt
[4449]	72
[1320]	73	USE kinds
[1]	74
[3083]	75	USE pegrid, &
	76	ONLY: myid
[2338]	77
[4449]	78
[1]	79	IMPLICIT NONE
	80
[2338]	81	INTEGER(iwp) :: current_timestep_number_1d = 0 !< current timestep number (1d-model)
[2696]	82	INTEGER(iwp) :: damp_level_ind_1d !< lower grid index of damping layer (1d-model)
[2338]	83
	84	LOGICAL :: run_control_header_1d = .FALSE. !< flag for output of run control header (1d-model)
	85	LOGICAL :: stop_dt_1d = .FALSE. !< termination flag, used in case of too small timestep (1d-model)
	86
[3083]	87	REAL(wp) :: alpha_buoyancy !< model constant according to Koblitz (2013)
[3126]	88	REAL(wp) :: c_0 = 0.416179145_wp !< = 0.03^0.25; model constant according to Koblitz (2013)
[3083]	89	REAL(wp) :: c_1 = 1.52_wp !< model constant according to Koblitz (2013)
	90	REAL(wp) :: c_2 = 1.83_wp !< model constant according to Koblitz (2013)
	91	REAL(wp) :: c_3 !< model constant
	92	REAL(wp) :: c_mu !< model constant
[2696]	93	REAL(wp) :: damp_level_1d = -1.0_wp !< namelist parameter
[2338]	94	REAL(wp) :: dt_1d = 60.0_wp !< dynamic timestep (1d-model)
	95	REAL(wp) :: dt_max_1d = 300.0_wp !< timestep limit (1d-model)
[2696]	96	REAL(wp) :: dt_pr_1d = 9999999.9_wp !< namelist parameter
	97	REAL(wp) :: dt_run_control_1d = 60.0_wp !< namelist parameter
	98	REAL(wp) :: end_time_1d = 864000.0_wp !< namelist parameter
[3083]	99	REAL(wp) :: lambda !< maximum mixing length
[2338]	100	REAL(wp) :: qs1d !< characteristic humidity scale (1d-model)
	101	REAL(wp) :: simulated_time_1d = 0.0_wp !< updated simulated time (1d-model)
[3083]	102	REAL(wp) :: sig_diss = 2.95_wp !< model constant according to Koblitz (2013)
	103	REAL(wp) :: sig_e = 2.95_wp !< model constant according to Koblitz (2013)
[2338]	104	REAL(wp) :: time_pr_1d = 0.0_wp !< updated simulated time for profile output (1d-model)
	105	REAL(wp) :: time_run_control_1d = 0.0_wp !< updated simulated time for run-control output (1d-model)
	106	REAL(wp) :: ts1d !< characteristic temperature scale (1d-model)
[2696]	107	REAL(wp) :: us1d !< friction velocity (1d-model)
	108	REAL(wp) :: usws1d !< u-component of the momentum flux (1d-model)
	109	REAL(wp) :: vsws1d !< v-component of the momentum flux (1d-model)
[2338]	110	REAL(wp) :: z01d !< roughness length for momentum (1d-model)
	111	REAL(wp) :: z0h1d !< roughness length for scalars (1d-model)
	112
[2696]	113	REAL(wp), DIMENSION(:), ALLOCATABLE :: diss1d !< tke dissipation rate (1d-model)
	114	REAL(wp), DIMENSION(:), ALLOCATABLE :: diss1d_p !< prognostic value of tke dissipation rate (1d-model)
	115	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d !< tke (1d-model)
[2338]	116	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d_p !< prognostic value of tke (1d-model)
[2696]	117	REAL(wp), DIMENSION(:), ALLOCATABLE :: kh1d !< turbulent diffusion coefficient for heat (1d-model)
	118	REAL(wp), DIMENSION(:), ALLOCATABLE :: km1d !< turbulent diffusion coefficient for momentum (1d-model)
	119	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d !< mixing length for turbulent diffusion coefficients (1d-model)
[2918]	120	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_init !< initial mixing length (1d-model)
[2338]	121	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_diss !< mixing length for dissipation (1d-model)
[2696]	122	REAL(wp), DIMENSION(:), ALLOCATABLE :: rif1d !< Richardson flux number (1d-model)
	123	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_diss !< tendency of diss (1d-model)
	124	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model)
[2338]	125	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_e !< tendency of e (1d-model)
	126	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_em !< weighted tendency of e for previous sub-timestep (1d-model)
	127	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_u !< tendency of u (1d-model)
	128	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_um !< weighted tendency of u for previous sub-timestep (1d-model)
	129	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_v !< tendency of v (1d-model)
	130	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_vm !< weighted tendency of v for previous sub-timestep (1d-model)
[2696]	131	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d !< u-velocity component (1d-model)
[2338]	132	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d_p !< prognostic value of u-velocity component (1d-model)
[2696]	133	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d !< v-velocity component (1d-model)
[2338]	134	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d_p !< prognostic value of v-velocity component (1d-model)
	135
	136	!
	137	!-- Initialize 1D model
	138	INTERFACE init_1d_model
	139	MODULE PROCEDURE init_1d_model
	140	END INTERFACE init_1d_model
	141
	142	!
	143	!-- Print profiles
	144	INTERFACE print_1d_model
	145	MODULE PROCEDURE print_1d_model
	146	END INTERFACE print_1d_model
	147
	148	!
	149	!-- Print run control information
	150	INTERFACE run_control_1d
	151	MODULE PROCEDURE run_control_1d
	152	END INTERFACE run_control_1d
	153
	154	!
	155	!-- Main procedure
	156	INTERFACE time_integration_1d
	157	MODULE PROCEDURE time_integration_1d
	158	END INTERFACE time_integration_1d
	159
	160	!
	161	!-- Calculate time step
	162	INTERFACE timestep_1d
	163	MODULE PROCEDURE timestep_1d
	164	END INTERFACE timestep_1d
	165
	166	SAVE
	167
	168	PRIVATE
	169	!
	170	!-- Public interfaces
	171	PUBLIC init_1d_model
	172
	173	!
	174	!-- Public variables
[2696]	175	PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, &
	176	dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, &
	177	us1d, usws1d, v1d, vsws1d
[2338]	178
	179
	180	CONTAINS
	181
	182	SUBROUTINE init_1d_model
[4449]	183
[2918]	184	USE grid_variables, &
	185	ONLY: dx, dy
	186
[2338]	187	IMPLICIT NONE
	188
[2696]	189	CHARACTER (LEN=9) :: time_to_string !< function to transform time from real to character string
	190
	191	INTEGER(iwp) :: k !< loop index
[1]	192
	193	!
	194	!-- Allocate required 1D-arrays
[2696]	195	ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1), &
	196	e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), &
[2918]	197	km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), &
[2696]	198	l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1), &
	199	te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), &
[2338]	200	te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &
	201	te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &
	202	u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
[1]	203
	204	!
	205	!-- Initialize arrays
	206	IF ( constant_diffusion ) THEN
[1001]	207	km1d = km_constant
	208	kh1d = km_constant / prandtl_number
[1]	209	ELSE
[2696]	210	diss1d = 0.0_wp; diss1d_p = 0.0_wp
[1353]	211	e1d = 0.0_wp; e1d_p = 0.0_wp
	212	kh1d = 0.0_wp; km1d = 0.0_wp
	213	rif1d = 0.0_wp
[1]	214	!
	215	!-- Compute the mixing length
[2918]	216	l1d_init(nzb) = 0.0_wp
[1]	217
	218	IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN
	219	!
	220	!-- Blackadar mixing length
[1353]	221	IF ( f /= 0.0_wp ) THEN
	222	lambda = 2.7E-4_wp * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
	223	ABS( f ) + 1E-10_wp
[1]	224	ELSE
[1353]	225	lambda = 30.0_wp
[1]	226	ENDIF
	227
	228	DO k = nzb+1, nzt+1
[2918]	229	l1d_init(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
[1]	230	ENDDO
	231
	232	ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN
	233	!
[2918]	234	!-- Use the same mixing length as in 3D model (LES-mode)
[3083]	235	!> @todo rename (delete?) this option
	236	!> As the mixing length is different between RANS and LES mode, it
	237	!> must be distinguished here between these modes. For RANS mode,
	238	!> the mixing length is calculated accoding to Blackadar, which is
	239	!> the other option at this point.
	240	!> Maybe delete this option entirely (not appropriate in LES case)
	241	!> 2018-03-20, gronemeier
[2918]	242	DO k = nzb+1, nzt
	243	l1d_init(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
	244	ENDDO
	245	l1d_init(nzt+1) = l1d_init(nzt)
[1]	246
	247	ENDIF
	248	ENDIF
[2918]	249	l1d = l1d_init
	250	l1d_diss = l1d_init
[2337]	251	u1d = u_init
	252	u1d_p = u_init
	253	v1d = v_init
	254	v1d_p = v_init
[1]	255
	256	!
	257	!-- Set initial horizontal velocities at the lowest grid levels to a very small
	258	!-- value in order to avoid too small time steps caused by the diffusion limit
	259	!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
[1353]	260	u1d(0:1) = 0.1_wp
	261	u1d_p(0:1) = 0.1_wp
	262	v1d(0:1) = 0.1_wp
	263	v1d_p(0:1) = 0.1_wp
[1]	264
	265	!
	266	!-- For u, theta and the momentum fluxes plausible values are set
[1691]	267	IF ( constant_flux_layer ) THEN
[1353]	268	us1d = 0.1_wp ! without initial friction the flow would not change
[1]	269	ELSE
[3083]	270	diss1d(nzb+1) = 0.001_wp
[1353]	271	e1d(nzb+1) = 1.0_wp
	272	km1d(nzb+1) = 1.0_wp
	273	us1d = 0.0_wp
[1]	274	ENDIF
[1353]	275	ts1d = 0.0_wp
	276	usws1d = 0.0_wp
	277	vsws1d = 0.0_wp
[996]	278	z01d = roughness_length
[4449]	279	z0h1d = z0h_factor * z01d
[1960]	280	IF ( humidity ) qs1d = 0.0_wp
[1]	281
	282	!
[3083]	283	!-- Tendencies must be preset in order to avoid runtime errors
	284	te_diss = 0.0_wp
[2696]	285	te_dissm = 0.0_wp
[3083]	286	te_e = 0.0_wp
[1353]	287	te_em = 0.0_wp
	288	te_um = 0.0_wp
	289	te_vm = 0.0_wp
[46]	290
	291	!
[2338]	292	!-- Set model constant
[3083]	293	IF ( dissipation_1d == 'as_in_3d_model' ) c_0 = 0.1_wp
	294	c_mu = c_0**4
[2338]	295
	296	!
[1]	297	!-- Set start time in hh:mm:ss - format
	298	simulated_time_chr = time_to_string( simulated_time_1d )
	299
	300	!
[2337]	301	!-- Integrate the 1D-model equations using the Runge-Kutta scheme
[1]	302	CALL time_integration_1d
	303
	304
	305	END SUBROUTINE init_1d_model
	306
	307
	308
	309	!------------------------------------------------------------------------------!
	310	! Description:
	311	! ------------
[2338]	312	!> Runge-Kutta time differencing scheme for the 1D-model.
[1]	313	!------------------------------------------------------------------------------!
[4449]	314
[1682]	315	SUBROUTINE time_integration_1d
[1]	316
	317	IMPLICIT NONE
	318
[2696]	319	CHARACTER (LEN=9) :: time_to_string !< function to transform time from real to character string
	320
[2338]	321	INTEGER(iwp) :: k !< loop index
[2696]	322
[2338]	323	REAL(wp) :: a !< auxiliary variable
	324	REAL(wp) :: b !< auxiliary variable
	325	REAL(wp) :: dpt_dz !< vertical temperature gradient
	326	REAL(wp) :: flux !< vertical temperature gradient
	327	REAL(wp) :: kmzm !< Km(z-dz/2)
	328	REAL(wp) :: kmzp !< Km(z+dz/2)
	329	REAL(wp) :: l_stable !< mixing length for stable case
	330	REAL(wp) :: pt_0 !< reference temperature
	331	REAL(wp) :: uv_total !< horizontal wind speed
[1]	332
	333	!
	334	!-- Determine the time step at the start of a 1D-simulation and
	335	!-- determine and printout quantities used for run control
[3083]	336	dt_1d = 0.01_wp
[1]	337	CALL run_control_1d
	338
	339	!
	340	!-- Start of time loop
	341	DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d )
	342
	343	!
	344	!-- Depending on the timestep scheme, carry out one or more intermediate
	345	!-- timesteps
	346
	347	intermediate_timestep_count = 0
	348	DO WHILE ( intermediate_timestep_count < &
	349	intermediate_timestep_count_max )
	350
	351	intermediate_timestep_count = intermediate_timestep_count + 1
	352
	353	CALL timestep_scheme_steering
	354
	355	!
[4449]	356	!-- Compute all tendency terms. If a constant-flux layer is simulated,
[2696]	357	!-- k starts at nzb+2.
[1]	358	DO k = nzb_diff, nzt
	359
[1353]	360	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	361	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[1]	362	!
	363	!-- u-component
	364	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
[1001]	365	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
	366	- kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) &
	367	) * ddzw(k)
[1]	368	!
	369	!-- v-component
[1001]	370	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	371	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
	372	- kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) &
	373	) * ddzw(k)
[1]	374	ENDDO
	375	IF ( .NOT. constant_diffusion ) THEN
	376	DO k = nzb_diff, nzt
	377	!
[2696]	378	!-- TKE and dissipation rate
[1353]	379	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	380	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[75]	381	IF ( .NOT. humidity ) THEN
[1]	382	pt_0 = pt_init(k)
	383	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	384	ELSE
[1353]	385	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
	386	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
[2337]	387	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
	388	pt_init(k-1) * q_init(k-1) ) &
	389	) * dd2zu(k)
[1]	390	ENDIF
	391
[2696]	392	!
	393	!-- Calculate dissipation rate if no prognostic equation is used for
	394	!-- dissipation rate
[1]	395	IF ( dissipation_1d == 'detering' ) THEN
[3083]	396	diss1d(k) = c_0*3 e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	397	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
[2918]	398	diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
[2696]	399	) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	400	ENDIF
[2696]	401	!
	402	!-- TKE
[1]	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) &
[3083]	410	) * ddzw(k) / sig_e &
[2696]	411	- diss1d(k)
	412
	413	IF ( dissipation_1d == 'prognostic' ) THEN
	414	!
	415	!-- dissipation rate
[3083]	416	IF ( rif1d(k) >= 0.0_wp ) THEN
	417	alpha_buoyancy = 1.0_wp - l1d(k) / lambda
	418	ELSE
	419	alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp ) &
	420	/ ( c_2 - c_1 ) ) &
	421	* l1d(k) / lambda
	422	ENDIF
[3135]	423	c_3 = ( c_1 - c_2 ) * alpha_buoyancy + 1.0_wp
[3083]	424	te_diss(k) = ( km1d(k) * &
[2696]	425	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	426	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
[3083]	427	) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda ) &
[2696]	428	- g / pt_0 * kh1d(k) * flux * c_3 &
[3083]	429	- c_2 * diss1d(k) &
	430	) * diss1d(k) / ( e1d(k) + 1.0E-20_wp ) &
	431	+ ( kmzp * ( diss1d(k+1) - diss1d(k) ) &
[2696]	432	* ddzu(k+1) &
	433	- kmzm * ( diss1d(k) - diss1d(k-1) ) &
	434	* ddzu(k) &
[3083]	435	) * ddzw(k) / sig_diss
[2696]	436
	437	ENDIF
	438
[1]	439	ENDDO
	440	ENDIF
	441
	442	!
[2696]	443	!-- Tendency terms at the top of the constant-flux layer.
[4449]	444	!-- Finite differences of the momentum fluxes are computed using half the
	445	!-- normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
[1691]	446	IF ( constant_flux_layer ) THEN
[1]	447
	448	k = nzb+1
[1353]	449	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
	450	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
[75]	451	IF ( .NOT. humidity ) THEN
[1]	452	pt_0 = pt_init(k)
	453	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
	454	ELSE
[1353]	455	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
	456	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
[2337]	457	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
	458	pt_init(k-1) * q_init(k-1) ) &
[1]	459	) * dd2zu(k)
	460	ENDIF
	461
[2696]	462	!
	463	!-- Calculate dissipation rate if no prognostic equation is used for
	464	!-- dissipation rate
[1]	465	IF ( dissipation_1d == 'detering' ) THEN
[3083]	466	diss1d(k) = c_0*3 e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	467	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
[2918]	468	diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) &
[2696]	469	* e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
[1]	470	ENDIF
	471
	472	!
	473	!-- u-component
[1001]	474	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
	475	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d &
[1353]	476	) * 2.0_wp * ddzw(k)
[1]	477	!
	478	!-- v-component
[1001]	479	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
	480	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d &
[1353]	481	) * 2.0_wp * ddzw(k)
[1]	482	!
	483	!-- TKE
[2696]	484	IF ( .NOT. dissipation_1d == 'prognostic' ) THEN
[3083]	485	!> @query why integrate over 2dz
	486	!> Why is it allowed to integrate over two delta-z for e
	487	!> while for u and v it is not?
	488	!> 2018-04-23, gronemeier
[2696]	489	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
	490	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
	491	) &
	492	- g / pt_0 * kh1d(k) * flux &
	493	+ ( &
	494	kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &
	495	- kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &
[3083]	496	) * ddzw(k) / sig_e &
[2696]	497	- diss1d(k)
	498	ENDIF
	499
[1]	500	ENDIF
	501
	502	!
	503	!-- Prognostic equations for all 1D variables
	504	DO k = nzb+1, nzt
	505
[1001]	506	u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
	507	tsc(3) * te_um(k) )
	508	v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
	509	tsc(3) * te_vm(k) )
[1]	510
	511	ENDDO
	512	IF ( .NOT. constant_diffusion ) THEN
[2696]	513
[1]	514	DO k = nzb+1, nzt
[1001]	515	e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
	516	tsc(3) * te_em(k) )
[2696]	517	ENDDO
[1]	518
	519	!
	520	!-- Eliminate negative TKE values, which can result from the
	521	!-- integration due to numerical inaccuracies. In such cases the TKE
	522	!-- value is reduced to 10 percent of its old value.
[1353]	523	WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d
[3083]	524
	525	IF ( dissipation_1d == 'prognostic' ) THEN
	526	DO k = nzb+1, nzt
	527	diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
	528	tsc(3) * te_dissm(k) )
	529	ENDDO
	530	WHERE ( diss1d_p < 0.0_wp ) diss1d_p = 0.1_wp * diss1d
	531	ENDIF
[1]	532	ENDIF
	533
	534	!
	535	!-- Calculate tendencies for the next Runge-Kutta step
	536	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	537	IF ( intermediate_timestep_count == 1 ) THEN
	538
	539	DO k = nzb+1, nzt
	540	te_um(k) = te_u(k)
	541	te_vm(k) = te_v(k)
	542	ENDDO
	543
	544	IF ( .NOT. constant_diffusion ) THEN
	545	DO k = nzb+1, nzt
	546	te_em(k) = te_e(k)
	547	ENDDO
[2696]	548	IF ( dissipation_1d == 'prognostic' ) THEN
	549	DO k = nzb+1, nzt
	550	te_dissm(k) = te_diss(k)
	551	ENDDO
	552	ENDIF
[1]	553	ENDIF
	554
	555	ELSEIF ( intermediate_timestep_count < &
	556	intermediate_timestep_count_max ) THEN
	557
	558	DO k = nzb+1, nzt
[1353]	559	te_um(k) = -9.5625_wp * te_u(k) + 5.3125_wp * te_um(k)
	560	te_vm(k) = -9.5625_wp * te_v(k) + 5.3125_wp * te_vm(k)
[1]	561	ENDDO
	562
	563	IF ( .NOT. constant_diffusion ) THEN
	564	DO k = nzb+1, nzt
[1353]	565	te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
[1]	566	ENDDO
[2696]	567	IF ( dissipation_1d == 'prognostic' ) THEN
	568	DO k = nzb+1, nzt
[3083]	569	te_dissm(k) = -9.5625_wp * te_diss(k) &
	570	+ 5.3125_wp * te_dissm(k)
[2696]	571	ENDDO
	572	ENDIF
[1]	573	ENDIF
	574
	575	ENDIF
	576	ENDIF
	577
	578	!
	579	!-- Boundary conditions for the prognostic variables.
[4449]	580	!-- At the top boundary (nzt+1) u, v, e, and diss keep their initial
[2696]	581	!-- values (ug(nzt+1), vg(nzt+1), 0, 0).
[2334]	582	!-- At the bottom boundary, Dirichlet condition is used for u and v (0)
[2696]	583	!-- and Neumann condition for e and diss (e(nzb)=e(nzb+1)).
[1353]	584	u1d_p(nzb) = 0.0_wp
	585	v1d_p(nzb) = 0.0_wp
[667]	586
[1]	587	!
	588	!-- Swap the time levels in preparation for the next time step.
	589	u1d = u1d_p
	590	v1d = v1d_p
	591	IF ( .NOT. constant_diffusion ) THEN
	592	e1d = e1d_p
[2696]	593	IF ( dissipation_1d == 'prognostic' ) THEN
	594	diss1d = diss1d_p
	595	ENDIF
[1]	596	ENDIF
	597
	598	!
	599	!-- Compute diffusion quantities
	600	IF ( .NOT. constant_diffusion ) THEN
	601
	602	!
[2696]	603	!-- First compute the vertical fluxes in the constant-flux layer
[1691]	604	IF ( constant_flux_layer ) THEN
[1]	605	!
	606	!-- Compute theta* using Rif numbers of the previous time step
[2334]	607	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	608	!
	609	!-- Stable stratification
[1353]	610	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	611	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
	612	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
[1]	613	)
	614	ELSE
	615	!
	616	!-- Unstable stratification
[1353]	617	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
	618	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	619	zu(nzb+1) * z0h1d )
[2337]	620
	621	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
	622	LOG( (a-1.0_wp) / (a+1.0_wp) * &
	623	(b+1.0_wp) / (b-1.0_wp) )
[1]	624	ENDIF
	625
[1691]	626	ENDIF ! constant_flux_layer
[3083]	627	!> @todo combine if clauses
	628	!> The previous and following if clauses can be combined into a
	629	!> single clause
	630	!> 2018-04-23, gronemeier
[1]	631	!
	632	!-- Compute the Richardson-flux numbers,
[4449]	633	!-- first at the top of the constant-flux layer using u* of the
[2696]	634	!-- previous time step (+1E-30, if u* = 0), then in the remaining area.
	635	!-- There the rif-numbers of the previous time step are used.
[1]	636
[1691]	637	IF ( constant_flux_layer ) THEN
[75]	638	IF ( .NOT. humidity ) THEN
[1]	639	pt_0 = pt_init(nzb+1)
	640	flux = ts1d
	641	ELSE
[1353]	642	pt_0 = pt_init(nzb+1) * ( 1.0_wp + 0.61_wp * q_init(nzb+1) )
	643	flux = ts1d + 0.61_wp * pt_init(k) * qs1d
[1]	644	ENDIF
	645	rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
[1353]	646	( pt_0 * ( us1d**2 + 1E-30_wp ) )
[1]	647	ENDIF
	648
	649	DO k = nzb_diff, nzt
[75]	650	IF ( .NOT. humidity ) THEN
[1]	651	pt_0 = pt_init(k)
	652	flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
	653	ELSE
[1353]	654	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
[1]	655	flux = ( ( pt_init(k+1) - pt_init(k-1) ) &
[2337]	656	+ 0.61_wp &
	657	* ( pt_init(k+1) * q_init(k+1) &
	658	- pt_init(k-1) * q_init(k-1) ) &
[1]	659	) * dd2zu(k)
	660	ENDIF
[1353]	661	IF ( rif1d(k) >= 0.0_wp ) THEN
	662	rif1d(k) = g / pt_0 * flux / &
	663	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	664	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	665	+ 1E-30_wp &
[1]	666	)
	667	ELSE
[1353]	668	rif1d(k) = g / pt_0 * flux / &
	669	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
	670	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
	671	+ 1E-30_wp &
	672	) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
[1]	673	ENDIF
	674	ENDDO
	675	!
	676	!-- Richardson-numbers must remain restricted to a realistic value
	677	!-- range. It is exceeded excessively for very small velocities
	678	!-- (u,v --> 0).
[2059]	679	WHERE ( rif1d < -5.0_wp ) rif1d = -5.0_wp
	680	WHERE ( rif1d > 1.0_wp ) rif1d = 1.0_wp
[1]	681
	682	!
	683	!-- Compute u* from the absolute velocity value
[1691]	684	IF ( constant_flux_layer ) THEN
[1]	685	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
	686
[1353]	687	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	688	!
	689	!-- Stable stratification
	690	us1d = kappa * uv_total / ( &
[1353]	691	LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * &
[1]	692	( zu(nzb+1) - z01d ) / zu(nzb+1) &
	693	)
	694	ELSE
	695	!
	696	!-- Unstable stratification
[1353]	697	a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) )
	698	b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	699	zu(nzb+1) * z01d ) )
[2337]	700	us1d = kappa * uv_total / ( &
	701	LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / &
	702	(1.0_wp+a) ) + &
	703	2.0_wp * ( ATAN( b ) - ATAN( a ) ) &
	704	)
[1]	705	ENDIF
	706
	707	!
	708	!-- Compute the momentum fluxes for the diffusion terms
	709	usws1d = - u1d(nzb+1) / uv_total * us1d**2
	710	vsws1d = - v1d(nzb+1) / uv_total * us1d**2
	711
	712	!
[4449]	713	!-- Boundary condition for the turbulent kinetic energy and
[2696]	714	!-- dissipation rate at the top of the constant-flux layer.
[1]	715	!-- Additional Neumann condition de/dz = 0 at nzb is set to ensure
	716	!-- compatibility with the 3D model.
	717	IF ( ibc_e_b == 2 ) THEN
[3083]	718	e1d(nzb+1) = ( us1d / c_0 )**2
[1]	719	ENDIF
[2696]	720	IF ( dissipation_1d == 'prognostic' ) THEN
[3083]	721	e1d(nzb+1) = ( us1d / c_0 )**2
[2696]	722	diss1d(nzb+1) = us1d*3 / ( kappa zu(nzb+1) )
	723	diss1d(nzb) = diss1d(nzb+1)
	724	ENDIF
[1]	725	e1d(nzb) = e1d(nzb+1)
	726
[1960]	727	IF ( humidity ) THEN
[1]	728	!
[4449]	729	!-- Compute q*
[2334]	730	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
[1]	731	!
[1960]	732	!-- Stable stratification
	733	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
[1353]	734	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
	735	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
[1]	736	)
[1960]	737	ELSE
[1]	738	!
[1960]	739	!-- Unstable stratification
	740	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
	741	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
	742	zu(nzb+1) * z0h1d )
[2337]	743	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
	744	LOG( (a-1.0_wp) / (a+1.0_wp) * &
	745	(b+1.0_wp) / (b-1.0_wp) )
	746	ENDIF
[1]	747	ELSE
[1353]	748	qs1d = 0.0_wp
[2337]	749	ENDIF
[1]	750
[1691]	751	ENDIF ! constant_flux_layer
[1]	752
	753	!
[2337]	754	!-- Compute the diabatic mixing length. The unstable stratification
	755	!-- must not be considered for l1d (km1d) as it is already considered
	756	!-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
	757	!-- too large.
[3083]	758	IF ( dissipation_1d /= 'prognostic' ) THEN
	759	IF ( mixing_length_1d == 'blackadar' ) THEN
	760	DO k = nzb+1, nzt
	761	IF ( rif1d(k) >= 0.0_wp ) THEN
	762	l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
	763	l1d_diss(k) = l1d(k)
	764	ELSE
	765	l1d(k) = l1d_init(k)
	766	l1d_diss(k) = l1d_init(k) * &
	767	SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
	768	ENDIF
	769	ENDDO
	770	ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN
	771	DO k = nzb+1, nzt
	772	dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
	773	IF ( dpt_dz > 0.0_wp ) THEN
	774	l_stable = 0.76_wp * SQRT( e1d(k) ) &
	775	/ SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
	776	ELSE
	777	l_stable = l1d_init(k)
	778	ENDIF
	779	l1d(k) = MIN( l1d_init(k), l_stable )
[2337]	780	l1d_diss(k) = l1d(k)
[3083]	781	ENDDO
	782	ENDIF
	783	ELSE
[1]	784	DO k = nzb+1, nzt
[3083]	785	l1d(k) = c_0*3 e1d(k) * SQRT( e1d(k) ) &
	786	/ ( diss1d(k) + 1.0E-30_wp )
[1]	787	ENDDO
	788	ENDIF
	789
	790	!
	791	!-- Compute the diffusion coefficients for momentum via the
	792	!-- corresponding Prandtl-layer relationship and according to
[2337]	793	!-- Prandtl-Kolmogorov, respectively
[1691]	794	IF ( constant_flux_layer ) THEN
[1353]	795	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
	796	km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &
	797	( 1.0_wp + 5.0_wp * rif1d(nzb+1) )
[1]	798	ELSE
[1353]	799	km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &
	800	( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp
[1]	801	ENDIF
	802	ENDIF
[3083]	803
[2696]	804	IF ( dissipation_1d == 'prognostic' ) THEN
	805	DO k = nzb_diff, nzt
	806	km1d(k) = c_mu * e1d(k)**2 / ( diss1d(k) + 1.0E-30_wp )
	807	ENDDO
	808	ELSE
	809	DO k = nzb_diff, nzt
[3083]	810	km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k)
[2696]	811	ENDDO
	812	ENDIF
[1]	813
	814	!
	815	!-- Add damping layer
	816	DO k = damp_level_ind_1d+1, nzt+1
[1353]	817	km1d(k) = 1.1_wp * km1d(k-1)
[1346]	818	km1d(k) = MIN( km1d(k), 10.0_wp )
[1]	819	ENDDO
	820
	821	!
	822	!-- Compute the diffusion coefficient for heat via the relationship
	823	!-- kh = phim / phih * km
	824	DO k = nzb+1, nzt
[1353]	825	IF ( rif1d(k) >= 0.0_wp ) THEN
[1]	826	kh1d(k) = km1d(k)
	827	ELSE
[1353]	828	kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
[1]	829	ENDIF
	830	ENDDO
	831
	832	ENDIF ! .NOT. constant_diffusion
	833
	834	ENDDO ! intermediate step loop
	835
	836	!
	837	!-- Increment simulated time and output times
	838	current_timestep_number_1d = current_timestep_number_1d + 1
	839	simulated_time_1d = simulated_time_1d + dt_1d
	840	simulated_time_chr = time_to_string( simulated_time_1d )
	841	time_pr_1d = time_pr_1d + dt_1d
	842	time_run_control_1d = time_run_control_1d + dt_1d
	843
	844	!
	845	!-- Determine and print out quantities for run control
	846	IF ( time_run_control_1d >= dt_run_control_1d ) THEN
	847	CALL run_control_1d
	848	time_run_control_1d = time_run_control_1d - dt_run_control_1d
	849	ENDIF
	850
	851	!
	852	!-- Profile output on file
	853	IF ( time_pr_1d >= dt_pr_1d ) THEN
	854	CALL print_1d_model
	855	time_pr_1d = time_pr_1d - dt_pr_1d
	856	ENDIF
	857
	858	!
	859	!-- Determine size of next time step
	860	CALL timestep_1d
	861
	862	ENDDO ! time loop
[4449]	863	!
	864	!-- Set intermediate_timestep_count back to zero. This is required e.g. for
	865	!-- initial calls of calc_mean_profile.
	866	intermediate_timestep_count = 0
[1]	867
	868	END SUBROUTINE time_integration_1d
	869
	870
	871	!------------------------------------------------------------------------------!
	872	! Description:
	873	! ------------
[1682]	874	!> Compute and print out quantities for run control of the 1D model.
[1]	875	!------------------------------------------------------------------------------!
[4449]	876
[1682]	877	SUBROUTINE run_control_1d
[1]	878
[1682]	879
[1]	880	IMPLICIT NONE
	881
[2338]	882	INTEGER(iwp) :: k !< loop index
[4449]	883
[2338]	884	REAL(wp) :: alpha !< angle of wind vector at top of constant-flux layer
	885	REAL(wp) :: energy !< kinetic energy
	886	REAL(wp) :: umax !< maximum of u
	887	REAL(wp) :: uv_total !< horizontal wind speed
	888	REAL(wp) :: vmax !< maximum of v
[1]	889
	890	!
	891	!-- Output
	892	IF ( myid == 0 ) THEN
	893	!
	894	!-- If necessary, write header
	895	IF ( .NOT. run_control_header_1d ) THEN
[184]	896	CALL check_open( 15 )
[1]	897	WRITE ( 15, 100 )
	898	run_control_header_1d = .TRUE.
	899	ENDIF
	900
	901	!
	902	!-- Compute control quantities
	903	!-- grid level nzp is excluded due to mirror boundary condition
[1353]	904	umax = 0.0_wp; vmax = 0.0_wp; energy = 0.0_wp
[1]	905	DO k = nzb+1, nzt+1
	906	umax = MAX( ABS( umax ), ABS( u1d(k) ) )
	907	vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
[1353]	908	energy = energy + 0.5_wp * ( u1d(k)2 + v1d(k)2 )
[1]	909	ENDDO
[1322]	910	energy = energy / REAL( nzt - nzb + 1, KIND=wp )
[1]	911
	912	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
[1691]	913	IF ( ABS( v1d(nzb+1) ) < 1.0E-5_wp ) THEN
[1346]	914	alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) )
[1]	915	ELSE
	916	alpha = ACOS( u1d(nzb+1) / uv_total )
[1353]	917	IF ( v1d(nzb+1) <= 0.0_wp ) alpha = 2.0_wp * pi - alpha
[1]	918	ENDIF
[1353]	919	alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
[1]	920
	921	WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &
	922	dt_1d, umax, vmax, us1d, alpha, energy
	923	!
	924	!-- Write buffer contents to disc immediately
[1808]	925	FLUSH( 15 )
[1]	926
	927	ENDIF
	928
	929	!
	930	!-- formats
[2299]	931	100 FORMAT (///'1D run control output:'/ &
[1]	932	&'------------------------------'// &
[2965]	933	&'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ &
[1]	934	&'-------------------------------------------------------------')
[2965]	935	101 FORMAT (I7,1X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
[1]	936
	937
	938	END SUBROUTINE run_control_1d
	939
	940
	941
	942	!------------------------------------------------------------------------------!
	943	! Description:
	944	! ------------
[1682]	945	!> Compute the time step w.r.t. the diffusion criterion
[1]	946	!------------------------------------------------------------------------------!
[4449]	947
[1682]	948	SUBROUTINE timestep_1d
[1]	949
	950	IMPLICIT NONE
	951
[2338]	952	INTEGER(iwp) :: k !< loop index
[1]	953
[2338]	954	REAL(wp) :: dt_diff !< time step accorind to diffusion criterion
[3083]	955	REAL(wp) :: dt_old !< previous time step
[2338]	956	REAL(wp) :: fac !< factor of criterion
	957	REAL(wp) :: value !< auxiliary variable
[1]	958
	959	!
[3083]	960	!-- Save previous time step
	961	dt_old = dt_1d
	962
	963	!
[1]	964	!-- Compute the currently feasible time step according to the diffusion
	965	!-- criterion. At nzb+1 the half grid length is used.
[3083]	966	fac = 0.125
[1]	967	dt_diff = dt_max_1d
	968	DO k = nzb+2, nzt
[1353]	969	value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20_wp )
[1]	970	dt_diff = MIN( value, dt_diff )
	971	ENDDO
[1353]	972	value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20_wp )
[1]	973	dt_1d = MIN( value, dt_diff )
	974
	975	!
[3083]	976	!-- Limit the new time step to a maximum of 10 times the previous time step
	977	dt_1d = MIN( dt_old * 10.0_wp, dt_1d )
	978
	979	!
[1]	980	!-- Set flag when the time step becomes too small
[3083]	981	IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) ) THEN
[1]	982	stop_dt_1d = .TRUE.
[254]	983
[3046]	984	WRITE( message_string, * ) 'timestep has exceeded the lower limit&', &
[254]	985	'dt_1d = ',dt_1d,' s simulation stopped!'
	986	CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
[4449]	987
[1]	988	ENDIF
	989
	990	END SUBROUTINE timestep_1d
	991
	992
	993
	994	!------------------------------------------------------------------------------!
	995	! Description:
	996	! ------------
[1682]	997	!> List output of profiles from the 1D-model
[1]	998	!------------------------------------------------------------------------------!
[4449]	999
[1682]	1000	SUBROUTINE print_1d_model
[1]	1001
	1002	IMPLICIT NONE
	1003
[2338]	1004	INTEGER(iwp) :: k !< loop parameter
[1]	1005
[2338]	1006	LOGICAL, SAVE :: write_first = .TRUE. !< flag for writing header
[1]	1007
	1008
	1009	IF ( myid == 0 ) THEN
	1010	!
	1011	!-- Open list output file for profiles from the 1D-model
	1012	CALL check_open( 17 )
	1013
	1014	!
	1015	!-- Write Header
[2338]	1016	IF ( write_first ) THEN
	1017	WRITE ( 17, 100 ) TRIM( run_description_header )
	1018	write_first = .FALSE.
	1019	ENDIF
[1]	1020
	1021	!
	1022	!-- Write the values
[2338]	1023	WRITE ( 17, 104 ) TRIM( simulated_time_chr )
	1024	WRITE ( 17, 101 )
[1]	1025	WRITE ( 17, 102 )
	1026	WRITE ( 17, 101 )
	1027	DO k = nzt+1, nzb, -1
	1028	WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
[2696]	1029	rif1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
[1]	1030	ENDDO
	1031	WRITE ( 17, 101 )
	1032	WRITE ( 17, 102 )
	1033	WRITE ( 17, 101 )
	1034
	1035	!
	1036	!-- Write buffer contents to disc immediately
[1808]	1037	FLUSH( 17 )
[1]	1038
	1039	ENDIF
	1040
	1041	!
	1042	!-- Formats
[2338]	1043	100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
	1044	104 FORMAT (//'# Time: ',A)
[2696]	1045	101 FORMAT ('#',111('-'))
	1046	102 FORMAT ('# k zu u v pt e ', &
	1047	'rif Km Kh l diss ')
	1048	103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3))
[1]	1049
	1050
	1051	END SUBROUTINE print_1d_model
[2338]	1052
	1053
[3045]	1054	END MODULE

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |