Home

Context Navigation

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

Last change on this file since 4372 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 12.0 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-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: global_min_max.f90 4360 2020-01-07 11:25:50Z banzhafs $
27	! OpenACC support added
28	!
29	! 4182 2019-08-22 15:20:23Z scharf
30	! Corrected "Former revisions" section
31	!
32	! 3655 2019-01-07 16:51:22Z knoop
33	! Corrected "Former revisions" section
34	!
35	! Revision 1.1 1997/07/24 11:14:03 raasch
36	! Initial revision
37	!
38	!
39	! Description:
40	! ------------
41	!> Determine the array minimum/maximum and the corresponding indices.
42	!------------------------------------------------------------------------------!
43	SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
44	value_ijk, value1, value1_ijk )
45
46
47	USE indices, &
48	ONLY: nbgp, ny, nx
49
50	USE kinds
51
52	USE pegrid
53
54	IMPLICIT NONE
55
56	CHARACTER (LEN=*) :: mode !<
57
58	INTEGER(iwp) :: i !<
59	INTEGER(iwp) :: i1 !<
60	INTEGER(iwp) :: i2 !<
61	INTEGER(iwp) :: id_fmax !<
62	INTEGER(iwp) :: id_fmin !<
63	INTEGER(iwp) :: j !<
64	INTEGER(iwp) :: j1 !<
65	INTEGER(iwp) :: j2 !<
66	INTEGER(iwp) :: k !<
67	INTEGER(iwp) :: k1 !<
68	INTEGER(iwp) :: k2 !<
69	INTEGER(iwp) :: fmax_ijk(3) !<
70	INTEGER(iwp) :: fmax_ijk_l(3) !<
71	INTEGER(iwp) :: fmin_ijk(3) !<
72	INTEGER(iwp) :: fmin_ijk_l(3) !<
73	INTEGER(iwp) :: value_ijk(3) !<
74
75	INTEGER(iwp), OPTIONAL :: value1_ijk(3) !<
76
77	REAL(wp) :: offset !<
78	REAL(wp) :: value !<
79	REAL(wp) :: ar(i1:i2,j1:j2,k1:k2) !<
80
81	#if defined( __ibm )
82	REAL(sp) :: fmax(2) !<
83	REAL(sp) :: fmax_l(2) !<
84	REAL(sp) :: fmin(2) !<
85	REAL(sp) :: fmin_l(2) !<
86	! on 32bit-machines MPI_2REAL must not be replaced
87	! by MPI_2DOUBLE_PRECISION
88	#else
89	REAL(wp) :: fmax(2) !<
90	REAL(wp) :: fmax_l(2) !<
91	REAL(wp) :: fmin(2) !<
92	REAL(wp) :: fmin_l(2) !<
93	#endif
94	#if defined( _OPENACC )
95	REAL(wp) :: red !< scalar for reduction with OpenACC
96	INTEGER(iwp) :: count_eq !< counter for locations of maximum
97	#endif
98	REAL(wp), OPTIONAL :: value1 !<
99
100
101	!
102	!-- Determine array minimum
103	IF ( mode == 'min' .OR. mode == 'minmax' ) THEN
104
105	!
106	!-- Determine the local minimum
107	fmin_ijk_l = MINLOC( ar )
108	fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1 ! MINLOC assumes lowerbound = 1
109	fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
110	fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - nbgp
111	fmin_l(1) = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
112
113	#if defined( __parallel )
114	fmin_l(2) = myid
115	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
116	CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, &
117	ierr )
118
119	!
120	!-- Determine the global minimum. Result stored on PE0.
121	id_fmin = fmin(2)
122	IF ( id_fmin /= 0 ) THEN
123	IF ( myid == 0 ) THEN
124	CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, &
125	status, ierr )
126	ELSEIF ( myid == id_fmin ) THEN
127	CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
128	ENDIF
129	ELSE
130	fmin_ijk = fmin_ijk_l
131	ENDIF
132	!
133	!-- Send the indices of the just determined array minimum to other PEs
134	CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
135	#else
136	fmin(1) = fmin_l(1)
137	fmin_ijk = fmin_ijk_l
138	#endif
139
140	ENDIF
141
142	!
143	!-- Determine array maximum
144	IF ( mode == 'max' .OR. mode == 'minmax' ) THEN
145
146	!
147	!-- Determine the local maximum
148	fmax_ijk_l = MAXLOC( ar )
149	fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1 ! MAXLOC assumes lowerbound = 1
150	fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
151	fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - nbgp
152	fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
153
154	#if defined( __parallel )
155	fmax_l(2) = myid
156	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
157	CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
158	ierr )
159
160	!
161	!-- Determine the global maximum. Result stored on PE0.
162	id_fmax = fmax(2)
163	IF ( id_fmax /= 0 ) THEN
164	IF ( myid == 0 ) THEN
165	CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
166	status, ierr )
167	ELSEIF ( myid == id_fmax ) THEN
168	CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
169	ENDIF
170	ELSE
171	fmax_ijk = fmax_ijk_l
172	ENDIF
173	!
174	!-- send the indices of the just determined array maximum to other PEs
175	CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
176	#else
177	fmax(1) = fmax_l(1)
178	fmax_ijk = fmax_ijk_l
179	#endif
180
181	ENDIF
182
183	!
184	!-- Determine absolute array maximum
185	IF ( mode == 'abs' ) THEN
186
187	#if defined( _OPENACC )
188	red = 0.0_wp
189	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
190	!$ACC PRESENT(ar) COPY(red) REDUCTION(MAX: red)
191	DO k = k1, k2
192	DO j = j1, j2
193	DO i = i1, i2
194	IF ( ABS( ar(i,j,k) ) > red ) THEN
195	red = ABS( ar(i,j,k) )
196	ENDIF
197	ENDDO
198	ENDDO
199	ENDDO
200	fmax_l(1) = red
201
202	!
203	!-- Determine the maximum's position and count how often it is found.
204	count_eq = 0
205	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
206	!$ACC PRESENT(ar) COPY(fmax_ijk_l(1:3), count_eq) &
207	!$ACC REDUCTION(+:count_eq)
208	DO k = k1, k2
209	DO j = j1, j2
210	DO i = i1, i2
211	IF ( ABS( ar(i,j,k) ) == red ) THEN
212	fmax_ijk_l(1) = i
213	fmax_ijk_l(2) = j
214	fmax_ijk_l(3) = k
215	count_eq = count_eq + 1
216	ENDIF
217	ENDDO
218	ENDDO
219	ENDDO
220
221	IF ( count_eq == 1 ) THEN
222	!
223	!-- We found a single maximum element and correctly got its position. Transfer its
224	!-- value to handle the negative case correctly.
225	!$ACC UPDATE HOST(ar(fmax_ijk_l(1):fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)))
226	ELSE
227	!
228	!-- We found no maximum element (?) or multiple, so the position is not correct.
229	!-- Copy the full array to the host and determine the maximum sequentially...
230	!$ACC UPDATE HOST(ar(i1:i2,j1:j2,k1:k2))
231	#endif
232
233	!
234	!-- Determine the local absolut maximum
235	fmax_l(1) = 0.0_wp
236	fmax_ijk_l(1) = i1
237	fmax_ijk_l(2) = j1
238	fmax_ijk_l(3) = k1
239	DO k = k1, k2
240	DO j = j1, j2
241	DO i = i1, i2
242	IF ( ABS( ar(i,j,k) ) > fmax_l(1) ) THEN
243	fmax_l(1) = ABS( ar(i,j,k) )
244	fmax_ijk_l(1) = i
245	fmax_ijk_l(2) = j
246	fmax_ijk_l(3) = k
247	ENDIF
248	ENDDO
249	ENDDO
250	ENDDO
251
252	#if defined( _OPENACC )
253	!
254	!-- Close ELSE case from above
255	ENDIF
256	#endif
257
258	!
259	!-- Set a flag in case that the determined value is negative.
260	!-- A constant offset has to be subtracted in order to handle the special
261	!-- case i=0 correctly
262	IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp ) THEN
263	fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
264	ENDIF
265
266	#if defined( __parallel )
267	fmax_l(2) = myid
268	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
269	CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
270	ierr )
271
272	!
273	!-- Determine the global absolut maximum. Result stored on PE0.
274	id_fmax = fmax(2)
275	IF ( id_fmax /= 0 ) THEN
276	IF ( myid == 0 ) THEN
277	CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
278	status, ierr )
279	ELSEIF ( myid == id_fmax ) THEN
280	CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
281	ENDIF
282	ELSE
283	fmax_ijk = fmax_ijk_l
284	ENDIF
285	!
286	!-- Send the indices of the just determined absolut maximum to other PEs
287	CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
288	#else
289	fmax(1) = fmax_l(1)
290	fmax_ijk = fmax_ijk_l
291	#endif
292
293	ENDIF
294
295	!
296	!-- Determine absolute maximum of ( array - offset )
297	IF ( mode == 'absoff' ) THEN
298
299	!
300	!-- Determine the local absolut maximum
301	fmax_l(1) = 0.0_wp
302	fmax_ijk_l(1) = i1
303	fmax_ijk_l(2) = j1
304	fmax_ijk_l(3) = k1
305	DO k = k1, k2
306	DO j = j1, j2
307	!
308	!-- Attention: the lowest gridpoint is excluded here, because there
309	!-- --------- is no advection at nzb=0 and mode 'absoff' is only
310	!-- used for calculating u,v extrema for CFL-criteria
311	DO i = i1+1, i2
312	IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) ) THEN
313	fmax_l(1) = ABS( ar(i,j,k) - offset )
314	fmax_ijk_l(1) = i
315	fmax_ijk_l(2) = j
316	fmax_ijk_l(3) = k
317	ENDIF
318	ENDDO
319	ENDDO
320	ENDDO
321
322	!
323	!-- Set a flag in case that the determined value is negative.
324	!-- A constant offset has to be subtracted in order to handle the special
325	!-- case i=0 correctly
326	IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp ) THEN
327	fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
328	ENDIF
329
330	#if defined( __parallel )
331	fmax_l(2) = myid
332	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
333	CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
334	ierr )
335
336	!
337	!-- Determine the global absolut maximum. Result stored on PE0.
338	id_fmax = fmax(2)
339	IF ( id_fmax /= 0 ) THEN
340	IF ( myid == 0 ) THEN
341	CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
342	status, ierr )
343	ELSEIF ( myid == id_fmax ) THEN
344	CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
345	ENDIF
346	ELSE
347	fmax_ijk = fmax_ijk_l
348	ENDIF
349	!
350	!-- Send the indices of the just determined absolut maximum to other PEs
351	CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
352	#else
353	fmax(1) = fmax_l(1)
354	fmax_ijk = fmax_ijk_l
355	#endif
356
357	ENDIF
358
359	!
360	!-- Determine output parameters
361	SELECT CASE( mode )
362
363	CASE( 'min' )
364
365	value = fmin(1)
366	value_ijk = fmin_ijk
367
368	CASE( 'max' )
369
370	value = fmax(1)
371	value_ijk = fmax_ijk
372
373	CASE( 'minmax' )
374
375	value = fmin(1)
376	value_ijk = fmin_ijk
377	value1 = fmax(1)
378	value1_ijk = fmax_ijk
379
380	CASE( 'abs', 'absoff' )
381
382	value = fmax(1)
383	value_ijk = fmax_ijk
384	IF ( fmax_ijk(1) < 0 ) THEN
385	value = -value
386	value_ijk(1) = -value_ijk(1) - 10 !???
387	ENDIF
388
389	END SELECT
390
391	!
392	!-- Limit index values to the range 0..nx, 0..ny
393	IF ( value_ijk(3) < 0 ) value_ijk(3) = nx +1 + value_ijk(3)
394	IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
395	IF ( value_ijk(2) < 0 ) value_ijk(2) = ny +1 + value_ijk(2)
396	IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
397
398
399	END SUBROUTINE global_min_max

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |