Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |