source: palm/trunk/SOURCE/global_min_max.f90

Last change on this file was 4855, checked in by raasch, 3 years ago

bugfix: mean w removal not applied to ghost points of the total domain in case of non-cyclic setups (pres), bugfix for correct identification of indices of extreme values in case of non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 12.6 KB
Line 
1!--------------------------------------------------------------------------------------------------!
2! This file is part of the PALM model system.
3!
4! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
5! Public License as published by the Free Software Foundation, either version 3 of the License, or
6! (at your option) any later version.
7!
8! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
9! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
10! Public License for more details.
11!
12! You should have received a copy of the GNU General Public License along with PALM. If not, see
13! <http://www.gnu.org/licenses/>.
14!
15! Copyright 1997-2021 Leibniz Universitaet Hannover
16!--------------------------------------------------------------------------------------------------!
17!
18! Current revisions:
19! ------------------
20!
21!
22! Former revisions:
23! -----------------
24! $Id: global_min_max.f90 4855 2021-01-25 12:30:54Z banzhafs $
25! bugfix for correct identification of indices of extreme values in case of non-cyclic
26! boundary conditions
27!
28! 4828 2021-01-05 11:21:41Z Giersch
29! preprocessor branch for ibm removed
30!
31! 4646 2020-08-24 16:02:40Z raasch
32! file re-formatted to follow the PALM coding standard
33!
34! 4429 2020-02-27 15:24:30Z raasch
35! bugfix: cpp-directives added for serial mode
36!
37! 4360 2020-01-07 11:25:50Z suehring
38! OpenACC support added
39!
40! 4182 2019-08-22 15:20:23Z scharf
41! Corrected "Former revisions" section
42!
43! 3655 2019-01-07 16:51:22Z knoop
44! Corrected "Former revisions" section
45!
46! Revision 1.1  1997/07/24 11:14:03  raasch
47! Initial revision
48!
49!
50! Description:
51! ------------
52!> Determine the array minimum/maximum and the corresponding indices.
53!--------------------------------------------------------------------------------------------------!
54 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, value_ijk, value1,    &
55                            value1_ijk )
56
57
58    USE control_parameters,                                                                        &
59        ONLY:  bc_lr, bc_ns
60
61    USE indices,                                                                                   &
62        ONLY:  nbgp, ny, nx
63
64    USE kinds
65
66    USE pegrid
67
68    IMPLICIT NONE
69
70    CHARACTER (LEN=*) ::  mode  !<
71
72    INTEGER(iwp) ::  i              !<
73    INTEGER(iwp) ::  i1             !<
74    INTEGER(iwp) ::  i2             !<
75#if defined( __parallel )
76    INTEGER(iwp) ::  id_fmax        !<
77    INTEGER(iwp) ::  id_fmin        !<
78#endif
79    INTEGER(iwp) ::  j              !<
80    INTEGER(iwp) ::  j1             !<
81    INTEGER(iwp) ::  j2             !<
82    INTEGER(iwp) ::  k              !<
83    INTEGER(iwp) ::  k1             !<
84    INTEGER(iwp) ::  k2             !<
85    INTEGER(iwp) ::  value_ijk(3)   !<
86
87    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk    !<
88    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk_l  !<
89    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk    !<
90    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk_l  !<
91
92    INTEGER(iwp), DIMENSION(3), OPTIONAL ::  value1_ijk  !<
93
94    REAL(wp) ::  offset            !<
95    REAL(wp) ::  value             !<
96    REAL(wp), OPTIONAL ::  value1  !<
97
98    REAL(wp), DIMENSION(2) ::  fmax    !<
99    REAL(wp), DIMENSION(2) ::  fmax_l  !<
100    REAL(wp), DIMENSION(2) ::  fmin    !<
101    REAL(wp), DIMENSION(2) ::  fmin_l  !<
102
103    REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::  ar  !<
104
105#if defined( _OPENACC )
106    INTEGER(iwp) ::  count_eq   !< counter for locations of maximum
107    REAL(wp)     ::  red        !< scalar for reduction with OpenACC
108#endif
109
110
111!
112!-- Determine array minimum
113    IF ( mode == 'min'  .OR.  mode == 'minmax' )  THEN
114
115!
116!--    Determine the local minimum
117       fmin_ijk_l = MINLOC( ar )
118       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1 ! MINLOC assumes lowerbound = 1
119       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
120       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - nbgp
121       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
122
123#if defined( __parallel )
124       fmin_l(2)  = myid
125       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
126       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, ierr )
127
128!
129!--    Determine the global minimum. Result stored on PE0.
130       id_fmin = fmin(2)
131       IF ( id_fmin /= 0 )  THEN
132          IF ( myid == 0 )  THEN
133             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, status, ierr )
134          ELSEIF ( myid == id_fmin )  THEN
135             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
136          ENDIF
137       ELSE
138          fmin_ijk = fmin_ijk_l
139       ENDIF
140!
141!--    Send the indices of the just determined array minimum to other PEs
142       CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
143#else
144       fmin(1)  = fmin_l(1)
145       fmin_ijk = fmin_ijk_l
146#endif
147
148    ENDIF
149
150!
151!-- Determine array maximum
152    IF ( mode == 'max'  .OR.  mode == 'minmax' )  THEN
153
154!
155!--    Determine the local maximum
156       fmax_ijk_l = MAXLOC( ar )
157       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1 ! MAXLOC assumes lowerbound = 1
158       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
159       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - nbgp
160       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
161
162#if defined( __parallel )
163       fmax_l(2)  = myid
164       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
165       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
166
167!
168!--    Determine the global maximum. Result stored on PE0.
169       id_fmax = fmax(2)
170       IF ( id_fmax /= 0 )  THEN
171          IF ( myid == 0 )  THEN
172             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
173          ELSEIF ( myid == id_fmax )  THEN
174             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
175          ENDIF
176       ELSE
177          fmax_ijk = fmax_ijk_l
178       ENDIF
179!
180!--    Send the indices of the just determined array maximum to other PEs
181       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
182#else
183       fmax(1)  = fmax_l(1)
184       fmax_ijk = fmax_ijk_l
185#endif
186
187    ENDIF
188
189!
190!-- Determine absolute array maximum
191    IF ( mode == 'abs' )  THEN
192
193#if defined( _OPENACC )
194       red = 0.0_wp
195       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
196       !$ACC PRESENT(ar) COPY(red) REDUCTION(MAX: red)
197       DO  k = k1, k2
198          DO  j = j1, j2
199             DO  i = i1, i2
200                IF ( ABS( ar(i,j,k) ) > red )  THEN
201                   red = ABS( ar(i,j,k) )
202                ENDIF
203             ENDDO
204          ENDDO
205       ENDDO
206       fmax_l(1) = red
207
208!
209!--    Determine the maximum's position and count how often it is found.
210       count_eq = 0
211       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
212       !$ACC PRESENT(ar) COPY(fmax_ijk_l(1:3), count_eq) &
213       !$ACC REDUCTION(+:count_eq)
214       DO  k = k1, k2
215          DO  j = j1, j2
216             DO  i = i1, i2
217                IF ( ABS( ar(i,j,k) ) == red )  THEN
218                   fmax_ijk_l(1) = i
219                   fmax_ijk_l(2) = j
220                   fmax_ijk_l(3) = k
221                   count_eq = count_eq + 1
222                ENDIF
223             ENDDO
224          ENDDO
225       ENDDO
226
227       IF ( count_eq == 1 ) THEN
228!
229!--       We found a single maximum element and correctly got its position. Transfer its value to
230!--       handle the negative case correctly.
231          !$ACC UPDATE HOST(ar(fmax_ijk_l(1):fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)))
232       ELSE
233!
234!--       We found no maximum element (?) or multiple, so the position is not correct.
235!--       Copy the full array to the host and determine the maximum sequentially...
236          !$ACC UPDATE HOST(ar(i1:i2,j1:j2,k1:k2))
237#endif
238
239!
240!--    Determine the local absolut maximum
241       fmax_l(1)     = 0.0_wp
242       fmax_ijk_l(1) =  i1
243       fmax_ijk_l(2) =  j1
244       fmax_ijk_l(3) =  k1
245       DO  k = k1, k2
246          DO  j = j1, j2
247             DO  i = i1, i2
248                IF ( ABS( ar(i,j,k) ) > fmax_l(1) )  THEN
249                   fmax_l(1) = ABS( ar(i,j,k) )
250                   fmax_ijk_l(1) = i
251                   fmax_ijk_l(2) = j
252                   fmax_ijk_l(3) = k
253                ENDIF
254             ENDDO
255          ENDDO
256       ENDDO
257
258#if defined( _OPENACC )
259!
260!--    Close ELSE case from above
261       ENDIF
262#endif
263
264!
265!--    Set a flag in case that the determined value is negative.
266!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
267       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
268          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
269       ENDIF
270
271#if defined( __parallel )
272       fmax_l(2)  = myid
273       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
274       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
275
276!
277!--    Determine the global absolut maximum. Result stored on PE0.
278       id_fmax = fmax(2)
279       IF ( id_fmax /= 0 )  THEN
280          IF ( myid == 0 )  THEN
281             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
282          ELSEIF ( myid == id_fmax )  THEN
283             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
284          ENDIF
285       ELSE
286          fmax_ijk = fmax_ijk_l
287       ENDIF
288!
289!--    Send the indices of the just determined absolut maximum to other PEs
290       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
291#else
292       fmax(1)  = fmax_l(1)
293       fmax_ijk = fmax_ijk_l
294#endif
295
296    ENDIF
297
298!
299!-- Determine absolute maximum of ( array - offset )
300    IF ( mode == 'absoff' )  THEN
301
302!
303!--    Determine the local absolut maximum
304       fmax_l(1)     = 0.0_wp
305       fmax_ijk_l(1) =  i1
306       fmax_ijk_l(2) =  j1
307       fmax_ijk_l(3) =  k1
308       DO  k = k1, k2
309          DO  j = j1, j2
310!
311!--          Attention: the lowest gridpoint is excluded here, because there is no advection at
312!--          ---------- nzb=0 and mode 'absoff' is only used for calculating u,v extrema for
313!--                     CFL-criteria.
314             DO  i = i1+1, i2
315                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
316                   fmax_l(1) = ABS( ar(i,j,k) - offset )
317                   fmax_ijk_l(1) = i
318                   fmax_ijk_l(2) = j
319                   fmax_ijk_l(3) = k
320                ENDIF
321             ENDDO
322          ENDDO
323       ENDDO
324
325!
326!--    Set a flag in case that the determined value is negative.
327!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
328       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
329          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
330       ENDIF
331
332#if defined( __parallel )
333       fmax_l(2)  = myid
334       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
335       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
336
337!
338!--    Determine the global absolut maximum. Result stored on PE0.
339       id_fmax = fmax(2)
340       IF ( id_fmax /= 0 )  THEN
341          IF ( myid == 0 )  THEN
342             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
343          ELSEIF ( myid == id_fmax )  THEN
344             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
345          ENDIF
346       ELSE
347          fmax_ijk = fmax_ijk_l
348       ENDIF
349!
350!--    Send the indices of the just determined absolut maximum to other PEs
351       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
352#else
353       fmax(1)  = fmax_l(1)
354       fmax_ijk = fmax_ijk_l
355#endif
356
357    ENDIF
358
359!
360!-- Determine output parameters
361    SELECT CASE( mode )
362
363       CASE( 'min' )
364
365          value     = fmin(1)
366          value_ijk = fmin_ijk
367
368       CASE( 'max' )
369
370          value     = fmax(1)
371          value_ijk = fmax_ijk
372
373       CASE( 'minmax' )
374
375          value      = fmin(1)
376          value_ijk  = fmin_ijk
377          value1     = fmax(1)
378          value1_ijk = fmax_ijk
379
380       CASE( 'abs', 'absoff' )
381
382          value     = fmax(1)
383          value_ijk = fmax_ijk
384          IF ( fmax_ijk(1) <= -10 )  THEN
385!
386!--          Index needs to be corrected because it has been modified above to indicate negative
387!--          values
388             value_ijk(1) = -value_ijk(1) - 10
389!
390!--          For this reason also change the sign of the quantity
391             value        = -value
392          ENDIF
393
394    END SELECT
395
396!
397!-- Limit index values to the range 0..nx, 0..ny. Non-cyclic setups may have extrema at the
398!-- outer borders, which should be correctly identified.
399    IF ( bc_lr == 'cyclic' )  THEN
400       IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
401       IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
402    ENDIF
403    IF ( bc_ns == 'cyclic' )  THEN
404       IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
405       IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
406    ENDIF
407
408 END SUBROUTINE global_min_max
Note: See TracBrowser for help on using the repository browser.