Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |