Home

Context Navigation

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

Last change on this file since 3302 was 3274, checked in by knoop, 6 years ago
Modularization of all bulk cloud physics code components
Property svn:keywords set to `Id`
File size: 44.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |