Home

Context Navigation

vdi_internal_controls.f90 @ 4889

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

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

Context Navigation

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

Download in other formats: