Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |