Home

Context Navigation

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

Last change on this file since 4523 was 4497, checked in by raasch, 5 years ago
last bugfix deactivated because of compile problems, files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 43.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |