source: palm/trunk/SOURCE/global_min_max.f90 @ 4725

Last change on this file since 4725 was 4651, checked in by raasch, 4 years ago

preprocessor branch for ibm removed

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