Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |