Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |