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

Last change on this file since 4598 was 4429, checked in by raasch, 15 months ago

serial (non-MPI) test case added, several bugfixes for the serial mode

  • Property svn:keywords set to Id
File size: 12.2 KB
RevLine 
[1682]1!> @file global_min_max.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[4360]17! Copyright 1997-2020 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[866]21! ------------------
[1354]22!
[2001]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: global_min_max.f90 4429 2020-02-27 15:24:30Z suehring $
[4429]27! bugfix: cpp-directives added for serial mode
28!
29! 4360 2020-01-07 11:25:50Z suehring
[4233]30! OpenACC support added
31!
32! 4182 2019-08-22 15:20:23Z scharf
[2716]33! Corrected "Former revisions" section
34!
[4182]35! 3655 2019-01-07 16:51:22Z knoop
36! Corrected "Former revisions" section
[1321]37!
[4182]38! Revision 1.1  1997/07/24 11:14:03  raasch
39! Initial revision
40!
41!
[1]42! Description:
43! ------------
[1682]44!> Determine the array minimum/maximum and the corresponding indices.
[623]45!------------------------------------------------------------------------------!
[1682]46 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
47                            value_ijk, value1, value1_ijk )
48 
[1]49
[1320]50    USE indices,                                                               &
51        ONLY:  nbgp, ny, nx
52       
53    USE kinds
54   
[1]55    USE pegrid
56
57    IMPLICIT NONE
58
[1682]59    CHARACTER (LEN=*) ::  mode  !<
[1]60
[1682]61    INTEGER(iwp) ::  i              !<
62    INTEGER(iwp) ::  i1             !<
63    INTEGER(iwp) ::  i2             !<
[4429]64#if defined( __parallel )
[1682]65    INTEGER(iwp) ::  id_fmax        !<
66    INTEGER(iwp) ::  id_fmin        !<
[4429]67#endif
[1682]68    INTEGER(iwp) ::  j              !<
69    INTEGER(iwp) ::  j1             !<
70    INTEGER(iwp) ::  j2             !<
71    INTEGER(iwp) ::  k              !<
72    INTEGER(iwp) ::  k1             !<
73    INTEGER(iwp) ::  k2             !<
74    INTEGER(iwp) ::  fmax_ijk(3)    !<
75    INTEGER(iwp) ::  fmax_ijk_l(3)  !<
76    INTEGER(iwp) ::  fmin_ijk(3)    !<
77    INTEGER(iwp) ::  fmin_ijk_l(3)  !<
78    INTEGER(iwp) ::  value_ijk(3)   !<
[1320]79   
[1682]80    INTEGER(iwp), OPTIONAL ::  value1_ijk(3)  !<
[1320]81   
[1682]82    REAL(wp) ::  offset                 !<
83    REAL(wp) ::  value                  !<
84    REAL(wp) ::  ar(i1:i2,j1:j2,k1:k2)  !<
[1320]85   
[1]86#if defined( __ibm )
[1682]87    REAL(sp) ::  fmax(2)    !<
88    REAL(sp) ::  fmax_l(2)  !<
89    REAL(sp) ::  fmin(2)    !<
90    REAL(sp) ::  fmin_l(2)  !<
[1320]91             ! on 32bit-machines MPI_2REAL must not be replaced
92             ! by MPI_2DOUBLE_PRECISION
[1]93#else
[1682]94    REAL(wp) ::  fmax(2)    !<
95    REAL(wp) ::  fmax_l(2)  !<
96    REAL(wp) ::  fmin(2)    !<
97    REAL(wp) ::  fmin_l(2)  !<
[1]98#endif
[4233]99#if defined( _OPENACC )
100    REAL(wp)     ::  red        !< scalar for reduction with OpenACC
101    INTEGER(iwp) ::  count_eq   !< counter for locations of maximum
102#endif
[1682]103    REAL(wp), OPTIONAL ::  value1  !<
[1]104
105
106!
107!-- Determine array minimum
108    IF ( mode == 'min'  .OR.  mode == 'minmax' )  THEN
109
110!
111!--    Determine the local minimum
112       fmin_ijk_l = MINLOC( ar )
[1188]113       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1 ! MINLOC assumes lowerbound = 1
[667]114       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
[1188]115       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - nbgp
[1]116       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
117
118#if defined( __parallel )
119       fmin_l(2)  = myid
[622]120       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[623]121       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, &
122                           ierr )
[1]123
124!
125!--    Determine the global minimum. Result stored on PE0.
126       id_fmin = fmin(2)
127       IF ( id_fmin /= 0 )  THEN
128          IF ( myid == 0 )  THEN
129             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, &
130                            status, ierr )
131          ELSEIF ( myid == id_fmin )  THEN
132             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
133          ENDIF
134       ELSE
135          fmin_ijk = fmin_ijk_l
136       ENDIF
137!
138!--    Send the indices of the just determined array minimum to other PEs
139       CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
140#else
141       fmin(1)  = fmin_l(1)
142       fmin_ijk = fmin_ijk_l
143#endif
144
145    ENDIF
146
147!
148!-- Determine array maximum
149    IF ( mode == 'max'  .OR.  mode == 'minmax' )  THEN
150
151!
152!--    Determine the local maximum
153       fmax_ijk_l = MAXLOC( ar )
[1188]154       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1 ! MAXLOC assumes lowerbound = 1
[667]155       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
[1188]156       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - nbgp
[1]157       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
158
159#if defined( __parallel )
160       fmax_l(2)  = myid
[622]161       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[623]162       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
163                           ierr )
[1]164
165!
166!--    Determine the global maximum. Result stored on PE0.
167       id_fmax = fmax(2)
168       IF ( id_fmax /= 0 )  THEN
169          IF ( myid == 0 )  THEN
170             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
171                            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
[4233]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
[1]207!
[4233]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
229!--       value to 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!
[1]239!--    Determine the local absolut maximum
[1353]240       fmax_l(1)     = 0.0_wp
[1]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
[4233]257#if defined( _OPENACC )
[1]258!
[4233]259!--       Close ELSE case from above
260       ENDIF
261#endif
262
263!
[1]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
266!--    case i=0 correctly
[1353]267       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
[1]268          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
269       ENDIF
270
271#if defined( __parallel )
272       fmax_l(2)  = myid
[622]273       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]274       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
275                           ierr )
276
277!
278!--    Determine the global absolut maximum. Result stored on PE0.
279       id_fmax = fmax(2)
280       IF ( id_fmax /= 0 )  THEN
281          IF ( myid == 0 )  THEN
282             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
283                            status, ierr )
284          ELSEIF ( myid == id_fmax )  THEN
285             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
286          ENDIF
287       ELSE
288          fmax_ijk = fmax_ijk_l
289       ENDIF
290!
291!--    Send the indices of the just determined absolut maximum to other PEs
292       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
293#else
294       fmax(1)  = fmax_l(1)
295       fmax_ijk = fmax_ijk_l
296#endif
297
298    ENDIF
299
300!
[866]301!-- Determine absolute maximum of ( array - offset )
302    IF ( mode == 'absoff' )  THEN
303
304!
305!--    Determine the local absolut maximum
[1353]306       fmax_l(1)     = 0.0_wp
[866]307       fmax_ijk_l(1) =  i1
308       fmax_ijk_l(2) =  j1
309       fmax_ijk_l(3) =  k1
310       DO  k = k1, k2
311          DO  j = j1, j2
312!
313!--          Attention: the lowest gridpoint is excluded here, because there
314!--          ---------  is no advection at nzb=0 and mode 'absoff' is only
315!--                     used for calculating u,v extrema for CFL-criteria
316             DO  i = i1+1, i2
317                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
318                   fmax_l(1) = ABS( ar(i,j,k) - offset )
319                   fmax_ijk_l(1) = i
320                   fmax_ijk_l(2) = j
321                   fmax_ijk_l(3) = k
322                ENDIF
323             ENDDO
324          ENDDO
325       ENDDO
326
327!
328!--    Set a flag in case that the determined value is negative.
329!--    A constant offset has to be subtracted in order to handle the special
330!--    case i=0 correctly
[1353]331       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
[866]332          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
333       ENDIF
334
335#if defined( __parallel )
336       fmax_l(2)  = myid
337       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
338       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
339                           ierr )
340
341!
342!--    Determine the global absolut maximum. Result stored on PE0.
343       id_fmax = fmax(2)
344       IF ( id_fmax /= 0 )  THEN
345          IF ( myid == 0 )  THEN
346             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
347                            status, ierr )
348          ELSEIF ( myid == id_fmax )  THEN
349             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
350          ENDIF
351       ELSE
352          fmax_ijk = fmax_ijk_l
353       ENDIF
354!
355!--    Send the indices of the just determined absolut maximum to other PEs
356       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
357#else
358       fmax(1)  = fmax_l(1)
359       fmax_ijk = fmax_ijk_l
360#endif
361
362    ENDIF
363
364!
[1]365!-- Determine output parameters
366    SELECT CASE( mode )
367
368       CASE( 'min' )
369
370          value     = fmin(1)
371          value_ijk = fmin_ijk
372
373       CASE( 'max' )
374
375          value     = fmax(1)
376          value_ijk = fmax_ijk
377
378       CASE( 'minmax' )
379
380          value      = fmin(1)
381          value_ijk  = fmin_ijk
382          value1     = fmax(1)
383          value1_ijk = fmax_ijk
384
[866]385       CASE( 'abs', 'absoff' )
[1]386
387          value     = fmax(1)
388          value_ijk = fmax_ijk
389          IF ( fmax_ijk(1) < 0 )  THEN
390             value        = -value
[667]391             value_ijk(1) = -value_ijk(1) - 10         !???
[1]392          ENDIF
393
394    END SELECT
395
396!
397!-- Limit index values to the range 0..nx, 0..ny
[667]398    IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
399    IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
400    IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
401    IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
[1]402
403
404 END SUBROUTINE global_min_max
Note: See TracBrowser for help on using the repository browser.