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

Last change on this file since 4395 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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