Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |