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

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