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

Last change on this file since 4647 was 4646, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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