Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |