Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |