Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |