Home

Context Navigation

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

Last change on this file since 4222 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 39.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |