Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |