Home

Context Navigation

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

Last change on this file since 3120 was 3083, checked in by gronemeier, 6 years ago
merge with branch rans
Property svn:keywords set to `Id`
File size: 44.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |