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

Last change on this file since 4329 was 4329, checked in by motisi, 20 months ago

Renamed wall_flags_0 to wall_flags_static_0

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