Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |