Home

Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

model_1d_mod.f90

Last change on this file was 4828, checked in by Giersch, 3 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 42.5 KB

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

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

Context Navigation

source: palm/trunk/SOURCE/model_1d_mod.f90

Download in other formats: