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

Last change on this file since 3933 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 11.1 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 3655 2019-01-07 16:51:22Z kanani $
27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
31!
32! 2101 2017-01-05 16:42:31Z suehring
33!
34! 2000 2016-08-20 18:09:15Z knoop
35! Forced header and separation lines into 80 columns
36!
37! 1682 2015-10-07 23:56:08Z knoop
38! Code annotations made doxygen readable
39!
40! 1353 2014-04-08 15:21:23Z heinze
41! REAL constants provided with KIND-attribute
42!
43! 1320 2014-03-20 08:40:49Z raasch
44! ONLY-attribute added to USE-statements,
45! kind-parameters added to all INTEGER and REAL declaration statements,
46! kinds are defined in new module kinds,
47! revision history before 2012 removed,
48! comment fields (!:) to be used for variable explanations added to
49! all variable declaration statements
50!
51! 1188 2013-06-20 12:00:08Z heinze
52! Bugfix in modes 'min' and 'max': x and z component were interchanged
53!
54! 1036 2012-10-22 13:43:42Z raasch
55! code put under GPL (PALM 3.9)
56!
57! 866 2012-03-28 06:44:41Z raasch
58! new mode "absoff" accounts for an offset in the respective array
59!
60! Revision 1.1  1997/07/24 11:14:03  raasch
61! Initial revision
62!
63!
64! Description:
65! ------------
66!> Determine the array minimum/maximum and the corresponding indices.
67!------------------------------------------------------------------------------!
68 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
69                            value_ijk, value1, value1_ijk )
70 
71
72    USE indices,                                                               &
73        ONLY:  nbgp, ny, nx
74       
75    USE kinds
76   
77    USE pegrid
78
79    IMPLICIT NONE
80
81    CHARACTER (LEN=*) ::  mode  !<
82
83    INTEGER(iwp) ::  i              !<
84    INTEGER(iwp) ::  i1             !<
85    INTEGER(iwp) ::  i2             !<
86    INTEGER(iwp) ::  id_fmax        !<
87    INTEGER(iwp) ::  id_fmin        !<
88    INTEGER(iwp) ::  j              !<
89    INTEGER(iwp) ::  j1             !<
90    INTEGER(iwp) ::  j2             !<
91    INTEGER(iwp) ::  k              !<
92    INTEGER(iwp) ::  k1             !<
93    INTEGER(iwp) ::  k2             !<
94    INTEGER(iwp) ::  fmax_ijk(3)    !<
95    INTEGER(iwp) ::  fmax_ijk_l(3)  !<
96    INTEGER(iwp) ::  fmin_ijk(3)    !<
97    INTEGER(iwp) ::  fmin_ijk_l(3)  !<
98    INTEGER(iwp) ::  value_ijk(3)   !<
99   
100    INTEGER(iwp), OPTIONAL ::  value1_ijk(3)  !<
101   
102    REAL(wp) ::  offset                 !<
103    REAL(wp) ::  value                  !<
104    REAL(wp) ::  ar(i1:i2,j1:j2,k1:k2)  !<
105   
106#if defined( __ibm )
107    REAL(sp) ::  fmax(2)    !<
108    REAL(sp) ::  fmax_l(2)  !<
109    REAL(sp) ::  fmin(2)    !<
110    REAL(sp) ::  fmin_l(2)  !<
111             ! on 32bit-machines MPI_2REAL must not be replaced
112             ! by MPI_2DOUBLE_PRECISION
113#else
114    REAL(wp) ::  fmax(2)    !<
115    REAL(wp) ::  fmax_l(2)  !<
116    REAL(wp) ::  fmin(2)    !<
117    REAL(wp) ::  fmin_l(2)  !<
118#endif
119    REAL(wp), OPTIONAL ::  value1  !<
120
121
122!
123!-- Determine array minimum
124    IF ( mode == 'min'  .OR.  mode == 'minmax' )  THEN
125
126!
127!--    Determine the local minimum
128       fmin_ijk_l = MINLOC( ar )
129       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1 ! MINLOC assumes lowerbound = 1
130       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
131       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - nbgp
132       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
133
134#if defined( __parallel )
135       fmin_l(2)  = myid
136       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
137       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, &
138                           ierr )
139
140!
141!--    Determine the global minimum. Result stored on PE0.
142       id_fmin = fmin(2)
143       IF ( id_fmin /= 0 )  THEN
144          IF ( myid == 0 )  THEN
145             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, &
146                            status, ierr )
147          ELSEIF ( myid == id_fmin )  THEN
148             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
149          ENDIF
150       ELSE
151          fmin_ijk = fmin_ijk_l
152       ENDIF
153!
154!--    Send the indices of the just determined array minimum to other PEs
155       CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
156#else
157       fmin(1)  = fmin_l(1)
158       fmin_ijk = fmin_ijk_l
159#endif
160
161    ENDIF
162
163!
164!-- Determine array maximum
165    IF ( mode == 'max'  .OR.  mode == 'minmax' )  THEN
166
167!
168!--    Determine the local maximum
169       fmax_ijk_l = MAXLOC( ar )
170       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1 ! MAXLOC assumes lowerbound = 1
171       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
172       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - nbgp
173       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
174
175#if defined( __parallel )
176       fmax_l(2)  = myid
177       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
178       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
179                           ierr )
180
181!
182!--    Determine the global maximum. Result stored on PE0.
183       id_fmax = fmax(2)
184       IF ( id_fmax /= 0 )  THEN
185          IF ( myid == 0 )  THEN
186             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
187                            status, ierr )
188          ELSEIF ( myid == id_fmax )  THEN
189             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
190          ENDIF
191       ELSE
192          fmax_ijk = fmax_ijk_l
193       ENDIF
194!
195!--    send the indices of the just determined array maximum to other PEs
196       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
197#else
198       fmax(1)  = fmax_l(1)
199       fmax_ijk = fmax_ijk_l
200#endif
201
202    ENDIF
203
204!
205!-- Determine absolute array maximum
206    IF ( mode == 'abs' )  THEN
207
208!
209!--    Determine the local absolut maximum
210       fmax_l(1)     = 0.0_wp
211       fmax_ijk_l(1) =  i1
212       fmax_ijk_l(2) =  j1
213       fmax_ijk_l(3) =  k1
214       DO  k = k1, k2
215          DO  j = j1, j2
216             DO  i = i1, i2
217                IF ( ABS( ar(i,j,k) ) > fmax_l(1) )  THEN
218                   fmax_l(1) = ABS( ar(i,j,k) )
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_wp )  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 absolute maximum of ( array - offset )
266    IF ( mode == 'absoff' )  THEN
267
268!
269!--    Determine the local absolut maximum
270       fmax_l(1)     = 0.0_wp
271       fmax_ijk_l(1) =  i1
272       fmax_ijk_l(2) =  j1
273       fmax_ijk_l(3) =  k1
274       DO  k = k1, k2
275          DO  j = j1, j2
276!
277!--          Attention: the lowest gridpoint is excluded here, because there
278!--          ---------  is no advection at nzb=0 and mode 'absoff' is only
279!--                     used for calculating u,v extrema for CFL-criteria
280             DO  i = i1+1, i2
281                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
282                   fmax_l(1) = ABS( ar(i,j,k) - offset )
283                   fmax_ijk_l(1) = i
284                   fmax_ijk_l(2) = j
285                   fmax_ijk_l(3) = k
286                ENDIF
287             ENDDO
288          ENDDO
289       ENDDO
290
291!
292!--    Set a flag in case that the determined value is negative.
293!--    A constant offset has to be subtracted in order to handle the special
294!--    case i=0 correctly
295       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
296          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
297       ENDIF
298
299#if defined( __parallel )
300       fmax_l(2)  = myid
301       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
302       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
303                           ierr )
304
305!
306!--    Determine the global absolut maximum. Result stored on PE0.
307       id_fmax = fmax(2)
308       IF ( id_fmax /= 0 )  THEN
309          IF ( myid == 0 )  THEN
310             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
311                            status, ierr )
312          ELSEIF ( myid == id_fmax )  THEN
313             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
314          ENDIF
315       ELSE
316          fmax_ijk = fmax_ijk_l
317       ENDIF
318!
319!--    Send the indices of the just determined absolut maximum to other PEs
320       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
321#else
322       fmax(1)  = fmax_l(1)
323       fmax_ijk = fmax_ijk_l
324#endif
325
326    ENDIF
327
328!
329!-- Determine output parameters
330    SELECT CASE( mode )
331
332       CASE( 'min' )
333
334          value     = fmin(1)
335          value_ijk = fmin_ijk
336
337       CASE( 'max' )
338
339          value     = fmax(1)
340          value_ijk = fmax_ijk
341
342       CASE( 'minmax' )
343
344          value      = fmin(1)
345          value_ijk  = fmin_ijk
346          value1     = fmax(1)
347          value1_ijk = fmax_ijk
348
349       CASE( 'abs', 'absoff' )
350
351          value     = fmax(1)
352          value_ijk = fmax_ijk
353          IF ( fmax_ijk(1) < 0 )  THEN
354             value        = -value
355             value_ijk(1) = -value_ijk(1) - 10         !???
356          ENDIF
357
358    END SELECT
359
360!
361!-- Limit index values to the range 0..nx, 0..ny
362    IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
363    IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
364    IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
365    IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
366
367
368 END SUBROUTINE global_min_max
Note: See TracBrowser for help on using the repository browser.