Home

Context Navigation

source: palm/trunk/SOURCE/vdi_internal_controls.f90 @ 4343

Last change on this file since 4343 was 4329, checked in by motisi, 5 years ago
Renamed wall_flags_0 to wall_flags_static_0
Property svn:keywords set to `Id`
File size: 39.3 KB

Rev	Line
[4173]	1	!> @file vdi_internal_controls.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 2019-2019 Leibniz Universitaet Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
[4175]	22	!
	23	!
[4173]	24	! Former revisions:
	25	! -----------------
	26	! $Id: vdi_internal_controls.f90 4329 2019-12-10 15:46:36Z oliver.maas $
[4329]	27	! Renamed wall_flags_0 to wall_flags_static_0
	28	!
	29	! 4182 2019-08-22 15:20:23Z scharf
[4182]	30	! added "Authors" section
	31	!
	32	! 4175 2019-08-20 13:19:16Z gronemeier
[4175]	33	! bugfix: removed unused variables
	34	!
	35	! 4173 2019-08-20 12:04:06Z weniger
[4173]	36	! Initial version
	37	!
[4182]	38	! Authors:
	39	! --------
	40	! @author Viola Weniger
[4173]	41	!
	42	!
	43	! Description:
	44	! ------------
	45	!> According to VDI Guideline 3783 Part 9, internal assessment have to be
	46	!> carried out within the program for the model to be considered as evaluated.
	47	!------------------------------------------------------------------------------!
	48	MODULE vdi_internal_controls
	49
	50	USE arrays_3d, &
	51	ONLY: dzw, &
	52	pt, &
	53	q, &
	54	u, &
	55	u_p, &
	56	v, &
	57	w
	58
	59	USE control_parameters, &
	60	ONLY: bc_dirichlet_l, &
	61	bc_dirichlet_n, &
	62	bc_dirichlet_r, &
	63	bc_dirichlet_s, &
	64	bc_lr_cyc, &
	65	bc_ns_cyc, &
	66	humidity, &
	67	end_time, &
	68	message_string, &
	69	neutral, &
	70	time_since_reference_point
	71
	72	USE indices, &
	73	ONLY: nx, &
	74	nxl, &
	75	nxlg, &
	76	nxr, &
	77	nxrg, &
	78	ny, &
	79	nyn, &
	80	nyng, &
	81	nys, &
	82	nysg, &
	83	nzb, &
	84	nzt, &
[4329]	85	wall_flags_static_0
[4173]	86
	87	USE kinds
	88
	89	USE pegrid, &
	90	ONLY: collective_wait, &
	91	comm2d, &
	92	ierr, &
	93	MPI_DOUBLE_PRECISION, &
	94	MPI_INTEGER, &
	95	MPI_MAX, &
	96	MPI_SUM, &
	97	myid
	98
	99
	100	USE grid_variables, &
	101	ONLY: dx, &
	102	dy
	103
	104	USE pmc_interface, &
	105	ONLY: nested_run
	106
	107	IMPLICIT NONE
	108
	109	INTEGER(iwp) :: internal_count = 0 !< counts calls to this module
	110
	111	INTERFACE vdi_2_deltat_wave
	112	MODULE PROCEDURE vdi_2_deltat_wave
	113	END INTERFACE vdi_2_deltat_wave
	114
	115	INTERFACE vdi_standard_differences
	116	MODULE PROCEDURE vdi_standard_differences
	117	END INTERFACE vdi_standard_differences
	118
	119	INTERFACE vdi_domain_averages
	120	MODULE PROCEDURE vdi_domain_averages
	121	END INTERFACE vdi_domain_averages
	122
	123	INTERFACE vdi_plausible_values
	124	MODULE PROCEDURE vdi_plausible_values
	125	END INTERFACE vdi_plausible_values
	126
	127	INTERFACE vdi_actions
	128	MODULE PROCEDURE vdi_actions
	129	END INTERFACE vdi_actions
	130
	131	INTERFACE vdi_conservation_of_mass
	132	MODULE PROCEDURE vdi_conservation_of_mass
	133	END INTERFACE vdi_conservation_of_mass
	134
	135	SAVE
	136
	137	PRIVATE
	138
	139	!
	140	!-- Public functions
[4175]	141	PUBLIC &
[4173]	142	vdi_actions
	143
	144
	145	CONTAINS
	146
	147	!------------------------------------------------------------------------------!
	148	! Description:
	149	! ------------
	150	!> Call for all grid points
	151	!> @todo Add proper description
	152	!------------------------------------------------------------------------------!
	153	SUBROUTINE vdi_actions( location )
	154
	155	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
	156
	157
	158	SELECT CASE ( location )
	159
	160	CASE ( 'after_integration' )
	161
	162	internal_count = internal_count + 1
	163
	164	CALL vdi_2_deltat_wave
	165
	166	CALL vdi_standard_differences
	167
	168	CALL vdi_domain_averages
	169
	170	CALL vdi_conservation_of_mass
	171
	172	CALL vdi_plausible_values
	173
	174	CASE DEFAULT
	175
	176	CONTINUE
	177
	178	END SELECT
	179
	180	END SUBROUTINE vdi_actions
	181	!------------------------------------------------------------------------------!
	182	! Description:
	183	! ------------
	184	!> At a control grid point in the interior of the model domain,
	185	!> 2 * deltat waves are not to be generated with increasing simulation time.
	186	!------------------------------------------------------------------------------!
	187	SUBROUTINE vdi_2_deltat_wave
	188
	189	INTEGER(iwp) :: count_wave = 0 !< counts the number of consecutive waves
	190	INTEGER(iwp) :: count_time = 0 !< counter, so that the waves follow each other without gaps
	191	INTEGER(iwp) :: cgp_i = 0 !< x coordinate of the control grid point for testing 2deltat waves
	192	INTEGER(iwp) :: cgp_j = 0 !< y coordinate of the control grid point for testing 2deltat waves
	193	INTEGER(iwp) :: cgp_k = 0 !< z coordinate of the control grid point for testing 2deltat waves
	194
	195	INTEGER(iwp), DIMENSION(4) :: sig_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of u in the last four time steps
	196
	197	REAL(wp) :: random !< random number
	198
	199	!
	200	!-- Defining the control grid point
	201	IF ( internal_count == 1 ) THEN
	202	cgp_i = INT( nxl + ( nxr - nxl ) / 2 )
	203	cgp_j = INT( nys + ( nyn - nys ) / 2 )
	204	cgp_k = INT( nzt / 2 )
	205	!
	206	!-- If the grid point lies in a building, a new point is defined
[4329]	207	DO WHILE ( .NOT. BTEST( wall_flags_static_0(cgp_k,cgp_j,cgp_i), 1 ) )
[4173]	208	CALL RANDOM_NUMBER( random )
	209	cgp_k = cgp_k + FLOOR( ( nzt - cgp_k ) * random ) !< Random number upon cgp_k
	210	!
	211	!-- If there is topography in the entire grid column, a new x coordinate is chosen
	212	IF ( cgp_k >= nzt-1 ) THEN
	213	CALL RANDOM_NUMBER( random )
	214	cgp_i = nxl + FLOOR( ( nxr + 1 - nxl ) * random )
	215	cgp_k = INT( nzt / 2 )
	216	ENDIF
	217	ENDDO
	218	ENDIF
	219
	220	CALL testing_2_deltat_wave( u_p(cgp_k,cgp_j,cgp_i), u(cgp_k,cgp_j,cgp_i), &
	221	sig_arr, count_wave, count_time )
	222
	223	END SUBROUTINE vdi_2_deltat_wave
	224
	225
	226	!------------------------------------------------------------------------------!
	227	! Description:
	228	! ------------
	229	!> In this subroutine a quantity quant is tested for 2 delta t waves.
	230	!> For this, the size must have a wave-shaped course over 4*4 time steps
	231	!> and the amplitude of the wave has to be greater than the change of quant with
	232	!> increasing time.
	233	!------------------------------------------------------------------------------!
	234	SUBROUTINE testing_2_deltat_wave( quant_p_r, quant_r, sig_arr, count_wave, count_time )
	235
	236	INTEGER(iwp), INTENT(INOUT) :: count_wave !< counts the number of consecutive waves
	237	INTEGER(iwp), INTENT(INOUT) :: count_time !< counter, so that the waves follow each other without gaps
	238	INTEGER(iwp), PARAMETER :: number_wave = 10 !< number of consecutive waves that are not allowed
	239
	240	REAL(wp), INTENT(IN) :: quant_p_r !< quantity from the previous time step as a real
	241	REAL(wp), INTENT(IN) :: quant_r !< quantity as a real
	242	REAL(wp) :: quant_rel = 0.0_wp !< rel. change of the quantity to the previous time step
	243
	244	INTEGER(iwp), DIMENSION(4), INTENT(INOUT) :: sig_arr !< indicates an increase(1) or a decrease (0) of
	245	!> quantity quant in the last four time steps
	246
	247
	248	IF ( quant_p_r - quant_r > 0.0 ) THEN
	249	sig_arr(4) = 0
	250	ELSE
	251	sig_arr(4) = 1
	252	ENDIF
	253
	254	quant_rel = ABS( ( quant_p_r - quant_r ) / quant_p_r )
	255
	256	!
	257	!-- With this criterion 2 delta t waves are detected if the amplitude of
	258	!-- the wave is greater than the change of quant with increasing time
	259	IF ( ALL( sig_arr(1:4) == (/ 1, 0, 1, 0 /) ) .AND. quant_rel > 0.01 ) THEN
	260
	261	count_wave = count_wave + 1
	262
	263	IF ( count_wave == number_wave .AND. count_time == 4 ) THEN
	264	message_string = '2 deltat waves are generated '
	265	CALL message( 'vdi_2_deltat_wave', 'PA0669', 2, 2, myid, 6, 0 )
	266	ENDIF
	267
	268	count_time = 0
	269
	270	ELSE
	271
	272	IF ( count_time >= 4 ) THEN
	273	count_wave = 0
	274	ENDIF
	275
	276	ENDIF
	277
	278	sig_arr(1) = sig_arr(2)
	279	sig_arr(2) = sig_arr(3)
	280	sig_arr(3) = sig_arr(4)
	281
	282	count_time = count_time + 1
	283
	284	END SUBROUTINE testing_2_deltat_wave
	285
	286
	287	!------------------------------------------------------------------------------!
	288	! Description:
	289	! ------------
	290	!> In this internal assessment the maxima of standarddifferences of the
	291	!> meteorological variables, computed layer by layer will be checked.
	292	!> The maxima should not to remain at the open edges of the model or
	293	!> travel from there into the interior of the domain with increasing
	294	!> simulation time.
	295	!> @todo try to reduce repeating code.
	296	!------------------------------------------------------------------------------!
	297	SUBROUTINE vdi_standard_differences
	298
	299	INTEGER(iwp) :: position_u_deviation = 0 !< position of the maximum of the standard deviation of u
	300	INTEGER(iwp) :: position_u_deviation_p = 0 !< position of the maximum of the standard deviation of u to the previous time step
	301	INTEGER(iwp) :: position_u_deviation_pp = 0 !< position of the maximum of the standard deviation of u two time steps ago
	302	INTEGER(iwp) :: position_v_deviation = 0 !< position of the maximum of the standard deviation of v
	303	INTEGER(iwp) :: position_v_deviation_p = 0 !< position of the maximum of the standard deviation of v to the previous time step
	304	INTEGER(iwp) :: position_v_deviation_pp = 0 !< position of the maximum of the standard deviation of v two time steps ago
	305	INTEGER(iwp) :: position_w_deviation = 0 !< position of the maximum of the standard deviation of w
	306	INTEGER(iwp) :: position_w_deviation_p = 0 !< position of the maximum of the standard deviation of w to the previous time step
	307	INTEGER(iwp) :: position_w_deviation_pp = 0 !< position of the maximum of the standard deviation of w two time steps ago
	308	INTEGER(iwp) :: position_pt_deviation = 0 !< position of the maximum of the standard deviation of pt
	309	INTEGER(iwp) :: position_pt_deviation_p = 0 !< position of the maximum of the standard deviation of pt to the previous time step
	310	INTEGER(iwp) :: position_pt_deviation_pp = 0 !< position of the maximum of the standard deviation of pt two time steps ago
	311	INTEGER(iwp) :: position_q_deviation = 0 !< position of the maximum of the standard deviation of q
	312	INTEGER(iwp) :: position_q_deviation_p = 0 !< position of the maximum of the standard deviation of q to the previous time step
	313	INTEGER(iwp) :: position_q_deviation_pp = 0 !< position of the maximum of the standard deviation of q two time steps ago
	314
	315	REAL(wp), DIMENSION(nzb:nzt+1) :: u_deviation !< standard deviation of u depending on k
	316	REAL(wp), DIMENSION(nzb:nzt+1) :: v_deviation !< standard deviation of v depending on k
	317	REAL(wp), DIMENSION(nzb:nzt+1) :: w_deviation !< standard deviation of w depending on k
	318	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_deviation !< standard deviation of pt depending on k
	319	REAL(wp), DIMENSION(nzb:nzt+1) :: q_deviation !< standard deviation of q depending on k
	320
	321	!
	322	!-- Calculation of the standard deviation of u
	323	CALL calc_standard_deviation( u, u_deviation, 1 )
	324
	325	!
	326	!-- Determination of the position of the maximum
	327	position_u_deviation = MAXLOC( u_deviation, DIM=1 )
	328
	329	!
	330	!-- Check the position of the maximum of the standard deviation of u
	331	IF ( internal_count > 2 ) THEN
	332	CALL check_position( position_u_deviation, position_u_deviation_p, position_u_deviation_pp )
	333	ENDIF
	334
	335	position_u_deviation_pp = position_u_deviation_p
	336	position_u_deviation_p = position_u_deviation
	337
	338	!
	339	!-- Calculation of the standard deviation of v
	340	CALL calc_standard_deviation( v, v_deviation, 2 )
	341
	342	!
	343	!-- Determination of the position of the maximum
	344	position_v_deviation = MAXLOC( v_deviation, DIM=1 )
	345
	346	!
	347	!-- Check the position of the maximum of the standard deviation of v
	348	IF ( internal_count > 2 ) THEN
	349	CALL check_position( position_v_deviation, position_v_deviation_p, position_v_deviation_pp )
	350	ENDIF
	351
	352	position_v_deviation_pp = position_v_deviation_p
	353	position_v_deviation_p = position_v_deviation
	354
	355	!
	356	!-- Calculation of the standard deviation of w
	357	CALL calc_standard_deviation( w, w_deviation, 3 )
	358
	359	!
	360	!-- Determination of the position of the maximum
	361	position_w_deviation = MAXLOC( w_deviation, DIM=1 )
	362
	363	!
	364	!-- Check the position of the maximum of the standard deviation of w
	365	IF ( internal_count > 2 ) THEN
	366	CALL check_position( position_w_deviation, position_w_deviation_p, position_w_deviation_pp )
	367	ENDIF
	368
	369	position_w_deviation_pp = position_w_deviation_p
	370	position_w_deviation_p = position_w_deviation
	371
	372
	373	!
	374	!-- Calculation of the standard deviation of pt
	375	IF ( .NOT. neutral ) THEN
	376	CALL calc_standard_deviation( pt, pt_deviation, 0 )
	377	!
	378	!-- Determination of the position of the maximum
	379	position_pt_deviation = MAXLOC( pt_deviation, DIM=1 )
	380
	381	!
	382	!-- Check the position of the maximum of the standard deviation of pt
	383	IF ( internal_count > 2 ) THEN
	384	CALL check_position( position_pt_deviation, &
	385	position_pt_deviation_p, &
	386	position_pt_deviation_pp )
	387	ENDIF
	388
	389	position_pt_deviation_pp = position_pt_deviation_p
	390	position_pt_deviation_p = position_pt_deviation
	391
	392	ENDIF
	393
	394	!
	395	!-- Calculation of the standard deviation of q
	396	IF ( humidity ) THEN
	397	CALL calc_standard_deviation( q, q_deviation, 0 )
	398
	399	!
	400	!-- Determination of the position of the maximum
	401	position_q_deviation = MAXLOC( q_deviation, DIM=1 )
	402
	403	!
	404	!-- Check the position of the maximum of the standard deviation of q
	405	IF ( internal_count > 2 ) THEN
	406	CALL check_position( position_q_deviation, &
	407	position_q_deviation_p, &
	408	position_q_deviation_pp )
	409	ENDIF
	410
	411	position_q_deviation_pp = position_q_deviation_p
	412	position_q_deviation_p = position_q_deviation
	413
	414	ENDIF
	415
	416	END SUBROUTINE vdi_standard_differences
	417
	418
	419	!------------------------------------------------------------------------------!
	420	! Description:
	421	! ------------
	422	!> Calculation of the standard deviation
	423	!------------------------------------------------------------------------------!
	424	SUBROUTINE calc_standard_deviation( quant, std_deviation, quant_type )
	425
	426	INTEGER(iwp) :: i !< loop index
	427	INTEGER(iwp) :: j !< loop index
	428	INTEGER(iwp) :: k !< loop index
	429	INTEGER(iwp), INTENT(IN) :: quant_type !< bit position (1 for u, 2 for v, 3 for w and 0 for scalar)
	430
	431	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: count_2d_l !< counter for averaging (local)
	432	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: count_2d !< counter for averaging
	433
	434	REAL(wp) :: flag !< flag indicating atmosphere (1) or wall (0) grid point
	435
	436	REAL(wp), DIMENSION(nzb:nzt+1) :: quant_av_k_l !< Mean of the quantity quant depending on k (local)
	437	REAL(wp), DIMENSION(nzb:nzt+1) :: quant_av_k !< Mean of the quantity quant depending on k
	438	REAL(wp), DIMENSION(nzb:nzt+1), INTENT(OUT) :: std_deviation !< standard deviation of quant
	439	REAL(wp), DIMENSION(nzb:nzt+1) :: std_deviation_l !< standard deviation of quant (local)
	440
	441	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: quant !< Quantity
	442
	443	!
	444	!-- Calculation of the standard deviation
	445	quant_av_k_l = 0.0_wp
	446	quant_av_k = 0.0_wp
	447	std_deviation = 0.0_wp
	448	std_deviation_l = 0.0_wp
	449	!
	450	!-- Average
	451	count_2d_l = 0
	452	count_2d = 0
	453	DO i = nxl, nxr
	454	DO j = nys, nyn
	455	DO k = nzb+1, nzt+1
[4329]	456	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), quant_type ) )
[4173]	457	quant_av_k_l(k) = quant_av_k_l(k) + quant(k,j,i) * flag
	458	count_2d_l(k) = count_2d_l(k) + INT( flag, KIND=iwp )
	459	ENDDO
	460	ENDDO
	461	ENDDO
	462
	463	#if defined( __parallel )
	464	CALL MPI_ALLREDUCE( quant_av_k_l, quant_av_k, nzt+1-nzb+1, &
	465	MPI_REAL, MPI_SUM, comm2d, ierr )
	466
	467	CALL MPI_ALLREDUCE( count_2d_l, count_2d, nzt+1-nzb+1, &
	468	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	469	#else
	470	quant_av_k = quant_av_k_l
	471	count_2d = count_2d_l
	472	#endif
	473
	474	DO k = nzb+1, nzt+1
	475	quant_av_k(k) = quant_av_k(k) / REAL( count_2d(k), KIND=wp )
	476	ENDDO
	477
	478	DO i = nxl, nxr
	479	DO j = nys, nyn
	480	DO k = nzb+1, nzt+1
	481	std_deviation_l(k) = std_deviation_l(k) &
	482	+ ( quant(k,j,i) - quant_av_k(k) )**2 &
	483	* MERGE( 1.0_wp, 0.0_wp, &
[4329]	484	BTEST( wall_flags_static_0(k,j,i), quant_type ) )
[4173]	485	ENDDO
	486	ENDDO
	487	ENDDO
	488
	489
	490	#if defined( __parallel )
	491	CALL MPI_ALLREDUCE( std_deviation_l, std_deviation, nzt+1-nzb+1, &
	492	MPI_REAL, MPI_SUM, comm2d, ierr )
	493	#else
	494	std_deviation = std_deviation_l
	495	#endif
	496
	497	DO k = nzb+1, nzt+1
	498	std_deviation(k) = SQRT( std_deviation(k) / REAL( count_2d(k), KIND=wp ) )
	499	ENDDO
	500
	501	END SUBROUTINE calc_standard_deviation
	502
	503
	504	!------------------------------------------------------------------------------!
	505	! Description:
	506	! ------------
	507	!> Tests for the position of the maxima of the standard deviation.
	508	!> If the maxima remain at the open edges of the model or travel from
	509	!> the open edges into the interior of the domain with increasing
	510	!> simulation time, the simulation should be aborted.
	511	!------------------------------------------------------------------------------!
	512	SUBROUTINE check_position( position_std_deviation, position_std_deviation_p, &
	513	position_std_deviation_pp )
	514
	515	INTEGER(iwp), INTENT(IN) :: position_std_deviation !< position of the maximum of the std
	516	INTEGER(iwp), INTENT(IN) :: position_std_deviation_p !< previous position of std-max
	517	INTEGER(iwp), INTENT(IN) :: position_std_deviation_pp !< prev. prev. position of std-max
	518
	519
	520	IF ( position_std_deviation == nzt .AND. &
	521	position_std_deviation_p == nzt .AND. &
	522	position_std_deviation_pp == nzt ) THEN
	523	message_string = 'The maxima of the standard deviation remain ' // &
	524	'at the open edges of the model.'
	525	CALL message( 'vdi_standard_differences', 'PA0663', 1, 2, 0, 6, 0 )
	526	ENDIF
	527
	528	IF ( position_std_deviation == nzt-2 .AND. &
	529	position_std_deviation_p == nzt-1 .AND. &
	530	position_std_deviation_pp == nzt ) THEN
	531	message_string = 'The maxima of the standard deviation travel ' // &
	532	'from the open edges into the interior ' // &
	533	'of the domain with increasing simulation time.'
	534	CALL message( 'vdi_standard_differences', 'PA0664', 1, 2, 0, 6, 0 )
	535	ENDIF
	536
	537	END SUBROUTINE check_position
	538
	539
	540	!------------------------------------------------------------------------------!
	541	! Description:
	542	! ------------
	543	!> In this control it will be checked, if the means of the meteorological
	544	!> variables over the model grid are not to exhibit 2 deltat waves or
	545	!> monotonic increase or decrease with increasing simulation time.
	546	!------------------------------------------------------------------------------!
	547	SUBROUTINE vdi_domain_averages
	548
	549	INTEGER(iwp) :: mono_count_u = 0 !< counter for monotonic decrease or increase of u
	550	INTEGER(iwp) :: mono_count_v = 0 !< counter for monotonic decrease or increase of v
	551	INTEGER(iwp) :: mono_count_w = 0 !< counter for monotonic decrease or increase of w
	552	INTEGER(iwp) :: mono_count_q = 0 !< counter for monotonic decrease or increase of q
	553	INTEGER(iwp) :: mono_count_pt = 0 !< counter for monotonic decrease or increase of pt
	554	INTEGER(iwp) :: count_time_u = 0 !< counter, so that the waves of u follow each other without gaps
	555	INTEGER(iwp) :: count_time_v = 0 !< counter, so that the waves of v follow each other without gaps
	556	INTEGER(iwp) :: count_time_w = 0 !< counter, so that the waves of w follow each other without gaps
	557	INTEGER(iwp) :: count_time_q = 0 !< counter, so that the waves of q follow each other without gaps
	558	INTEGER(iwp) :: count_time_pt = 0 !< counter, so that the waves of pt follow each other without gaps
	559	INTEGER(iwp) :: count_wave_u = 0 !< counts the number of consecutive waves of u
	560	INTEGER(iwp) :: count_wave_v = 0 !< counts the number of consecutive waves of v
	561	INTEGER(iwp) :: count_wave_w = 0 !< counts the number of consecutive waves of w
	562	INTEGER(iwp) :: count_wave_q = 0 !< counts the number of consecutive waves of q
	563	INTEGER(iwp) :: count_wave_pt = 0 !< counts the number of consecutive waves of pt
	564
	565	INTEGER(iwp), DIMENSION(4) :: sig_u_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of u in the last four time steps
	566	INTEGER(iwp), DIMENSION(4) :: sig_v_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of v in the last four time steps
	567	INTEGER(iwp), DIMENSION(4) :: sig_w_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of w in the last four time steps
	568	INTEGER(iwp), DIMENSION(4) :: sig_q_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of q in the last four time steps
	569	INTEGER(iwp), DIMENSION(4) :: sig_pt_arr = (/ 0, 0, 0, 0/) !< indicates an increase(1) or a decrease (0) of pt in the last four time steps
	570
	571	REAL(wp) :: u_av = 0.0_wp !< Mean of u
	572	REAL(wp) :: u_av_p = 0.0_wp !< Mean of u at the previous time step
	573	REAL(wp) :: v_av = 0.0_wp !< Mean of v
	574	REAL(wp) :: v_av_p = 0.0_wp !< Mean of v at the previous time step
	575	REAL(wp) :: w_av = 0.0_wp !< Mean of w
	576	REAL(wp) :: w_av_p = 0.0_wp !< Mean of w at the previous time step
	577	REAL(wp) :: q_av = 0.0_wp !< Mean of q
	578	REAL(wp) :: q_av_p = 0.0_wp !< Mean of q at the previous time step
	579	REAL(wp) :: pt_av = 0.0_wp !< Mean of pt
	580	REAL(wp) :: pt_av_p = 0.0_wp !< Mean of pt at the previous time step
	581
	582	!
	583	!-- Averaging the meteorological variables over the model grid
	584	CALL calc_average( u, u_av, 1 )
	585	CALL calc_average( v, v_av, 2 )
	586	CALL calc_average( w, w_av, 3 )
	587	IF ( .NOT. neutral ) THEN
	588	CALL calc_average( pt, pt_av, 0 )
	589	ENDIF
	590	IF ( humidity ) THEN
	591	CALL calc_average( q, q_av, 0 )
	592	ENDIF
	593
	594	!
	595	!-- Testing the meteorological variables for 2 delta t waves
	596	IF ( internal_count > 1 ) THEN
	597	CALL testing_2_deltat_wave( u_av_p, u_av, sig_u_arr, count_wave_u, count_time_u )
	598	CALL testing_2_deltat_wave( v_av_p, v_av, sig_v_arr, count_wave_v, count_time_v )
	599	CALL testing_2_deltat_wave( w_av_p, w_av, sig_w_arr, count_wave_w, count_time_w )
	600	IF ( .NOT. neutral ) THEN
	601	CALL testing_2_deltat_wave( pt_av_p, pt_av, sig_pt_arr, count_wave_pt, count_time_pt )
	602	ENDIF
	603	IF ( humidity ) THEN
	604	CALL testing_2_deltat_wave( q_av_p, q_av, sig_q_arr, count_wave_q, count_time_q )
	605	ENDIF
	606	ENDIF
	607
	608	!
	609	!-- Testing if there is a monotonic increase or decrease with increasing simulation time
	610	IF ( sig_u_arr(2) /= sig_u_arr(3) ) THEN
	611	mono_count_u = 0
	612	ELSE
	613	mono_count_u = mono_count_u + 1
	614	ENDIF
	615
	616	IF ( time_since_reference_point >= end_time .AND. &
	617	mono_count_u > 0.9_wp * internal_count ) THEN
	618
	619	message_string = 'Monotonic decrease or increase with ' // &
	620	'increasing simulation time for u'
	621	CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
	622	ENDIF
	623
	624	IF ( sig_v_arr(2) /= sig_v_arr(3) ) THEN
	625	mono_count_v = 0
	626	ELSE
	627	mono_count_v = mono_count_v + 1
	628	ENDIF
	629
	630	IF ( time_since_reference_point >= end_time .AND. &
	631	mono_count_v > 0.9_wp * internal_count ) THEN
	632	message_string = 'Monotonic decrease or increase with ' // &
	633	'increasing simulation time for v'
	634	CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
	635	ENDIF
	636
	637	IF ( sig_w_arr(2) /= sig_w_arr(3) ) THEN
	638	mono_count_w = 0
	639	ELSE
	640	mono_count_w = mono_count_w + 1
	641	ENDIF
	642
	643	IF ( time_since_reference_point >= end_time .AND. &
	644	mono_count_w > 0.9_wp * internal_count ) THEN
	645	message_string = 'Monotonic decrease or increase with ' // &
	646	'increasing simulation time for w'
	647	CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
	648	ENDIF
	649
	650	IF ( .NOT. neutral ) THEN
	651	IF ( sig_pt_arr(2) /= sig_pt_arr(3) ) THEN
	652	mono_count_pt = 0
	653	ELSE
	654	mono_count_pt = mono_count_pt + 1
	655	ENDIF
	656
	657	IF ( time_since_reference_point >= end_time .AND. &
	658	mono_count_pt > 0.9_wp * internal_count ) THEN
	659	message_string = 'Monotonic decrease or increase with ' // &
	660	'increasing simulation time for pt'
	661	CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
	662	ENDIF
	663	ENDIF
	664
	665	IF ( humidity ) THEN
	666	IF ( sig_q_arr(2) /= sig_q_arr(3) ) THEN
	667	mono_count_q = 0
	668	ELSE
	669	mono_count_q = mono_count_q + 1
	670	ENDIF
	671
	672	IF ( time_since_reference_point >= end_time .AND. &
	673	mono_count_q > 0.9_wp * internal_count ) THEN
	674	message_string = 'Monotonic decrease or increase with ' // &
	675	'increasing simulation time for q'
	676	CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
	677	ENDIF
	678	ENDIF
	679
	680	!
	681	!-- Save the values from the previous time step
	682	u_av_p = u_av
	683	v_av_p = v_av
	684	w_av_p = w_av
	685
	686	IF ( .NOT. neutral ) THEN
	687	pt_av_p = pt_av
	688	ENDIF
	689
	690	IF ( humidity ) THEN
	691	q_av_p = q_av
	692	ENDIF
	693
	694	END SUBROUTINE vdi_domain_averages
	695
	696
	697	!------------------------------------------------------------------------------!
	698	! Description:
	699	! ------------
	700	!> Calculate the average of a quantity 'quant'.
	701	!------------------------------------------------------------------------------!
	702	SUBROUTINE calc_average( quant, quant_av, quant_type )
	703
	704	INTEGER(iwp) :: average_count = 0 !< counter for averaging
	705	INTEGER(iwp) :: average_count_l = 0 !< counter for averaging (local)
	706	INTEGER :: i !< loop index
	707	INTEGER :: j !< loop index
	708	INTEGER :: k !< loop index
	709	INTEGER(iwp) :: quant_type !< bit position (1 for u, 2 for v, 3 for w and 0 for scalar)
	710
	711	REAL(wp) :: flag !< flag indicating atmosphere (1) or wall (0) grid point
	712	REAL(wp) :: quant_av !< average of the quantity quant
	713	REAL(wp) :: quant_av_l = 0.0_wp !< average of the quantity quant (local)
	714
	715	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: quant
	716
	717	!
	718	!-- Averaging the quantity over the model grid
	719	average_count_l = 0
	720	quant_av_l = 0.0_wp
	721	DO i = nxl, nxr
	722	DO j = nys, nyn
	723	DO k = nzb, nzt+1
[4329]	724	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), quant_type ) )
[4173]	725	quant_av_l = quant_av_l + quant(k,j,i) * flag
	726	average_count_l = average_count_l + INT( flag, KIND=iwp )
	727	ENDDO
	728	ENDDO
	729	ENDDO
	730
	731	#if defined( __parallel )
	732	CALL MPI_ALLREDUCE( quant_av_l, quant_av, 1, &
	733	MPI_REAL, MPI_SUM, comm2d, ierr )
	734	CALL MPI_ALLREDUCE( average_count_l, average_count, 1, &
	735	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	736	#else
	737	quant_av = quant_av_l
	738	average_count = average_count_l
	739	#endif
	740
	741	quant_av = quant_av / REAL( average_count, KIND(wp) )
	742
	743	END SUBROUTINE calc_average
	744
	745
	746	!------------------------------------------------------------------------------!
	747	! Description:
	748	! ------------
	749	!> Testing for conservation of mass.
	750	!------------------------------------------------------------------------------!
	751	SUBROUTINE vdi_conservation_of_mass
	752
	753	INTEGER(iwp) :: i !< loop index
	754	INTEGER(iwp) :: j !< loop index
	755	INTEGER(iwp) :: k !< loop index
	756
	757	REAL(wp) :: sum_mass_flux !< sum of the mass flow
	758
	759	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< volume flow (local)
	760	REAL(wp), DIMENSION(1:3) :: volume_flow !< volume flow
	761
	762
	763	volume_flow = 0.0_wp
	764	volume_flow_l = 0.0_wp
	765
	766	!
	767	!-- Left/right:
	768	!-- Sum up the volume flow through the left boundary
	769	IF ( nxl == 0 ) THEN
	770	i = 0
	771	DO j = nys, nyn
	772	DO k = nzb+1, nzt
	773	volume_flow_l(1) = volume_flow_l(1) &
	774	+ u(k,j,i) * dzw(k) * dy &
	775	* MERGE( 1.0_wp, 0.0_wp, &
[4329]	776	BTEST( wall_flags_static_0(k,j,i), 1 ) &
[4173]	777	)
	778	ENDDO
	779	ENDDO
	780	ENDIF
	781	!
	782	!-- Sum up the volume flow through the right boundary
	783	IF ( nxr == nx ) THEN
	784	i = nx+1
	785	DO j = nys, nyn
	786	DO k = nzb+1, nzt
	787	volume_flow_l(1) = volume_flow_l(1) &
	788	- u(k,j,i) * dzw(k) * dy &
	789	* MERGE( 1.0_wp, 0.0_wp, &
[4329]	790	BTEST( wall_flags_static_0(k,j,i), 1 ) &
[4173]	791	)
	792	ENDDO
	793	ENDDO
	794	ENDIF
	795	!
	796	!-- South/north:
	797	!-- Sum up the volume flow through the south boundary
	798	IF ( nys == 0 ) THEN
	799	j = 0
	800	DO i = nxl, nxr
	801	DO k = nzb+1, nzt
	802	volume_flow_l(2) = volume_flow_l(2) &
	803	+ v(k,j,i) * dzw(k) * dx &
	804	* MERGE( 1.0_wp, 0.0_wp, &
[4329]	805	BTEST( wall_flags_static_0(k,j,i), 2 ) &
[4173]	806	)
	807	ENDDO
	808	ENDDO
	809	ENDIF
	810	!
	811	!-- Sum up the volume flow through the north boundary
	812	IF ( nyn == ny ) THEN
	813	j = ny+1
	814	DO i = nxl, nxr
	815	DO k = nzb+1, nzt
	816	volume_flow_l(2) = volume_flow_l(2) &
	817	- v(k,j,i) * dzw(k) * dx &
	818	* MERGE( 1.0_wp, 0.0_wp, &
[4329]	819	BTEST( wall_flags_static_0(k,j,i), 2 ) &
[4173]	820	)
	821	ENDDO
	822	ENDDO
	823	ENDIF
	824	!
	825	!-- Top boundary
	826	k = nzt
	827	DO i = nxl, nxr
	828	DO j = nys, nyn
	829	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
	830	ENDDO
	831	ENDDO
	832
	833	#if defined( __parallel )
	834	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	835	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, comm2d, ierr )
	836	#else
	837	volume_flow = volume_flow_l
	838	#endif
	839
	840	sum_mass_flux = SUM( volume_flow ) / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
	841
	842	IF ( ABS( sum_mass_flux ) > 0.001 ) THEN
	843	message_string = 'The mass is not conserved. '
	844	CALL message( 'vdi_conservation_of_mass', 'PA0666', 1, 2, 0, 6, 0 )
	845	ENDIF
	846
	847	END SUBROUTINE vdi_conservation_of_mass
	848
	849
	850	!------------------------------------------------------------------------------!
	851	! Description:
	852	! ------------
	853	!> The results will be checked for exceedance the specified limits.
	854	!> The controls are performed at every time step and at every grid point.
	855	!> No wind component is allowed to have a magnitude greater than ten times
	856	!> the maximum wind velocity at the approach flow profile (Vdi 3783 part 9).
	857	!> Note, that the supersaturation can not be higher than 10%. Therefore, no
	858	!> test is required.
	859	!------------------------------------------------------------------------------!
	860	SUBROUTINE vdi_plausible_values
	861
	862	INTEGER(iwp) :: i !< loop index
	863	INTEGER(iwp) :: j !< loop index
	864	INTEGER(iwp) :: k !< loop index
	865
	866	REAL(wp) :: max_uv_l_l !< maximum speed at the left edge (local)
	867	REAL(wp) :: max_uv_l !< maximum speed at the left edge
	868	REAL(wp) :: max_uv_r_l !< maximum speed at the right edge (local)
	869	REAL(wp) :: max_uv_r !< maximum speed at the right edge
	870	REAL(wp) :: max_uv_s_l !< maximum speed at the south edge (local)
	871	REAL(wp) :: max_uv_s !< maximum speed at the south edge
	872	REAL(wp) :: max_uv_n_l !< maximum speed at the north edge (local)
	873	REAL(wp) :: max_uv_n !< maximum speed at the north edge
	874	REAL(wp) :: max_uv !< maximum speed of all edges
	875
	876	REAL(wp), DIMENSION(4) :: max_arr !<
	877	REAL(wp), DIMENSION(:), ALLOCATABLE :: uv !< wind velocity at the approach flow
	878	REAL(wp), DIMENSION(:), ALLOCATABLE :: uv_l !< wind velocity at the approach flow (local)
	879
	880	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn) :: uv_l_nest !< wind profile at the left edge (nesting)
	881	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn) :: uv_r_nest !< wind profile at the right edge (nesting)
	882	REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr) :: uv_s_nest !< wind profile at the south edge (nesting)
	883	REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr) :: uv_n_nest !< wind profile at the north edge (nesting)
	884
	885
	886	IF ( .NOT. ALLOCATED( uv ) ) THEN
	887	ALLOCATE( uv(nzb:nzt+1) )
	888	ALLOCATE( uv_l(nzb:nzt+1) )
	889
	890	uv = 0.0_wp
	891	uv_l = 0.0_wp
	892	ENDIF
	893
	894	!
	895	!-- Determination of the approach flow profile
	896	IF ( nested_run ) THEN
	897
	898	uv_l_nest = 0.0_wp
	899	uv_r_nest = 0.0_wp
	900	uv_s_nest = 0.0_wp
	901	uv_n_nest = 0.0_wp
	902	!
	903	!-- Left boundary
	904	IF ( nxl == 0 ) THEN
	905	i = nxl
	906	DO j = nys, nyn
	907	DO k = nzb, nzt+1
	908	uv_l_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2 &
	909	+ ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
	910	ENDDO
	911	ENDDO
	912	max_uv_l_l = MAXVAL(uv_l_nest)
	913	ENDIF
	914	!
	915	!-- Right boundary
	916	IF( nxr == nx ) THEN
	917	i = nxr
	918	DO j = nys, nyn
	919	DO k = nzb, nzt+1
	920	uv_r_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2 &
	921	+ ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
	922
	923	ENDDO
	924	ENDDO
	925	max_uv_r_l = MAXVAL(uv_r_nest)
	926	ENDIF
	927	!
	928	!-- South boundary
	929	IF ( nys == 0 ) THEN
	930	j = nys
	931	DO i = nxl, nxr
	932	DO k = nzb, nzt+1
	933	uv_s_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2 &
	934	+ ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
	935	ENDDO
	936	ENDDO
	937	max_uv_s_l = MAXVAL(uv_s_nest)
	938	ENDIF
	939	!
	940	!-- North boundary
	941	IF ( nyn == ny ) THEN
	942	j = nyn
	943	DO i = nxl, nxr
	944	DO k = nzb, nzt+1
	945	uv_n_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2 &
	946	+ ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
	947
	948	ENDDO
	949	ENDDO
	950	max_uv_n_l = MAXVAL(uv_n_nest)
	951	ENDIF
	952
	953	#if defined( __parallel )
	954	CALL MPI_ALLREDUCE( max_uv_l_l, max_uv_l, 1, MPI_REAL, MPI_MAX, comm2d, ierr )
	955	CALL MPI_ALLREDUCE( max_uv_r_l, max_uv_r, 1, MPI_REAL, MPI_MAX, comm2d, ierr )
	956	CALL MPI_ALLREDUCE( max_uv_s_l, max_uv_s, 1, MPI_REAL, MPI_MAX, comm2d, ierr )
	957	CALL MPI_ALLREDUCE( max_uv_n_l, max_uv_n, 1, MPI_REAL, MPI_MAX, comm2d, ierr )
	958	#else
	959	max_uv_l = max_uv_l_l
	960	max_uv_r = max_uv_r_l
	961	max_uv_s = max_uv_s_l
	962	max_uv_n = max_uv_n_l
	963	#endif
	964
	965	max_arr = (/ max_uv_r, max_uv_l, max_uv_s, max_uv_n /)
	966	max_uv = MAXVAl( max_arr )
	967
	968	ELSE ! non-nested run
	969
	970	IF ( bc_lr_cyc .AND. bc_ns_cyc ) THEN
	971	IF ( nxl == 0 .AND. nys == 0 ) THEN
	972	DO k = nzb, nzt+1
	973	uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2 &
	974	+ ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
	975	ENDDO
	976	ENDIF
	977	ENDIF
	978
	979
	980	IF ( bc_dirichlet_l ) THEN
	981	IF ( nxl == 0 .AND. nys == 0 ) THEN
	982	DO k = nzb, nzt+1
	983	uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2 &
	984	+ ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
	985	ENDDO
	986	ENDIF
	987
	988	ELSEIF (bc_dirichlet_r ) THEN
	989	IF ( nxr == nx .AND. nys == 0 ) THEN
	990	DO k = nzb, nzt+1
	991	uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,nxr) + u(k,0,nxr+1) ) )**2 &
	992	+ ( 0.5_wp * ( v(k,-1,nxr) + v(k,0,nxr) ) )**2 )
	993	ENDDO
	994	ENDIF
	995	ENDIF
	996
	997	IF ( bc_dirichlet_n ) THEN
	998	IF ( nxl == 0 .AND. nyn == ny ) THEN
	999	DO k = nzb, nzt+1
	1000	uv_l(k) = SQRT( ( 0.5_wp * ( u(k,nyn,-1) + u(k,nyn,0) ) )**2 &
	1001	+ ( 0.5_wp * ( v(k,nyn+1,0) + v(k,nyn,0) ) )**2 )
	1002	ENDDO
	1003	ENDIF
	1004
	1005	ELSEIF ( bc_dirichlet_s ) THEN
	1006	IF ( nxl == 0 .AND. nys == 0 ) THEN
	1007	DO k = nzb, nzt+1
	1008	uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2 &
	1009	+ ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
	1010	ENDDO
	1011	ENDIF
	1012	ENDIF
	1013
	1014	#if defined( __parallel )
	1015	CALL MPI_ALLREDUCE( uv_l, uv, nzt+1-nzb+1, MPI_REAL, MPI_MAX, comm2d, ierr )
	1016	#else
	1017	uv = uv_l
	1018	#endif
	1019
	1020	max_uv = MAXVAL( uv )
	1021
	1022	ENDIF
	1023
	1024	!
	1025	!-- Test for exceedance the specified limits
	1026	message_string = 'A wind component have a magnitude greater ' // &
	1027	'than ten times the maximum wind velocity ' // &
	1028	'at the approach flow profile.'
	1029
	1030	IF ( MAXVAL( ABS( u ) ) > 10.0_wp * max_uv ) THEN
	1031	CALL message( 'vdi_plausible_values', 'PA0667', 2, 2, myid, 6, 0 )
	1032	ENDIF
	1033
	1034	IF ( MAXVAL( ABS( v ) ) > 10.0_wp * max_uv ) THEN
	1035	CALL message( 'vdi_plausible_values', 'PA0667', 2, 2, myid, 6, 0 )
	1036	ENDIF
	1037
	1038	IF ( MAXVAL( ABS( w ) ) > 10.0_wp * max_uv ) THEN
	1039	CALL message( 'vdi_plausible_values', 'PA0667', 2, 2, myid, 6, 0 )
	1040	ENDIF
	1041
	1042	!
	1043	!-- Test if the potential temperature lies between 220 K and 330 K
	1044	IF ( MAXVAL( pt ) > 330.0_wp .OR. MAXVAL( pt ) < 220.0_wp ) THEN
	1045	message_string = 'The potential temperature does not lie ' // &
	1046	'between 220 K and 330 K.'
	1047	CALL message( 'vdi_plausible_values', 'PA0668', 2, 2, myid, 6, 0 )
	1048	ENDIF
	1049
	1050	END SUBROUTINE vdi_plausible_values
	1051
	1052	END MODULE vdi_internal_controls

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |