Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |