Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |