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

Last change on this file since 1188 was 1188, checked in by heinze, 11 years ago

Bugfix in modes 'min' and 'max': x and z component were interchange

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