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