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

Last change on this file since 4274 was 4233, checked in by knoop, 5 years ago

OpenACC support added to global_min_max.f90

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