Home

Context Navigation

vdi_internal_controls.f90 @ 4883

Last change on this file since 4883 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 43.2 KB

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

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

Context Navigation

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

Download in other formats: