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

Last change on this file since 866 was 866, checked in by raasch, 12 years ago

bugfix for timestep calculation in case of Galilei transformation

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