Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |