Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |