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

Last change on this file since 4739 was 4557, checked in by raasch, 4 years ago

bugfix: MPI_DOUBLE_PRECISION replaced by MPI_REAL

  • Property svn:keywords set to Id
File size: 43.2 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 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 maronga $
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
Note: See TracBrowser for help on using the repository browser.