Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |