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

Last change on this file since 4466 was 4415, checked in by raasch, 4 years ago

bugfix for misplaced/missing preprocessor directives

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