Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |