Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |