Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |