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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

  • Property svn:keywords set to Id
File size: 9.9 KB
Line 
1 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
2                            value_ijk, value1, value1_ijk )
3
4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2012  Leibniz University Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: global_min_max.f90 1036 2012-10-22 13:43:42Z raasch $
28!
29! 866 2012-03-28 06:44:41Z raasch
30! new mode "absoff" accounts for an offset in the respective array
31!
32! 667 2010-12-23 12:06:00Z suehring/gryschka
33! Adapting of the index arrays, because MINLOC assumes lowerbound at 1 and not
34! at nbgp.
35!
36! 622 2010-12-10 08:08:13Z raasch
37! optional barriers included in order to speed up collective operations
38!
39! Feb. 2007
40! RCS Log replace by Id keyword, revision history cleaned up
41!
42! Revision 1.11  2003/04/16 12:56:58  raasch
43! Index values of the extrema are limited to the range 0..nx, 0..ny
44!
45! Revision 1.1  1997/07/24 11:14:03  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51! Determine the array minimum/maximum and the corresponding indices.
52!------------------------------------------------------------------------------!
53
54    USE indices
55    USE pegrid
56
57    IMPLICIT NONE
58
59    CHARACTER (LEN=*) ::  mode
60
61    INTEGER           ::  i, i1, i2, id_fmax, id_fmin, j, j1, j2, k, k1, k2, &
62                          fmax_ijk(3), fmax_ijk_l(3), fmin_ijk(3), &
63                          fmin_ijk_l(3), value_ijk(3)
64    INTEGER, OPTIONAL ::  value1_ijk(3)
65    REAL              ::  offset, value, &
66                          ar(i1:i2,j1:j2,k1:k2)
67#if defined( __ibm )
68    REAL (KIND=4)     ::  fmax(2), fmax_l(2), fmin(2), fmin_l(2)  ! on 32bit-
69                          ! machines MPI_2REAL must not be replaced by
70                          ! MPI_2DOUBLE_PRECISION
71#else
72    REAL              ::  fmax(2), fmax_l(2), fmin(2), fmin_l(2)
73#endif
74    REAL, OPTIONAL    ::  value1
75
76
77!
78!-- Determine array minimum
79    IF ( mode == 'min'  .OR.  mode == 'minmax' )  THEN
80
81!
82!--    Determine the local minimum
83       fmin_ijk_l = MINLOC( ar )
84       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp ! MINLOC assumes lowerbound = 1
85       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
86       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1
87       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
88
89#if defined( __parallel )
90       fmin_l(2)  = myid
91       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
92       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, &
93                           ierr )
94
95!
96!--    Determine the global minimum. Result stored on PE0.
97       id_fmin = fmin(2)
98       IF ( id_fmin /= 0 )  THEN
99          IF ( myid == 0 )  THEN
100             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, &
101                            status, ierr )
102          ELSEIF ( myid == id_fmin )  THEN
103             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
104          ENDIF
105       ELSE
106          fmin_ijk = fmin_ijk_l
107       ENDIF
108!
109!--    Send the indices of the just determined array minimum to other PEs
110       CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
111#else
112       fmin(1)  = fmin_l(1)
113       fmin_ijk = fmin_ijk_l
114#endif
115
116    ENDIF
117
118!
119!-- Determine array maximum
120    IF ( mode == 'max'  .OR.  mode == 'minmax' )  THEN
121
122!
123!--    Determine the local maximum
124       fmax_ijk_l = MAXLOC( ar )
125       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp ! MAXLOC assumes lowerbound = 1
126       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
127       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1
128       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
129
130#if defined( __parallel )
131       fmax_l(2)  = myid
132       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
133       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
134                           ierr )
135
136!
137!--    Determine the global maximum. Result stored on PE0.
138       id_fmax = fmax(2)
139       IF ( id_fmax /= 0 )  THEN
140          IF ( myid == 0 )  THEN
141             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
142                            status, ierr )
143          ELSEIF ( myid == id_fmax )  THEN
144             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
145          ENDIF
146       ELSE
147          fmax_ijk = fmax_ijk_l
148       ENDIF
149!
150!--    send the indices of the just determined array maximum to other PEs
151       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
152#else
153       fmax(1)  = fmax_l(1)
154       fmax_ijk = fmax_ijk_l
155#endif
156
157    ENDIF
158
159!
160!-- Determine absolute array maximum
161    IF ( mode == 'abs' )  THEN
162
163!
164!--    Determine the local absolut maximum
165       fmax_l(1)     = 0.0
166       fmax_ijk_l(1) =  i1
167       fmax_ijk_l(2) =  j1
168       fmax_ijk_l(3) =  k1
169       DO  k = k1, k2
170          DO  j = j1, j2
171             DO  i = i1, i2
172                IF ( ABS( ar(i,j,k) ) > fmax_l(1) )  THEN
173                   fmax_l(1) = ABS( ar(i,j,k) )
174                   fmax_ijk_l(1) = i
175                   fmax_ijk_l(2) = j
176                   fmax_ijk_l(3) = k
177                ENDIF
178             ENDDO
179          ENDDO
180       ENDDO
181
182!
183!--    Set a flag in case that the determined value is negative.
184!--    A constant offset has to be subtracted in order to handle the special
185!--    case i=0 correctly
186       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0 )  THEN
187          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
188       ENDIF
189
190#if defined( __parallel )
191       fmax_l(2)  = myid
192       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
193       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
194                           ierr )
195
196!
197!--    Determine the global absolut maximum. Result stored on PE0.
198       id_fmax = fmax(2)
199       IF ( id_fmax /= 0 )  THEN
200          IF ( myid == 0 )  THEN
201             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
202                            status, ierr )
203          ELSEIF ( myid == id_fmax )  THEN
204             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
205          ENDIF
206       ELSE
207          fmax_ijk = fmax_ijk_l
208       ENDIF
209!
210!--    Send the indices of the just determined absolut maximum to other PEs
211       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
212#else
213       fmax(1)  = fmax_l(1)
214       fmax_ijk = fmax_ijk_l
215#endif
216
217    ENDIF
218
219!
220!-- Determine absolute maximum of ( array - offset )
221    IF ( mode == 'absoff' )  THEN
222
223!
224!--    Determine the local absolut maximum
225       fmax_l(1)     = 0.0
226       fmax_ijk_l(1) =  i1
227       fmax_ijk_l(2) =  j1
228       fmax_ijk_l(3) =  k1
229       DO  k = k1, k2
230          DO  j = j1, j2
231!
232!--          Attention: the lowest gridpoint is excluded here, because there
233!--          ---------  is no advection at nzb=0 and mode 'absoff' is only
234!--                     used for calculating u,v extrema for CFL-criteria
235             DO  i = i1+1, i2
236                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
237                   fmax_l(1) = ABS( ar(i,j,k) - offset )
238                   fmax_ijk_l(1) = i
239                   fmax_ijk_l(2) = j
240                   fmax_ijk_l(3) = k
241                ENDIF
242             ENDDO
243          ENDDO
244       ENDDO
245
246!
247!--    Set a flag in case that the determined value is negative.
248!--    A constant offset has to be subtracted in order to handle the special
249!--    case i=0 correctly
250       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0 )  THEN
251          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
252       ENDIF
253
254#if defined( __parallel )
255       fmax_l(2)  = myid
256       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
257       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
258                           ierr )
259
260!
261!--    Determine the global absolut maximum. Result stored on PE0.
262       id_fmax = fmax(2)
263       IF ( id_fmax /= 0 )  THEN
264          IF ( myid == 0 )  THEN
265             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
266                            status, ierr )
267          ELSEIF ( myid == id_fmax )  THEN
268             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
269          ENDIF
270       ELSE
271          fmax_ijk = fmax_ijk_l
272       ENDIF
273!
274!--    Send the indices of the just determined absolut maximum to other PEs
275       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
276#else
277       fmax(1)  = fmax_l(1)
278       fmax_ijk = fmax_ijk_l
279#endif
280
281    ENDIF
282
283!
284!-- Determine output parameters
285    SELECT CASE( mode )
286
287       CASE( 'min' )
288
289          value     = fmin(1)
290          value_ijk = fmin_ijk
291
292       CASE( 'max' )
293
294          value     = fmax(1)
295          value_ijk = fmax_ijk
296
297       CASE( 'minmax' )
298
299          value      = fmin(1)
300          value_ijk  = fmin_ijk
301          value1     = fmax(1)
302          value1_ijk = fmax_ijk
303
304       CASE( 'abs', 'absoff' )
305
306          value     = fmax(1)
307          value_ijk = fmax_ijk
308          IF ( fmax_ijk(1) < 0 )  THEN
309             value        = -value
310             value_ijk(1) = -value_ijk(1) - 10         !???
311          ENDIF
312
313    END SELECT
314
315!
316!-- Limit index values to the range 0..nx, 0..ny
317    IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
318    IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
319    IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
320    IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
321
322
323 END SUBROUTINE global_min_max
Note: See TracBrowser for help on using the repository browser.