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

Last change on this file since 4178 was 4175, checked in by gronemeier, 2 years ago

bugfix: removed unused variables

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