Line data Source code
1 : !> \file mo_utils.f90
2 : !> \brief \copybrief mo_utils
3 : !> \details \copydetails mo_utils
4 :
5 : !> \brief General utilities for the CHS library
6 : !> \details This module provides general utilities such as comparisons of two reals.
7 : !> \changelog
8 : !! - Matthias Cuntz, Juliane Mai, Feb 2014
9 : !! - written
10 : !! - Matthias Cuntz, Juliane Mai, Feb 2014
11 : !! - equal, notequal
12 : !! - Matthias Cuntz, May 2014
13 : !! - swap
14 : !! - Matthias Cuntz, May 2014
15 : !! - is_finite, is_nan, is_normal, special_value
16 : !> \authors Matthias Cuntz, Juliane Mai
17 : !> \date 2014 - 2016
18 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
19 : !! FORCES is released under the LGPLv3+ license \license_note
20 : MODULE mo_utils
21 :
22 : USE mo_kind, only : sp, dp, i1, i4, i8, spc, dpc
23 : USE mo_string_utils, only : toupper
24 :
25 : IMPLICIT NONE
26 :
27 : PUBLIC :: arange ! Natural numbers within interval
28 : PUBLIC :: cumsum ! Cumulative sum
29 : PUBLIC :: imaxloc ! maxloc(arr)(1)
30 : PUBLIC :: iminloc ! maxloc(arr)(1)
31 : PUBLIC :: linspace ! numpy like linspace
32 : PUBLIC :: equal ! a == b, a .eq. b
33 : PUBLIC :: greaterequal ! a >= b, a .ge. b
34 : PUBLIC :: lesserequal ! a <= b, a .le. b
35 : PUBLIC :: notequal ! a /= b, a .ne. b
36 : PUBLIC :: eq ! a == b, a .eq. b
37 : PUBLIC :: is_close ! a =~ b, a .eq. b in defined precision
38 : PUBLIC :: ge ! a >= b, a .ge. b
39 : PUBLIC :: le ! a <= b, a .le. b
40 : PUBLIC :: ne ! a /= b, a .ne. b
41 : PUBLIC :: is_finite ! .true. if not IEEE Inf and not IEEE NaN
42 : PUBLIC :: is_nan ! .true. if IEEE NaN
43 : PUBLIC :: is_normal ! .true. if not IEEE Inf and not IEEE NaN
44 : PUBLIC :: locate ! Find closest values in a monotonic series
45 : PUBLIC :: swap ! swaps arrays or elements of an array
46 : PUBLIC :: special_value ! Special IEEE values
47 : PUBLIC :: relational_operator_dp, relational_operator_sp ! abstract interface for relational operators
48 :
49 : public :: flip ! flips a dimension of an array
50 : public :: unpack_chunkwise ! flips a dimension of an array
51 :
52 : !> \brief flip an array at a certain dimension
53 : interface flip
54 : procedure flip_1D_dp, flip_2D_dp, flip_3D_dp, flip_4D_dp, flip_1D_i4, flip_2D_i4, flip_3D_i4, flip_4D_i4
55 : end interface
56 :
57 : !> \brief chunk version of the unpack operation
58 : interface unpack_chunkwise
59 : procedure unpack_chunkwise_i1, unpack_chunkwise_dp
60 : end interface
61 :
62 : ! ------------------------------------------------------------------
63 : !> \brief Numbers within a given range.
64 : !> \details Gives array with numbers in a given interval, i.e.
65 : !! \f[ arange(1) = lower \f]
66 : !! \f[ arange(2) = lower+1 \f]
67 : !! ...
68 : !! \f[ arange(n) = upper \f]
69 : !!
70 : !! Default is lower=1.
71 : !! Output array has kind of lower.
72 : !!
73 : !! \b Example
74 : !! \code{.f90}
75 : !! rr = arange(100._dp)
76 : !! \endcode
77 : !! -> see also example in test directory
78 :
79 : !> \param[in] "integer(i4/i8)/real(sp/dp) :: lower" Start of interval if upper is given,
80 : !! Otherwise end of interval and start of interval is 1.
81 : !> \param[in] "integer(i4/i8)/real(sp/dp) :: upper End of interval"
82 : !> \return kind(arr) :: arange(upper-lower+1) — 1D array with values within given interval.
83 :
84 : !> \authors Matthias Cuntz
85 : !> \date Jun 2016
86 : INTERFACE arange
87 : MODULE PROCEDURE arange_i4, arange_i8, arange_dp, arange_sp
88 : END INTERFACE arange
89 :
90 : ! ------------------------------------------------------------------
91 : !> \brief Cumulative sum.
92 : !> \details The cumulative sum of the elements of an array
93 : !! \f[ cumsum(i) = \sum_{j=1}^i array(j) \f]
94 : !!
95 : !! \b Example
96 : !! \code{.f90}
97 : !! vec = (/ 1., 2., 3., 4., 5., 6. /)
98 : !! cum = cumsum(vec)
99 : !! \endcode
100 : !! -> see also example in test directory
101 :
102 : !> \param[in] "integer(i4/i8)/real(sp/dp)/complex(spc/dpc) :: arr(:)" 1D array
103 : !> \return kind(arr) :: cumsum(size(arr)) — Cumulative sum
104 :
105 : !> \authors Matthias Cuntz
106 : !> \date Jun 2016
107 : INTERFACE cumsum
108 : MODULE PROCEDURE cumsum_i4, cumsum_i8, cumsum_dp, cumsum_sp, cumsum_dpc, cumsum_spc
109 : END INTERFACE cumsum
110 :
111 : ! ------------------------------------------------------------------
112 : !> \brief First location in array of element with the maximum value.
113 : !> \details Fortran intrinsic maxloc return arrays with all subsripts
114 : !! corresponding to the maximum value in the array.\n
115 : !! This routine returns only the first entry as scalar integer.
116 : !!
117 : !! \b Example
118 : !! \code{.f90}
119 : !! integer(i4) :: imax
120 : !! imax = imaxloc(vec, mask=mask)
121 : !! \endcode
122 : !! -> see also example in test directory
123 :
124 : !> \param[in] "integer(i4/i8)/real(sp/dp)/complex(spc/dpc) :: array(:)" Input array
125 : !> \param[in] "logical :: mask(:)" If present, only those locations in array corresponding to
126 : !! the true values in mask are searched for the maximum value.
127 : !> \return integer(i4) :: imaxloc — First location of maximum
128 :
129 : !> \authors Matthias Cuntz, Juliane Mai
130 : !> \date Feb 2014
131 : INTERFACE imaxloc
132 : MODULE PROCEDURE imaxloc_i4, imaxloc_i8, imaxloc_sp, imaxloc_dp
133 : END INTERFACE imaxloc
134 :
135 : ! ------------------------------------------------------------------
136 : !> \brief First location in array of element with the minimum value.
137 : !> \details Fortran intrinsic minloc return arrays with all subsripts
138 : !! corresponding to the minimum value in the array.\n
139 : !! This routine returns only the first entry as scalar integer.
140 : !!
141 : !! \b Example
142 : !! \code{.f90}
143 : !! integer(i4) :: imin
144 : !! imin = iminloc(vec, mask=mask)
145 : !! \endcode
146 : !! -> see also example in test directory
147 :
148 : !> \param[in] "integer(i4/i8)/real(sp/dp)/complex(spc/dpc) :: array(:)" Input array
149 : !> \param[in] "logical :: mask(:)" If present, only those locations in array corresponding to
150 : !! the true values in mask are searched for the maximum value.
151 : !> \return integer(i4) :: iminloc — First location of minimum
152 :
153 : !> \authors Matthias Cuntz, Juliane Mai
154 : !> \date Feb 2014
155 : INTERFACE iminloc
156 : MODULE PROCEDURE iminloc_i4, iminloc_i8, iminloc_sp, iminloc_dp
157 : END INTERFACE iminloc
158 :
159 : ! ------------------------------------------------------------------
160 : !> \brief Evenly spaced numbers in interval.
161 : !> \details Return N evenly spaced numbers over a specified interval [lower,upper].
162 : !! \f[ linspace(lower,upper,N) = lower + arange(0,N-1)/(N-1) * (upper-lower) \f]
163 : !> Output array has kind of lower.
164 : !!
165 : !! \b Example
166 : !! \code{.f90}
167 : !! rr = linspace(1.0_dp,11._dp,101)
168 : !! \endcode
169 : !! -> see also example in test directory
170 :
171 : !> \param[in] "integer(i4/i8)/real(sp/dp) :: lower" Start of interval.
172 : !> \param[in] "integer(i4/i8)/real(sp/dp) :: upper" End of interval.
173 : !> \param[in] "integer(i4) :: nstep" Number of steps.
174 : !> \return kind(lower) :: linspace(N) — 1D array with evenly spaced numbers between lower and upper.
175 :
176 : !> \authors Matthias Cuntz
177 : !> \date Jun 2016
178 : INTERFACE linspace
179 : MODULE PROCEDURE linspace_i4, linspace_i8, linspace_dp, linspace_sp
180 : END INTERFACE linspace
181 : ! ------------------------------------------------------------------
182 :
183 : !> \brief Comparison of real values.
184 : !> \details Compares two reals if they are numerically equal or not, with option for precision, i.e.
185 : !! equal: \f[ |a-b| < \mathrm{max} \{ \epsilon_\mathrm{rel} \mathrm{max} \{ |a|,|b| \}, \epsilon_\mathrm{abs} \} \f]
186 : !!
187 : !! \b Example
188 : !!
189 : !! Returns ´.false.´ in 5th element of isequal
190 : !! \code{.f90}
191 : !! vec1 = (/ 1., 2., 3., -999., 5., 6. /)
192 : !! vec2 = (/ 1., 1., 3., -999., 10., 6. /)
193 : !! isequal = is_close(vec1, vec2)
194 : !! \endcode
195 : !! Returns ´.true.´ in all elements of isequal
196 : !! \code{.f90}
197 : !! vec1 = (/ 1., 2., 3., -999., 5., 6. /)
198 : !! vec2 = (/ 1., 1., 3., -999., 10., 6. /)
199 : !! isequal = is_close(vec1, vec2, atol = 5.0_dp)
200 : !! \endcode
201 : !! Returns ´.false.´ in 6th element of isequal
202 : !! \code{.f90}
203 : !! vec1 = (/ 1., 2., 3., -999., 5., NaN /)
204 : !! vec2 = (/ 1., 1., 3., -999., 10., NaN /)
205 : !! isequal = is_close(vec1, vec2, atol = 5.0_dp)
206 : !! \endcode
207 : !! Returns ´.true.´ in all elements of isequal
208 : !! \code{.f90}
209 : !! vec1 = (/ 1., 2., 3., -999., 5., NaN /)
210 : !! vec2 = (/ 1., 1., 3., -999., 10., NaN /)
211 : !! isequal = is_close(vec1, vec2, equal_nan = .true.)
212 : !! \endcode
213 :
214 : !> \param[in] "real(sp/dp) :: a" First number to compare
215 : !> \param[in] "real(sp/dp) :: b" Second number to compare
216 : !> \param[in] "real(sp/dp), optional :: rtol" Relative tolerance, will scale with a (DEFAULT : 1.0E-5).
217 : !> \param[in] "real(sp/dp), optional :: atol" Absolute tolerance (DEFAULT : 1.0E-8).
218 : !> \param[in] "logical, optional :: equal_nan" If \f$.true.\f$, NaN values will between a and b are considered
219 : !! equal.
220 : !> \retval "real(sp/dp) :: is_close" \f$ a == b \f$ logically true or false
221 :
222 : !> \authors Arya Prasetya
223 : !> \date Jan 2022
224 : !! - add precision arguments and is_close
225 : INTERFACE is_close
226 : MODULE PROCEDURE is_close_sp, is_close_dp
227 : END INTERFACE is_close
228 :
229 : !> \brief Comparison of real values.
230 : !> \details Compares two reals if they are exactly numerically equal or not, i.e.
231 : !! equal: \f[ |\frac{a-b}{b}| < \epsilon \f]
232 : !!
233 : !! \b Example
234 : !!
235 : !! Returns ´.false.´ in 5th element of isequal
236 : !! \code{.f90}
237 : !! vec1 = (/ 1., 2., 3., -999., 5., 6. /)
238 : !! vec2 = (/ 1., 1., 3., -999., 10., 6. /)
239 : !! isequal = equal(vec1, vec2)
240 : !! \endcode
241 : !! Returns ´.true.´ in all elements of isequal
242 : !! \code{.f90}
243 : !! vec1 = (/ 1., 2., 3., -999., 5., 6. /)
244 : !! vec2 = (/ 1., 1., 3., -999., 10., 6. /)
245 : !! isequal = equal(vec1, vec2)
246 : !! \endcode
247 : !! Returns ´.false.´ in 6th element of isequal
248 : !! \code{.f90}
249 : !! vec1 = (/ 1., 2., 3., -999., 5., NaN /)
250 : !! vec2 = (/ 1., 1., 3., -999., 10., NaN /)
251 : !! isequal = equal(vec1, vec2)
252 : !! \endcode
253 :
254 : !> \param[in] "real(sp/dp) :: a" First number to compare
255 : !> \param[in] "real(sp/dp) :: b" Second number to compare
256 : !> \retval "real(sp/dp) :: equal" \f$ a == b \f$ logically true or false
257 :
258 : !> \authors Matthias Cuntz, Juliane Mai
259 : !> \date Feb 2014
260 : !! - sp, dp
261 : INTERFACE equal
262 : MODULE PROCEDURE equal_sp, equal_dp
263 : END INTERFACE equal
264 :
265 : !> \brief Comparison of real values for inequality.
266 : !> \see equal
267 : INTERFACE notequal
268 : MODULE PROCEDURE notequal_sp, notequal_dp
269 : END INTERFACE notequal
270 :
271 : !> \brief Comparison of real values: `a >= b`.
272 : !> \see equal
273 : INTERFACE greaterequal
274 : MODULE PROCEDURE greaterequal_sp, greaterequal_dp
275 : END INTERFACE greaterequal
276 :
277 : !> \brief Comparison of real values: `a <= b`.
278 : !> \see equal
279 : INTERFACE lesserequal
280 : MODULE PROCEDURE lesserequal_sp, lesserequal_dp
281 : END INTERFACE lesserequal
282 :
283 : !> \copydoc equal
284 : INTERFACE eq
285 : MODULE PROCEDURE equal_sp, equal_dp
286 : END INTERFACE eq
287 :
288 : !> \copydoc notequal
289 : INTERFACE ne
290 : MODULE PROCEDURE notequal_sp, notequal_dp
291 : END INTERFACE ne
292 :
293 : !> \copydoc greaterequal
294 : INTERFACE ge
295 : MODULE PROCEDURE greaterequal_sp, greaterequal_dp
296 : END INTERFACE ge
297 :
298 : !> \copydoc lesserequal
299 : INTERFACE le
300 : MODULE PROCEDURE lesserequal_sp, lesserequal_dp
301 : END INTERFACE le
302 :
303 :
304 : ! ------------------------------------------------------------------
305 :
306 : !> \brief .true. if not IEEE Inf.
307 : !> \details
308 : !! Checks for IEEE Inf, i.e. Infinity.\n
309 : !! Wraps to functions of the intrinsic module ieee_arithmetic.
310 : !!
311 : !! \b Example
312 : !!
313 : !! Returns `.false.` in 1st and 4th element.
314 : !! \code{.f90}
315 : !! vec1 = (/ NaN, 2., 3., Inf, 5., 6. /)
316 : !! isfinite = is_finite(vec1)
317 : !! \endcode
318 :
319 : !> \param[in] "real(sp/dp) :: a" Number to be evaluated.
320 : !> \retval "logical :: is_finite" \f$ a \neq \infty \f$, logically true or false.
321 :
322 : !> \authors Matthias Cuntz
323 : !> \date Mar 2015
324 : INTERFACE is_finite
325 : MODULE PROCEDURE is_finite_sp, is_finite_dp
326 : END INTERFACE is_finite
327 :
328 : !> \brief .true. if IEEE NaN.
329 : !> \details
330 : !! Checks for IEEE NaN, i.e. Not-a-Number.\n
331 : !! Wraps to functions of the intrinsic module ieee_arithmetic.
332 : !!
333 : !! \b Example
334 : !!
335 : !! Returns `.false.` in all but 1st element
336 : !! \code{.f90}
337 : !! vec1 = (/ NaN, 2., 3., Inf, 5., 6. /)
338 : !! isnan = is_nan(vec1)
339 : !! \endcode
340 :
341 : !> \param[in] "real(sp/dp) :: a" Number to be evaluated.
342 : !> \retval "logical :: is_nan" \f$ a = NaN \f$, logically true or false.
343 :
344 : INTERFACE is_nan
345 : MODULE PROCEDURE is_nan_sp, is_nan_dp
346 : END INTERFACE is_nan
347 :
348 : !> \brief .true. if nor IEEE Inf nor IEEE NaN.
349 : !> \details
350 : !! Checks if IEEE Inf and IEEE NaN, i.e. Infinity and Not-a-Number.\n
351 : !! Wraps to functions of the intrinsic module ieee_arithmetic.
352 : !!
353 : !! \b Example
354 : !!
355 : !! Returns `.true.` in all but 1st and 4th element.
356 : !! \code{.f90}
357 : !! vec1 = (/ NaN, 2., 3., Inf, 5., 6. /)
358 : !! isnormal = is_normal(vec1)
359 : !! \endcode
360 :
361 : !> \param[in] "real(sp/dp) :: a" Number to be evaluated.
362 : !> \retval "logical :: is_normal" \f$ a \neq \infty \land a = NaN \f$, logically true or false.
363 :
364 : INTERFACE is_normal
365 : MODULE PROCEDURE is_normal_sp, is_normal_dp
366 : END INTERFACE is_normal
367 :
368 :
369 : ! ------------------------------------------------------------------
370 :
371 : !> \brief Find closest values in a monotonic series, returns the indexes.
372 : !> \details
373 : !! Given an array x(1:n), and given a value y,
374 : !! returns a value j such that y is between
375 : !! x(j) and x(j+1).\n
376 : !!
377 : !! x must be monotonically increasing.\n
378 : !! j=0 or j=N is returned to indicate that x is out of range.
379 : !!
380 : !! \b Example
381 : !!
382 : !! Returns `ii = (/1, 5/)`
383 : !! \code{.f90}
384 : !! x = (/ 1., 2., 3., -999., 5., 6. /)
385 : !! y = (/ 1.1, 5.6 /)
386 : !! ii = locate(x, y)
387 : !! \endcode
388 : !!
389 : !! Returns `ii = 1`
390 : !! \code{.f90}
391 : !! y = 1.1
392 : !! ii = locate(x, y)
393 : !! \endcode
394 :
395 : !> \param[in] "real(dp/sp) :: x(:)" Sorted array
396 : !> \param[in] "real(dp/sp) :: y[(:)]" Value(s) of which the closest match in x(:) is wanted
397 : !> \retval "integer(i4) :: index[(:)]" Index(es) of x so that y is between x(index) and x(index+1)
398 :
399 : !> \note x must be monotonically increasing.\n
400 :
401 : !> \author Matthias Cuntz
402 : !> \date May 2014
403 : INTERFACE locate
404 : MODULE PROCEDURE locate_0d_dp, locate_0d_sp, locate_1d_dp, locate_1d_sp
405 : END INTERFACE locate
406 :
407 :
408 : ! ------------------------------------------------------------------
409 :
410 : !> \brief Swap to values or two elements in array.
411 : !> \details
412 : !! Swaps either two entities, i.e. scalars, vectors, matrices,
413 : !! or two elements in a vector.
414 : !! The call is either \n
415 : !! call swap(x,y) \n
416 : !! or \n
417 : !! call swap(vec,i,j)
418 : !!
419 : !! \b Example
420 : !!
421 : !! \code{.f90}
422 : !! vec1 = (/ 1., 2., 3., -999., 5., 6. /)
423 : !! vec2 = (/ 1., 1., 3., -999., 10., 6. /)
424 : !! \endcode
425 : !!
426 : !! Swaps elements in vec1 and vec2
427 : !! \code{.f90}
428 : !! call swap(vec1, vec2)
429 : !! \endcode
430 : !!
431 : !! Swaps 1st and 3rd element of vec1
432 : !! \code{.f90}
433 : !! call swap(vec1, 1, 3)
434 : !! \endcode
435 :
436 :
437 : !> \param[in] "integer(i4) :: i" Index of first element to be swapped with second [case swap(vec,i,j)]
438 : !> \param[in] "integer(i4) :: j" Index of second element to be swapped with first [case swap(vec,i,j)]
439 : !> \param[inout] "real(sp/dp/i4) :: x[(:,...)]" First scalar or array to swap with second [case swap(x,y)]
440 : !> \param[inout] "real(sp/dp/i4) :: y[(:[,:])]" Second scalar or array to swap with first [case swap(x,y)]
441 : !> \param[inout] "real(sp/dp/i4) :: x(:)" Vector of which to elements are swapped [case swap(vec,i,j)]
442 :
443 : !> \note No mask or undef.
444 :
445 : !> \author Matthias Cuntz
446 : !> \date May 2014
447 : INTERFACE swap
448 : MODULE PROCEDURE &
449 : swap_xy_dp, swap_xy_sp, swap_xy_i4, &
450 : swap_vec_dp, swap_vec_sp, swap_vec_i4
451 : END INTERFACE swap
452 :
453 :
454 : ! ------------------------------------------------------------------
455 :
456 : !> \brief Special IEEE values.
457 : !> \details
458 : !! Returns special IEEE values such as Infinity or Not-a-Number.\n
459 : !! Wraps to function ieee_value of the intrinsic module ieee_arithmetic.\n
460 : !! Current special values are:\n
461 : !! - IEEE_SIGNALING_NAN
462 : !! - IEEE_QUIET_NAN
463 : !! - IEEE_NEGATIVE_INF
464 : !! - IEEE_POSITIVE_INF
465 : !! - IEEE_NEGATIVE_DENORMAL
466 : !! - IEEE_POSITIVE_DENORMAL
467 : !! - IEEE_NEGATIVE_NORMAL
468 : !! - IEEE_POSITIVE_NORMAL
469 : !! - IEEE_NEGATIVE_ZERO
470 : !! - IEEE_POSITIVE_ZERO
471 : !!
472 : !! \b Example
473 : !!
474 : !! Returns NaN
475 : !! \code{.f90}
476 : !! NaN = special_value(1.0, 'IEEE_QUIET_NAN')
477 : !! nan = special_value(1.0_dp, 'ieee_quiet_nan')
478 : !! \endcode
479 :
480 : !> \param[in] "real(sp/dp) :: x" dummy for kind of output
481 : !> \param[in] "character(le=*) :: ieee" ieee signal nanme
482 : !> \retval "real(sp/dp) :: special_value" IEEE special value,
483 : !! IEEE_SIGNALING_NAN,
484 : !! IEEE_QUIET_NAN (==IEEE_SIGNALING_NAN for gfortran),
485 : !! IEEE_NEGATIVE_INF,
486 : !! IEEE_POSITIVE_INF,
487 : !! IEEE_NEGATIVE_DENORMAL (==-0.0 for gfortran),
488 : !! IEEE_POSITIVE_DENORMAL (==0.0 for gfortran),
489 : !! IEEE_NEGATIVE_NORMAL (==-1.0 for gfortran),
490 : !! IEEE_POSITIVE_NORMAL (==1.0 for gfortran),
491 : !! IEEE_NEGATIVE_ZERO,
492 : !! IEEE_POSITIVE_ZERO,
493 :
494 : !> \authors Matthias Cuntz
495 : !> \date Mar 2015
496 : INTERFACE special_value
497 : MODULE PROCEDURE special_value_sp, special_value_dp
498 : END INTERFACE special_value
499 :
500 : !> \brief abstract interface for a relational operator on double precision arguments
501 : abstract interface
502 : logical pure function relational_operator_dp(a, b) result(boolean)
503 : import dp
504 : real(dp), intent(in) :: a, b
505 : end function relational_operator_dp
506 : end interface
507 :
508 : !> \brief abstract interface for a relational operator on single precision arguments
509 : abstract interface
510 : logical pure function relational_operator_sp(a, b) result(boolean)
511 : import sp
512 : real(sp), intent(in) :: a, b
513 : end function relational_operator_sp
514 : end interface
515 :
516 : ! ------------------------------------------------------------------
517 :
518 : PRIVATE
519 :
520 : ! ------------------------------------------------------------------
521 :
522 : CONTAINS
523 :
524 : ! ------------------------------------------------------------------
525 :
526 1 : function arange_i4(lower, upper)
527 :
528 : implicit none
529 :
530 : integer(i4), intent(in) :: lower
531 : integer(i4), intent(in), optional :: upper
532 : integer(i4), dimension(:), allocatable :: arange_i4
533 :
534 : integer(i4) :: istart, istop
535 : integer(i4) :: i
536 :
537 1 : if (present(upper)) then
538 1 : istart = lower
539 1 : istop = upper
540 : else
541 0 : istart = 1_i4
542 0 : istop = lower
543 : endif
544 :
545 3 : allocate(arange_i4(istop-istart+1_i4))
546 :
547 101 : forall(i=istart:istop) arange_i4(i-istart+1) = i
548 :
549 1 : end function arange_i4
550 :
551 1 : function arange_i8(lower, upper)
552 :
553 : implicit none
554 :
555 : integer(i8), intent(in) :: lower
556 : integer(i8), intent(in), optional :: upper
557 : integer(i8), dimension(:), allocatable :: arange_i8
558 :
559 : integer(i8) :: istart, istop
560 : integer(i8) :: i
561 :
562 0 : if (present(upper)) then
563 0 : istart = lower
564 0 : istop = upper
565 : else
566 0 : istart = 1_i8
567 0 : istop = lower
568 : endif
569 :
570 0 : allocate(arange_i8(istop-istart+1_i8))
571 :
572 0 : forall(i=istart:istop) arange_i8(i-istart+1) = i
573 :
574 0 : end function arange_i8
575 :
576 25 : function arange_dp(lower, upper)
577 :
578 : implicit none
579 :
580 : real(dp), intent(in) :: lower
581 : real(dp), intent(in), optional :: upper
582 : real(dp), dimension(:), allocatable :: arange_dp
583 :
584 : integer(i8) :: istart, istop
585 : integer(i8) :: i
586 :
587 25 : if (present(upper)) then
588 22 : istart = int(lower,i8)
589 22 : istop = int(upper,i8)
590 : else
591 3 : istart = 1_i8
592 3 : istop = int(lower,i8)
593 : endif
594 :
595 75 : allocate(arange_dp(istop-istart+1_i8))
596 :
597 2545 : forall(i=istart:istop) arange_dp(i-istart+1) = real(i,dp)
598 :
599 25 : end function arange_dp
600 :
601 28 : function arange_sp(lower, upper)
602 :
603 : implicit none
604 :
605 : real(sp), intent(in) :: lower
606 : real(sp), intent(in), optional :: upper
607 : real(sp), dimension(:), allocatable :: arange_sp
608 :
609 : integer(i8) :: istart, istop
610 : integer(i8) :: i
611 :
612 3 : if (present(upper)) then
613 1 : istart = int(lower,i8)
614 1 : istop = int(upper,i8)
615 : else
616 2 : istart = 1_i8
617 2 : istop = int(lower,i8)
618 : endif
619 :
620 9 : allocate(arange_sp(istop-istart+1_i8))
621 :
622 303 : forall(i=istart:istop) arange_sp(i-istart+1) = real(i,sp)
623 :
624 3 : end function arange_sp
625 :
626 : ! ------------------------------------------------------------------
627 :
628 9 : function cumsum_i4(arr)
629 :
630 : implicit none
631 :
632 : integer(i4), dimension(:), intent(in) :: arr
633 : integer(i4), dimension(size(arr,1)) :: cumsum_i4
634 :
635 : integer(i4) :: i
636 :
637 3 : cumsum_i4(1) = arr(1)
638 300 : do i=2, size(arr)
639 300 : cumsum_i4(i) = cumsum_i4(i-1) + arr(i)
640 : end do
641 :
642 3 : end function cumsum_i4
643 :
644 3 : function cumsum_i8(arr)
645 :
646 : implicit none
647 :
648 : integer(i8), dimension(:), intent(in) :: arr
649 : integer(i8), dimension(size(arr,1)) :: cumsum_i8
650 :
651 : integer(i4) :: i
652 :
653 0 : cumsum_i8(1) = arr(1)
654 0 : do i=2, size(arr)
655 0 : cumsum_i8(i) = cumsum_i8(i-1) + arr(i)
656 : end do
657 :
658 0 : end function cumsum_i8
659 :
660 4 : function cumsum_dp(arr)
661 :
662 : implicit none
663 :
664 : real(dp), dimension(:), intent(in) :: arr
665 : real(dp), dimension(size(arr,1)) :: cumsum_dp
666 :
667 : integer(i4) :: i
668 :
669 2 : cumsum_dp(1) = arr(1)
670 200 : do i=2, size(arr)
671 200 : cumsum_dp(i) = cumsum_dp(i-1) + arr(i)
672 : end do
673 :
674 2 : end function cumsum_dp
675 :
676 4 : function cumsum_dpc(arr)
677 :
678 : implicit none
679 :
680 : complex(dpc), dimension(:), intent(in) :: arr
681 : complex(dpc), dimension(size(arr,1)) :: cumsum_dpc
682 :
683 : integer(i4) :: i
684 :
685 1 : cumsum_dpc(1) = arr(1)
686 100 : do i=2, size(arr)
687 100 : cumsum_dpc(i) = cumsum_dpc(i-1) + arr(i)
688 : end do
689 :
690 1 : end function cumsum_dpc
691 :
692 5 : function cumsum_sp(arr)
693 :
694 : implicit none
695 :
696 : real(sp), dimension(:), intent(in) :: arr
697 : real(sp), dimension(size(arr,1)) :: cumsum_sp
698 :
699 : integer(i4) :: i
700 :
701 2 : cumsum_sp(1) = arr(1)
702 200 : do i=2, size(arr)
703 200 : cumsum_sp(i) = cumsum_sp(i-1) + arr(i)
704 : end do
705 :
706 2 : end function cumsum_sp
707 :
708 2 : function cumsum_spc(arr)
709 :
710 : implicit none
711 :
712 : complex(spc), dimension(:), intent(in) :: arr
713 : complex(spc), dimension(size(arr,1)) :: cumsum_spc
714 :
715 : integer(i4) :: i
716 :
717 0 : cumsum_spc(1) = arr(1)
718 0 : do i=2, size(arr)
719 0 : cumsum_spc(i) = cumsum_spc(i-1) + arr(i)
720 : end do
721 :
722 0 : end function cumsum_spc
723 :
724 : ! ------------------------------------------------------------------
725 :
726 0 : function imaxloc_i4(arr, mask)
727 :
728 : implicit none
729 :
730 : integer(i4), dimension(:), intent(in) :: arr
731 : logical, dimension(:), intent(in), optional :: mask
732 : integer(i4) :: imaxloc_i4
733 :
734 : integer(i4), dimension(1) :: imax
735 :
736 0 : if (present(mask)) then
737 0 : imax = maxloc(arr, 1, mask)
738 : else
739 0 : imax = maxloc(arr, 1)
740 : endif
741 0 : imaxloc_i4 = imax(1)
742 :
743 0 : end function imaxloc_i4
744 :
745 0 : function imaxloc_i8(arr, mask)
746 :
747 : implicit none
748 :
749 : integer(i8), dimension(:), intent(in) :: arr
750 : logical, dimension(:), intent(in), optional :: mask
751 : integer(i4) :: imaxloc_i8
752 :
753 : integer(i4), dimension(1) :: imax
754 :
755 0 : if (present(mask)) then
756 0 : imax = maxloc(arr, 1, mask)
757 : else
758 0 : imax = maxloc(arr, 1)
759 : endif
760 0 : imaxloc_i8 = imax(1)
761 :
762 0 : end function imaxloc_i8
763 :
764 0 : function imaxloc_dp(arr, mask)
765 :
766 : implicit none
767 :
768 : real(dp), dimension(:), intent(in) :: arr
769 : logical, dimension(:), intent(in), optional :: mask
770 : integer(i4) :: imaxloc_dp
771 :
772 : integer(i4), dimension(1) :: imax
773 :
774 0 : if (present(mask)) then
775 0 : imax = maxloc(arr, 1, mask)
776 : else
777 0 : imax = maxloc(arr, 1)
778 : endif
779 0 : imaxloc_dp = imax(1)
780 :
781 0 : end function imaxloc_dp
782 :
783 0 : function imaxloc_sp(arr, mask)
784 :
785 : implicit none
786 :
787 : real(sp), dimension(:), intent(in) :: arr
788 : logical, dimension(:), intent(in), optional :: mask
789 : integer(i4) :: imaxloc_sp
790 :
791 : integer(i4), dimension(1) :: imax
792 :
793 0 : if (present(mask)) then
794 0 : imax = maxloc(arr, 1, mask)
795 : else
796 0 : imax = maxloc(arr, 1)
797 : endif
798 0 : imaxloc_sp = imax(1)
799 :
800 0 : end function imaxloc_sp
801 :
802 : ! ------------------------------------------------------------------
803 :
804 0 : function iminloc_i4(arr, mask)
805 :
806 : implicit none
807 :
808 : integer(i4), dimension(:), intent(in) :: arr
809 : logical, dimension(:), intent(in), optional :: mask
810 : integer(i4) :: iminloc_i4
811 :
812 : integer(i4), dimension(1) :: imin
813 :
814 0 : if (present(mask)) then
815 0 : imin = minloc(arr, 1, mask)
816 : else
817 0 : imin = minloc(arr, 1)
818 : endif
819 0 : iminloc_i4 = imin(1)
820 :
821 0 : end function iminloc_i4
822 :
823 0 : function iminloc_i8(arr, mask)
824 :
825 : implicit none
826 :
827 : integer(i8), dimension(:), intent(in) :: arr
828 : logical, dimension(:), intent(in), optional :: mask
829 : integer(i4) :: iminloc_i8
830 :
831 : integer(i4), dimension(1) :: imin
832 :
833 0 : if (present(mask)) then
834 0 : imin = minloc(arr, 1, mask)
835 : else
836 0 : imin = minloc(arr, 1)
837 : endif
838 0 : iminloc_i8 = imin(1)
839 :
840 0 : end function iminloc_i8
841 :
842 0 : function iminloc_dp(arr, mask)
843 :
844 : implicit none
845 :
846 : real(dp), dimension(:), intent(in) :: arr
847 : logical, dimension(:), intent(in), optional :: mask
848 : integer(i4) :: iminloc_dp
849 :
850 : integer(i4), dimension(1) :: imin
851 :
852 0 : if (present(mask)) then
853 0 : imin = minloc(arr, 1, mask)
854 : else
855 0 : imin = minloc(arr, 1)
856 : endif
857 0 : iminloc_dp = imin(1)
858 :
859 0 : end function iminloc_dp
860 :
861 0 : function iminloc_sp(arr, mask)
862 :
863 : implicit none
864 :
865 : real(sp), dimension(:), intent(in) :: arr
866 : logical, dimension(:), intent(in), optional :: mask
867 : integer(i4) :: iminloc_sp
868 :
869 : integer(i4), dimension(1) :: imin
870 :
871 0 : if (present(mask)) then
872 0 : imin = minloc(arr, 1, mask)
873 : else
874 0 : imin = minloc(arr, 1)
875 : endif
876 0 : iminloc_sp = imin(1)
877 :
878 0 : end function iminloc_sp
879 :
880 : ! ------------------------------------------------------------------
881 :
882 2 : function linspace_i4(lower, upper, nstep)
883 :
884 : implicit none
885 :
886 : integer(i4), intent(in) :: lower
887 : integer(i4), intent(in) :: upper
888 : integer(i4), intent(in) :: nstep
889 : integer(i4), dimension(nstep) :: linspace_i4
890 :
891 101 : linspace_i4 = lower + nint(arange(0.0_dp,real(nstep-1_i4,dp))/real(nstep-1_i4,dp) * real(upper-lower,dp), i4)
892 :
893 1 : end function linspace_i4
894 :
895 1 : function linspace_i8(lower, upper, nstep)
896 :
897 : implicit none
898 :
899 : integer(i8), intent(in) :: lower
900 : integer(i8), intent(in) :: upper
901 : integer(i4), intent(in) :: nstep
902 : integer(i8), dimension(nstep) :: linspace_i8
903 :
904 0 : linspace_i8 = lower + nint(arange(0.0_dp,real(nstep-1_i4,dp))/real(nstep-1_i4,dp) * real(upper-lower,dp), i8)
905 :
906 0 : end function linspace_i8
907 :
908 42 : function linspace_dp(lower, upper, nstep)
909 :
910 : implicit none
911 :
912 : real(dp), intent(in) :: lower
913 : real(dp), intent(in) :: upper
914 : integer(i4), intent(in) :: nstep
915 : real(dp), dimension(nstep) :: linspace_dp
916 :
917 2141 : linspace_dp = lower + arange(0.0_dp,real(nstep-1_i4,dp))/real(nstep-1_i4,dp) * (upper-lower)
918 :
919 21 : end function linspace_dp
920 :
921 23 : function linspace_sp(lower, upper, nstep)
922 :
923 : implicit none
924 :
925 : real(sp), intent(in) :: lower
926 : real(sp), intent(in) :: upper
927 : integer(i4), intent(in) :: nstep
928 : real(sp), dimension(nstep) :: linspace_sp
929 :
930 101 : linspace_sp = lower + arange(0.0_sp,real(nstep-1_i4,sp))/real(nstep-1_i4,sp) * (upper-lower)
931 :
932 1 : end function linspace_sp
933 :
934 : ! ------------------------------------------------------------------
935 :
936 6 : logical elemental pure function is_close_dp(a, b, rtol, atol, equal_nan) result(boolean)
937 :
938 : real(dp), intent(in) :: a
939 : real(dp), intent(in) :: b
940 : real(dp), intent(in), optional :: rtol, atol
941 : logical, intent(in), optional :: equal_nan
942 :
943 : real(dp) :: rt, at
944 : logical :: n
945 :
946 6 : rt = 1.0E-05_dp
947 6 : at = 1.0E-08_dp
948 6 : n = .false.
949 4 : if (present(rtol)) rt = rtol
950 6 : if (present(atol)) at = atol
951 6 : if (present(equal_nan)) n = equal_nan
952 :
953 6 : if ((rt < 0._dp).or.(at < 0._dp)) error stop
954 6 : boolean = (a == b)
955 6 : if (boolean) return
956 6 : boolean = (n.and.(is_nan_dp(a).or.is_nan_dp(b)))
957 6 : if (boolean) return
958 6 : if (.not.is_finite_dp(a) .or. .not.is_finite_dp(b)) return
959 :
960 6 : boolean = abs(a - b) <= max(rt * max(abs(a),abs(b)), at)
961 :
962 7 : end function is_close_dp
963 :
964 :
965 :
966 6 : logical elemental pure function is_close_sp(a, b, rtol, atol, equal_nan) result(boolean)
967 :
968 : real(sp), intent(in) :: a
969 : real(sp), intent(in) :: b
970 : real(sp), intent(in), optional :: rtol, atol
971 : logical, intent(in), optional :: equal_nan
972 :
973 : real(sp) :: rt, at
974 : logical :: n
975 :
976 6 : rt = 1.0E-05_sp
977 6 : at = 1.0E-08_sp
978 6 : n = .false.
979 4 : if (present(rtol)) rt = rtol
980 6 : if (present(atol)) at = atol
981 6 : if (present(equal_nan)) n = equal_nan
982 :
983 6 : if ((rt < 0._sp).or.(at < 0._sp)) error stop
984 6 : boolean = (a == b)
985 6 : if (boolean) return
986 5 : boolean = (n.and.(is_nan_sp(a).or.is_nan_sp(b)))
987 5 : if (boolean) return
988 5 : if (.not.is_finite_sp(a) .or. .not.is_finite_sp(b)) return
989 :
990 5 : boolean = abs(a - b) <= max(rt * max(abs(a),abs(b)), at)
991 :
992 11 : end function is_close_sp
993 :
994 : ! ------------------------------------------------------------------
995 :
996 2574989 : logical elemental pure function equal_dp(a, b) result(boolean)
997 :
998 : real(dp), intent(in) :: a
999 : real(dp), intent(in) :: b
1000 :
1001 2574989 : boolean = (a == b)
1002 2574989 : if (boolean) return
1003 2574933 : boolean = .not. ((epsilon(1.0_dp) * abs(b) - abs(a - b)) < 0.0_dp)
1004 :
1005 2574939 : end function equal_dp
1006 :
1007 :
1008 298 : logical elemental pure function equal_sp(a, b) result(boolean)
1009 :
1010 : real(sp), intent(in) :: a
1011 : real(sp), intent(in) :: b
1012 :
1013 298 : boolean = (a == b)
1014 298 : if (boolean) return
1015 280 : boolean = .not. ((epsilon(1.0_sp) * abs(b) - abs(a - b)) < 0.0_sp)
1016 :
1017 2575269 : end function equal_sp
1018 :
1019 : ! ------------------------------------------------------------------
1020 :
1021 2574564 : logical elemental pure function greaterequal_dp(a, b) result(boolean)
1022 :
1023 : real(dp), intent(in) :: a
1024 : real(dp), intent(in) :: b
1025 :
1026 2574564 : boolean = equal_dp(a, b).or.(a > b)
1027 :
1028 298 : end function greaterequal_dp
1029 :
1030 :
1031 14 : logical elemental pure function greaterequal_sp(a, b) result(boolean)
1032 :
1033 : real(sp), intent(in) :: a
1034 : real(sp), intent(in) :: b
1035 :
1036 14 : boolean = equal_sp(a, b).or.(a > b)
1037 :
1038 2574578 : end function greaterequal_sp
1039 :
1040 : ! ------------------------------------------------------------------
1041 :
1042 90 : logical elemental pure function lesserequal_dp(a, b) result(boolean)
1043 :
1044 : real(dp), intent(in) :: a
1045 : real(dp), intent(in) :: b
1046 :
1047 90 : boolean = equal_dp(a, b).or.(a < b)
1048 :
1049 104 : end function lesserequal_dp
1050 :
1051 :
1052 80 : logical elemental pure function lesserequal_sp(a, b) result(boolean)
1053 :
1054 : real(sp), intent(in) :: a
1055 : real(sp), intent(in) :: b
1056 :
1057 80 : boolean = equal_sp(a, b).or.(a < b)
1058 :
1059 170 : end function lesserequal_sp
1060 :
1061 : ! ------------------------------------------------------------------
1062 :
1063 242 : logical elemental pure function notequal_dp(a, b) result(boolean)
1064 :
1065 : real(dp), intent(in) :: a
1066 : real(dp), intent(in) :: b
1067 :
1068 242 : boolean = .not.equal_dp(a, b)
1069 :
1070 322 : end function notequal_dp
1071 :
1072 :
1073 202 : logical elemental pure function notequal_sp(a, b) result(boolean)
1074 :
1075 : real(sp), intent(in) :: a
1076 : real(sp), intent(in) :: b
1077 :
1078 202 : boolean = .not.equal_sp(a, b)
1079 :
1080 444 : end function notequal_sp
1081 :
1082 : ! ------------------------------------------------------------------
1083 :
1084 266903 : ELEMENTAL PURE FUNCTION is_finite_dp(a)
1085 :
1086 :
1087 202 : use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
1088 :
1089 : IMPLICIT NONE
1090 :
1091 : REAL(dp), INTENT(IN) :: a !< Number to be evaluated.
1092 : LOGICAL :: is_finite_dp !< logical :: is_finite — \f$ a \neq \infty \f$, logically true or false.
1093 :
1094 266903 : is_finite_dp = ieee_is_finite(a)
1095 :
1096 266903 : END FUNCTION is_finite_dp
1097 :
1098 118 : ELEMENTAL PURE FUNCTION is_finite_sp(a)
1099 :
1100 266903 : use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
1101 :
1102 : IMPLICIT NONE
1103 :
1104 : REAL(sp), INTENT(IN) :: a !< Number to be evaluated.
1105 : LOGICAL :: is_finite_sp !< logical :: is_finite — \f$ a \neq \infty \f$, logically true or false.
1106 :
1107 118 : is_finite_sp = ieee_is_finite(a)
1108 :
1109 118 : END FUNCTION is_finite_sp
1110 :
1111 :
1112 114 : ELEMENTAL PURE FUNCTION is_nan_dp(a)
1113 :
1114 118 : use, intrinsic :: ieee_arithmetic, only : isnan => ieee_is_nan
1115 :
1116 : IMPLICIT NONE
1117 :
1118 : REAL(dp), INTENT(IN) :: a !< Number to be evaluated.
1119 : LOGICAL :: is_nan_dp !< logical :: is_nan — \f$ a = NaN \f$, logically true or false.
1120 :
1121 114 : is_nan_dp = isnan(a)
1122 :
1123 114 : END FUNCTION is_nan_dp
1124 :
1125 112 : ELEMENTAL PURE FUNCTION is_nan_sp(a)
1126 :
1127 114 : use, intrinsic :: ieee_arithmetic, only : isnan => ieee_is_nan
1128 :
1129 : IMPLICIT NONE
1130 :
1131 : REAL(sp), INTENT(IN) :: a !< Number to be evaluated.
1132 : LOGICAL :: is_nan_sp !< logical :: is_nan — \f$ a = NaN \f$, logically true or false.
1133 :
1134 112 : is_nan_sp = isnan(a)
1135 :
1136 112 : END FUNCTION is_nan_sp
1137 :
1138 :
1139 100 : ELEMENTAL PURE FUNCTION is_normal_dp(a)
1140 :
1141 112 : use, intrinsic :: ieee_arithmetic, only : ieee_is_normal
1142 :
1143 : IMPLICIT NONE
1144 :
1145 : REAL(dp), INTENT(IN) :: a !< Number to be evaluated.
1146 : LOGICAL :: is_normal_dp !< logical :: is_normal — \f$ a \neq \infty \land a = NaN \f$, logically true or false.
1147 :
1148 100 : is_normal_dp = ieee_is_normal(a)
1149 :
1150 100 : END FUNCTION is_normal_dp
1151 :
1152 100 : ELEMENTAL PURE FUNCTION is_normal_sp(a)
1153 :
1154 100 : use, intrinsic :: ieee_arithmetic, only : ieee_is_normal
1155 :
1156 : IMPLICIT NONE
1157 :
1158 : REAL(sp), INTENT(IN) :: a !< Number to be evaluated.
1159 : LOGICAL :: is_normal_sp !< logical :: is_normal — \f$ a \neq \infty \land a = NaN \f$, logically true or false.
1160 :
1161 100 : is_normal_sp = ieee_is_normal(a)
1162 :
1163 100 : END FUNCTION is_normal_sp
1164 :
1165 : ! ------------------------------------------------------------------
1166 :
1167 : ! Given an array x(1:N), and given a value y, returns a value j such that y is between
1168 : ! x(j) and x(j+1). x must be monotonically increasing.
1169 : ! j=0 or j=N is returned to indicate that x is out of range.
1170 :
1171 1 : FUNCTION locate_0d_dp(x, y)
1172 :
1173 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
1174 : REAL(dp), INTENT(IN) :: y
1175 : INTEGER(i4) :: locate_0d_dp
1176 :
1177 : INTEGER(i4), dimension(1) :: c
1178 :
1179 101 : c = minloc(abs(x - y))
1180 1 : if (le(x(c(1)), y)) then
1181 : locate_0d_dp = c(1)
1182 : else
1183 0 : locate_0d_dp = c(1) - 1
1184 : end if
1185 :
1186 100 : END FUNCTION locate_0d_dp
1187 :
1188 1 : FUNCTION locate_0d_sp(x, y)
1189 :
1190 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
1191 : REAL(sp), INTENT(IN) :: y
1192 : INTEGER(i4) :: locate_0d_sp
1193 :
1194 : INTEGER(i4), dimension(1) :: c
1195 :
1196 101 : c = minloc(abs(x - y))
1197 1 : if (le(x(c(1)), y)) then
1198 : locate_0d_sp = c(1)
1199 : else
1200 0 : locate_0d_sp = c(1) - 1
1201 : end if
1202 :
1203 1 : END FUNCTION locate_0d_sp
1204 :
1205 3 : FUNCTION locate_1d_dp(x, y)
1206 :
1207 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
1208 : REAL(dp), DIMENSION(:), INTENT(IN) :: y
1209 : INTEGER(i4), DIMENSION(:), allocatable :: locate_1d_dp
1210 :
1211 : INTEGER(i4) :: ny, i
1212 : INTEGER(i4), dimension(1) :: c
1213 :
1214 :
1215 1 : ny = size(y)
1216 3 : if (.not. allocated(locate_1d_dp)) allocate(locate_1d_dp(ny))
1217 :
1218 6 : do i = 1, ny
1219 505 : c = minloc(abs(x - y(i)))
1220 6 : if (le(x(c(1)), y(i))) then
1221 4 : locate_1d_dp(i) = c(1)
1222 : else
1223 1 : locate_1d_dp(i) = c(1) - 1
1224 : end if
1225 : end do
1226 :
1227 1 : END FUNCTION locate_1d_dp
1228 :
1229 3 : FUNCTION locate_1d_sp(x, y)
1230 :
1231 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
1232 : REAL(sp), DIMENSION(:), INTENT(IN) :: y
1233 : INTEGER(i4), DIMENSION(:), allocatable :: locate_1d_sp
1234 :
1235 : INTEGER(i4) :: ny, i
1236 : INTEGER(i4), dimension(1) :: c
1237 :
1238 :
1239 1 : ny = size(y)
1240 3 : if (.not. allocated(locate_1d_sp)) allocate(locate_1d_sp(ny))
1241 :
1242 6 : do i = 1, ny
1243 505 : c = minloc(abs(x - y(i)))
1244 6 : if (le(x(c(1)), y(i))) then
1245 4 : locate_1d_sp(i) = c(1)
1246 : else
1247 1 : locate_1d_sp(i) = c(1) - 1
1248 : end if
1249 : end do
1250 :
1251 1 : END FUNCTION locate_1d_sp
1252 :
1253 : ! ------------------------------------------------------------------
1254 :
1255 100 : elemental pure subroutine swap_xy_dp(x, y)
1256 :
1257 : real(dp), intent(inout) :: x
1258 : real(dp), intent(inout) :: y
1259 :
1260 100 : real(dp) :: z
1261 :
1262 100 : z = x
1263 100 : x = y
1264 100 : y = z
1265 :
1266 1 : end subroutine swap_xy_dp
1267 :
1268 100 : elemental pure subroutine swap_xy_sp(x, y)
1269 :
1270 : real(sp), intent(inout) :: x
1271 : real(sp), intent(inout) :: y
1272 :
1273 100 : real(sp) :: z
1274 :
1275 100 : z = x
1276 100 : x = y
1277 100 : y = z
1278 :
1279 200 : end subroutine swap_xy_sp
1280 :
1281 100 : elemental pure subroutine swap_xy_i4(x, y)
1282 :
1283 : integer(i4), intent(inout) :: x
1284 : integer(i4), intent(inout) :: y
1285 :
1286 : integer(i4) :: z
1287 :
1288 100 : z = x
1289 100 : x = y
1290 100 : y = z
1291 :
1292 200 : end subroutine swap_xy_i4
1293 :
1294 :
1295 2 : subroutine swap_vec_dp(x, i1, i2)
1296 :
1297 : real(dp), dimension(:), intent(inout) :: x
1298 : integer(i4), intent(in) :: i1
1299 : integer(i4), intent(in) :: i2
1300 :
1301 2 : real(dp) :: z
1302 :
1303 2 : z = x(i1)
1304 2 : x(i1) = x(i2)
1305 2 : x(i2) = z
1306 :
1307 102 : end subroutine swap_vec_dp
1308 :
1309 2 : subroutine swap_vec_sp(x, i1, i2)
1310 :
1311 : real(sp), dimension(:), intent(inout) :: x
1312 : integer(i4), intent(in) :: i1
1313 : integer(i4), intent(in) :: i2
1314 :
1315 2 : real(sp) :: z
1316 :
1317 2 : z = x(i1)
1318 2 : x(i1) = x(i2)
1319 2 : x(i2) = z
1320 :
1321 2 : end subroutine swap_vec_sp
1322 :
1323 2 : subroutine swap_vec_i4(x, i1, i2)
1324 :
1325 : integer(i4), dimension(:), intent(inout) :: x
1326 : integer(i4), intent(in) :: i1
1327 : integer(i4), intent(in) :: i2
1328 :
1329 : integer(i4) :: z
1330 :
1331 2 : z = x(i1)
1332 2 : x(i1) = x(i2)
1333 2 : x(i2) = z
1334 :
1335 2 : end subroutine swap_vec_i4
1336 :
1337 : ! ------------------------------------------------------------------
1338 :
1339 20 : function special_value_dp(x, ieee)
1340 :
1341 2 : use, intrinsic :: ieee_arithmetic, only : ieee_value, &
1342 : IEEE_SIGNALING_NAN, &
1343 : IEEE_QUIET_NAN, &
1344 : IEEE_NEGATIVE_INF, &
1345 : IEEE_POSITIVE_INF, &
1346 : IEEE_NEGATIVE_DENORMAL, &
1347 : IEEE_POSITIVE_DENORMAL, &
1348 : IEEE_NEGATIVE_NORMAL, &
1349 : IEEE_POSITIVE_NORMAL, &
1350 : IEEE_NEGATIVE_ZERO, &
1351 : IEEE_POSITIVE_ZERO
1352 :
1353 : implicit none
1354 :
1355 : real(dp), intent(in) :: x !< dummy for kind of output.
1356 : character(len = *), intent(in) :: ieee !< ieee signal name.
1357 : real(dp) :: special_value_dp !< real(dp) :: special_value — IEEE special value,
1358 : !! IEEE_SIGNALING_NAN,
1359 : !! IEEE_QUIET_NAN,
1360 : !! IEEE_NEGATIVE_INF,
1361 : !! IEEE_POSITIVE_INF,
1362 : !! IEEE_NEGATIVE_DENORMAL,
1363 : !! IEEE_POSITIVE_DENORMAL,
1364 : !! IEEE_NEGATIVE_NORMAL,
1365 : !! IEEE_POSITIVE_NORMAL,
1366 : !! IEEE_NEGATIVE_ZERO,
1367 : !! IEEE_POSITIVE_ZERO,
1368 :
1369 : ! local
1370 : character(len = 21) :: ieee_up
1371 :
1372 20 : ieee_up = toupper(ieee)
1373 40 : select case(trim(ieee_up))
1374 : case('IEEE_SIGNALING_NAN')
1375 1 : special_value_dp = ieee_value(x, IEEE_SIGNALING_NAN)
1376 : case('IEEE_QUIET_NAN')
1377 1 : special_value_dp = ieee_value(x, IEEE_QUIET_NAN)
1378 : case('IEEE_NEGATIVE_INF')
1379 2 : special_value_dp = ieee_value(x, IEEE_NEGATIVE_INF)
1380 : case('IEEE_POSITIVE_INF')
1381 2 : special_value_dp = ieee_value(x, IEEE_POSITIVE_INF)
1382 : case('IEEE_NEGATIVE_DENORMAL')
1383 0 : special_value_dp = ieee_value(x, IEEE_NEGATIVE_DENORMAL)
1384 : case('IEEE_POSITIVE_DENORMAL')
1385 0 : special_value_dp = ieee_value(x, IEEE_POSITIVE_DENORMAL)
1386 : case('IEEE_NEGATIVE_NORMAL')
1387 2 : special_value_dp = ieee_value(x, IEEE_NEGATIVE_NORMAL)
1388 : case('IEEE_POSITIVE_NORMAL')
1389 2 : special_value_dp = ieee_value(x, IEEE_POSITIVE_NORMAL)
1390 : case('IEEE_NEGATIVE_ZERO')
1391 3 : special_value_dp = ieee_value(x, IEEE_NEGATIVE_ZERO)
1392 : case('IEEE_POSITIVE_ZERO')
1393 3 : special_value_dp = ieee_value(x, IEEE_POSITIVE_ZERO)
1394 : case default
1395 40 : special_value_dp = 0.0_dp
1396 : end select
1397 :
1398 20 : end function special_value_dp
1399 :
1400 20 : function special_value_sp(x, ieee)
1401 :
1402 20 : use, intrinsic :: ieee_arithmetic, only : ieee_value, &
1403 : IEEE_SIGNALING_NAN, &
1404 : IEEE_QUIET_NAN, &
1405 : IEEE_NEGATIVE_INF, &
1406 : IEEE_POSITIVE_INF, &
1407 : IEEE_NEGATIVE_DENORMAL, &
1408 : IEEE_POSITIVE_DENORMAL, &
1409 : IEEE_NEGATIVE_NORMAL, &
1410 : IEEE_POSITIVE_NORMAL, &
1411 : IEEE_NEGATIVE_ZERO, &
1412 : IEEE_POSITIVE_ZERO
1413 :
1414 : implicit none
1415 :
1416 : real(sp), intent(in) :: x !< dummy for kind of output.
1417 : character(len = *), intent(in) :: ieee !< ieee signal name.
1418 : real(sp) :: special_value_sp !< IEEE special value,
1419 : !! IEEE_SIGNALING_NAN,
1420 : !! IEEE_QUIET_NAN,
1421 : !! IEEE_NEGATIVE_INF,
1422 : !! IEEE_POSITIVE_INF,
1423 : !! IEEE_NEGATIVE_DENORMAL,
1424 : !! IEEE_POSITIVE_DENORMAL,
1425 : !! IEEE_NEGATIVE_NORMAL,
1426 : !! IEEE_POSITIVE_NORMAL,
1427 : !! IEEE_NEGATIVE_ZERO,
1428 : !! IEEE_POSITIVE_ZERO,
1429 :
1430 : ! local
1431 : character(len = 21) :: ieee_up
1432 :
1433 20 : ieee_up = toupper(ieee)
1434 40 : select case(trim(ieee_up))
1435 : case('IEEE_SIGNALING_NAN')
1436 1 : special_value_sp = ieee_value(x, IEEE_SIGNALING_NAN)
1437 : case('IEEE_QUIET_NAN')
1438 1 : special_value_sp = ieee_value(x, IEEE_QUIET_NAN)
1439 : case('IEEE_NEGATIVE_INF')
1440 2 : special_value_sp = ieee_value(x, IEEE_NEGATIVE_INF)
1441 : case('IEEE_POSITIVE_INF')
1442 2 : special_value_sp = ieee_value(x, IEEE_POSITIVE_INF)
1443 : case('IEEE_NEGATIVE_DENORMAL')
1444 0 : special_value_sp = ieee_value(x, IEEE_NEGATIVE_DENORMAL)
1445 : case('IEEE_POSITIVE_DENORMAL')
1446 0 : special_value_sp = ieee_value(x, IEEE_POSITIVE_DENORMAL)
1447 : case('IEEE_NEGATIVE_NORMAL')
1448 2 : special_value_sp = ieee_value(x, IEEE_NEGATIVE_NORMAL)
1449 : case('IEEE_POSITIVE_NORMAL')
1450 2 : special_value_sp = ieee_value(x, IEEE_POSITIVE_NORMAL)
1451 : case('IEEE_NEGATIVE_ZERO')
1452 3 : special_value_sp = ieee_value(x, IEEE_NEGATIVE_ZERO)
1453 : case('IEEE_POSITIVE_ZERO')
1454 3 : special_value_sp = ieee_value(x, IEEE_POSITIVE_ZERO)
1455 : case default
1456 40 : special_value_sp = 0.0_sp
1457 : end select
1458 :
1459 20 : end function special_value_sp
1460 :
1461 0 : subroutine flip_1D_dp(data, iDim)
1462 20 : use mo_string_utils, only: compress, num2str
1463 : use mo_message, only: message
1464 : real(dp), dimension(:), allocatable, intent(inout) :: data
1465 : integer(i4), intent(in) :: iDim
1466 :
1467 0 : real(dp), dimension(:), allocatable :: temp_data
1468 : integer(i4) :: iDim1
1469 :
1470 0 : if (iDim > 1_i4) then
1471 0 : call message('Cannot flip 1D-array at dimension ', compress(trim(num2str(iDim))))
1472 0 : stop 1
1473 : end if
1474 0 : allocate(temp_data(size(data, 1)))
1475 :
1476 0 : do iDim1 = 1, size(data, 1)
1477 0 : temp_data(size(data, 1) - iDim1 + 1) = data(iDim1)
1478 : end do
1479 0 : call move_alloc(temp_data, data)
1480 0 : end subroutine flip_1D_dp
1481 :
1482 0 : subroutine flip_2D_dp(data, iDim)
1483 0 : use mo_string_utils, only: compress, num2str
1484 : use mo_message, only: message
1485 :
1486 : real(dp), dimension(:, :), allocatable, intent(inout) :: data
1487 : integer(i4), intent(in) :: iDim
1488 :
1489 0 : real(dp), dimension(:, :), allocatable :: temp_data
1490 : integer(i4) :: iDim2, iDim1
1491 :
1492 0 : if (iDim > 2_i4) then
1493 0 : call message('Cannot flip 2D-array at dimension ', compress(trim(num2str(iDim))))
1494 0 : stop 1
1495 : end if
1496 :
1497 0 : allocate(temp_data(size(data, 1), size(data, 2)))
1498 :
1499 0 : if (iDim == 1_i4) then
1500 0 : do iDim2 = 1, size(data, 2)
1501 0 : do iDim1 = 1, size(data, 1)
1502 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2) = data(iDim1, iDim2)
1503 : end do
1504 : end do
1505 0 : else if (iDim == 2_i4) then
1506 0 : do iDim2 = 1, size(data, 2)
1507 0 : temp_data(:, size(data, 2) - iDim2 + 1) = data(:, iDim2)
1508 : end do
1509 : end if
1510 0 : call move_alloc(temp_data, data)
1511 0 : end subroutine flip_2D_dp
1512 :
1513 0 : subroutine flip_3D_dp(data, iDim)
1514 0 : use mo_string_utils, only: compress, num2str
1515 : use mo_message, only: message
1516 : real(dp), dimension(:, :, :), allocatable, intent(inout) :: data
1517 : integer(i4), intent(in) :: iDim
1518 :
1519 0 : real(dp), dimension(:, :, :), allocatable :: temp_data
1520 : integer(i4) :: iDim3, iDim2, iDim1
1521 :
1522 0 : if (iDim > 3_i4) then
1523 0 : call message('Cannot flip 3D-array at dimension ', compress(trim(num2str(iDim))))
1524 0 : stop 1
1525 : end if
1526 :
1527 0 : allocate(temp_data(size(data, 1), size(data, 2), size(data, 3)))
1528 :
1529 0 : if (iDim == 1_i4) then
1530 0 : do iDim3 = 1, size(data, 3)
1531 0 : do iDim2 = 1, size(data, 2)
1532 0 : do iDim1 = 1, size(data, 1)
1533 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2, iDim3) = data(iDim1, iDim2, iDim3)
1534 : end do
1535 : end do
1536 : end do
1537 0 : else if (iDim == 2_i4) then
1538 0 : do iDim3 = 1, size(data, 3)
1539 0 : do iDim2 = 1, size(data, 2)
1540 0 : temp_data(:, size(data, 2) - iDim2 + 1, iDim3) = data(:, iDim2, iDim3)
1541 : end do
1542 : end do
1543 0 : else if (iDim == 3_i4) then
1544 0 : do iDim3 = 1, size(data, 3)
1545 0 : temp_data(:, :, size(data, 3) - iDim3 + 1) = data(:, :, iDim3)
1546 : end do
1547 : end if
1548 0 : call move_alloc(temp_data, data)
1549 0 : end subroutine flip_3D_dp
1550 :
1551 0 : subroutine flip_4D_dp(data, iDim)
1552 0 : use mo_string_utils, only: compress, num2str
1553 : use mo_message, only: message
1554 : real(dp), dimension(:, :, :, :), allocatable, intent(inout) :: data
1555 : integer(i4), intent(in) :: iDim
1556 :
1557 0 : real(dp), dimension(:, :, :, :), allocatable :: temp_data
1558 : integer(i4) :: iDim4, iDim3, iDim2, iDim1
1559 :
1560 0 : if (iDim > 4_i4) then
1561 0 : call message('Cannot flip 4D-array at dimension ', compress(trim(num2str(iDim))))
1562 0 : stop 1
1563 : end if
1564 :
1565 0 : allocate(temp_data(size(data, 1), size(data, 2), size(data, 3), size(data, 4)))
1566 :
1567 0 : if (iDim == 1_i4) then
1568 0 : do iDim4 = 1, size(data, 4)
1569 0 : do iDim3 = 1, size(data, 3)
1570 0 : do iDim2 = 1, size(data, 2)
1571 0 : do iDim1 = 1, size(data, 1)
1572 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2, iDim3, iDim4) = data(iDim1, iDim2, iDim3, iDim4)
1573 : end do
1574 : end do
1575 : end do
1576 : end do
1577 0 : else if (iDim == 2_i4) then
1578 0 : do iDim4 = 1, size(data, 4)
1579 0 : do iDim3 = 1, size(data, 3)
1580 0 : do iDim2 = 1, size(data, 2)
1581 0 : temp_data(:, size(data, 2) - iDim2 + 1, iDim3, iDim4) = data(:, iDim2, iDim3, iDim4)
1582 : end do
1583 : end do
1584 : end do
1585 0 : else if (iDim == 3_i4) then
1586 0 : do iDim4 = 1, size(data, 4)
1587 0 : do iDim3 = 1, size(data, 3)
1588 0 : temp_data(:, :, size(data, 3) - iDim3 + 1, iDim4) = data(:, :, iDim3, iDim4)
1589 : end do
1590 : end do
1591 0 : else if (iDim == 4_i4) then
1592 0 : do iDim4 = 1, size(data, 4)
1593 0 : temp_data(:, :, :, size(data, 4) - iDim4 + 1) = data(:, :, :, iDim4)
1594 : end do
1595 : end if
1596 0 : call move_alloc(temp_data, data)
1597 0 : end subroutine flip_4D_dp
1598 :
1599 0 : subroutine flip_1D_i4(data, iDim)
1600 0 : use mo_string_utils, only: compress, num2str
1601 : use mo_message, only: message
1602 : integer(i4), dimension(:), allocatable, intent(inout) :: data
1603 : integer(i4), intent(in) :: iDim
1604 :
1605 0 : integer(i4), dimension(:), allocatable :: temp_data
1606 : integer(i4) :: iDim1
1607 :
1608 0 : if (iDim > 1_i4) then
1609 0 : call message('Cannot flip 1D-array at dimension ', compress(trim(num2str(iDim))))
1610 0 : stop 1
1611 : end if
1612 0 : allocate(temp_data(size(data, 1)))
1613 :
1614 0 : do iDim1 = 1, size(data, 1)
1615 0 : temp_data(size(data, 1) - iDim1 + 1) = data(iDim1)
1616 : end do
1617 0 : call move_alloc(temp_data, data)
1618 0 : end subroutine flip_1D_i4
1619 :
1620 0 : subroutine flip_2D_i4(data, iDim)
1621 0 : use mo_string_utils, only: compress, num2str
1622 : use mo_message, only: message
1623 :
1624 : integer(i4), dimension(:, :), allocatable, intent(inout) :: data
1625 : integer(i4), intent(in) :: iDim
1626 :
1627 0 : integer(i4), dimension(:, :), allocatable :: temp_data
1628 : integer(i4) :: iDim2, iDim1
1629 :
1630 0 : if (iDim > 2_i4) then
1631 0 : call message('Cannot flip 2D-array at dimension ', compress(trim(num2str(iDim))))
1632 0 : stop 1
1633 : end if
1634 :
1635 0 : allocate(temp_data(size(data, 1), size(data, 2)))
1636 :
1637 0 : if (iDim == 1_i4) then
1638 0 : do iDim2 = 1, size(data, 2)
1639 0 : do iDim1 = 1, size(data, 1)
1640 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2) = data(iDim1, iDim2)
1641 : end do
1642 : end do
1643 0 : else if (iDim == 2_i4) then
1644 0 : do iDim2 = 1, size(data, 2)
1645 0 : temp_data(:, size(data, 2) - iDim2 + 1) = data(:, iDim2)
1646 : end do
1647 : end if
1648 0 : call move_alloc(temp_data, data)
1649 0 : end subroutine flip_2D_i4
1650 :
1651 0 : subroutine flip_3D_i4(data, iDim)
1652 0 : use mo_string_utils, only: compress, num2str
1653 : use mo_message, only: message
1654 : integer(i4), dimension(:, :, :), allocatable, intent(inout) :: data
1655 : integer(i4), intent(in) :: iDim
1656 :
1657 0 : integer(i4), dimension(:, :, :), allocatable :: temp_data
1658 : integer(i4) :: iDim3, iDim2, iDim1
1659 :
1660 0 : if (iDim > 3_i4) then
1661 0 : call message('Cannot flip 3D-array at dimension ', compress(trim(num2str(iDim))))
1662 0 : stop 1
1663 : end if
1664 :
1665 0 : allocate(temp_data(size(data, 1), size(data, 2), size(data, 3)))
1666 :
1667 0 : if (iDim == 1_i4) then
1668 0 : do iDim3 = 1, size(data, 3)
1669 0 : do iDim2 = 1, size(data, 2)
1670 0 : do iDim1 = 1, size(data, 1)
1671 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2, iDim3) = data(iDim1, iDim2, iDim3)
1672 : end do
1673 : end do
1674 : end do
1675 0 : else if (iDim == 2_i4) then
1676 0 : do iDim3 = 1, size(data, 3)
1677 0 : do iDim2 = 1, size(data, 2)
1678 0 : temp_data(:, size(data, 2) - iDim2 + 1, iDim3) = data(:, iDim2, iDim3)
1679 : end do
1680 : end do
1681 0 : else if (iDim == 3_i4) then
1682 0 : do iDim3 = 1, size(data, 3)
1683 0 : temp_data(:, :, size(data, 3) - iDim3 + 1) = data(:, :, iDim3)
1684 : end do
1685 : end if
1686 0 : call move_alloc(temp_data, data)
1687 0 : end subroutine flip_3D_i4
1688 :
1689 0 : subroutine flip_4D_i4(data, iDim)
1690 0 : use mo_string_utils, only: compress, num2str
1691 : use mo_message, only: message
1692 : integer(i4), dimension(:, :, :, :), allocatable, intent(inout) :: data
1693 : integer(i4), intent(in) :: iDim
1694 :
1695 0 : integer(i4), dimension(:, :, :, :), allocatable :: temp_data
1696 : integer(i4) :: iDim4, iDim3, iDim2, iDim1
1697 :
1698 0 : if (iDim > 4_i4) then
1699 0 : call message('Cannot flip 4D-array at dimension ', compress(trim(num2str(iDim))))
1700 0 : stop 1
1701 : end if
1702 :
1703 0 : allocate(temp_data(size(data, 1), size(data, 2), size(data, 3), size(data, 4)))
1704 :
1705 0 : if (iDim == 1_i4) then
1706 0 : do iDim4 = 1, size(data, 4)
1707 0 : do iDim3 = 1, size(data, 3)
1708 0 : do iDim2 = 1, size(data, 2)
1709 0 : do iDim1 = 1, size(data, 1)
1710 0 : temp_data(size(data, 1) - iDim1 + 1, iDim2, iDim3, iDim4) = data(iDim1, iDim2, iDim3, iDim4)
1711 : end do
1712 : end do
1713 : end do
1714 : end do
1715 0 : else if (iDim == 2_i4) then
1716 0 : do iDim4 = 1, size(data, 4)
1717 0 : do iDim3 = 1, size(data, 3)
1718 0 : do iDim2 = 1, size(data, 2)
1719 0 : temp_data(:, size(data, 2) - iDim2 + 1, iDim3, iDim4) = data(:, iDim2, iDim3, iDim4)
1720 : end do
1721 : end do
1722 : end do
1723 0 : else if (iDim == 3_i4) then
1724 0 : do iDim4 = 1, size(data, 4)
1725 0 : do iDim3 = 1, size(data, 3)
1726 0 : temp_data(:, :, size(data, 3) - iDim3 + 1, iDim4) = data(:, :, iDim3, iDim4)
1727 : end do
1728 : end do
1729 0 : else if (iDim == 4_i4) then
1730 0 : do iDim4 = 1, size(data, 4)
1731 0 : temp_data(:, :, :, size(data, 4) - iDim4 + 1) = data(:, :, :, iDim4)
1732 : end do
1733 : end if
1734 0 : call move_alloc(temp_data, data)
1735 0 : end subroutine flip_4D_i4
1736 :
1737 13 : function unpack_chunkwise_dp(vector, mask, field, chunksizeArg) result(unpacked)
1738 : !< this is a chunkwise application of the intrinsic unpack function
1739 : !! it became necessary as the unpack intrinsic can handle only arrays
1740 : !! with size smaller than huge(default_integer_kind)...
1741 : !! it has the following restrictions:\n
1742 : !! - vector must be of type dp
1743 : !! - mask must have rank 1
1744 : !! - field must be a scalar
1745 : real(dp), dimension(:), intent(in) :: vector
1746 : logical, dimension(:), intent(in) :: mask
1747 : real(dp), intent(in) :: field
1748 : real(dp), dimension(size(mask, kind=i8)) :: unpacked
1749 : integer(i8), intent(in), optional :: chunksizeArg
1750 :
1751 : integer(i8) :: i, chunksize, indexMin, indexMax, currentCounts, counts
1752 :
1753 1 : if (present(chunksizeArg)) then
1754 1 : chunksize = chunksizeArg
1755 : else
1756 : chunksize = int(huge(0_i4), i8)
1757 : end if
1758 : ! init some values
1759 1 : i = 1_i8
1760 1 : indexMax = i * chunksize
1761 1 : currentCounts = 1_i8
1762 6 : do while (indexMax < size(mask, kind=i8))
1763 : ! get the indices for the mask
1764 5 : indexMin = (i-1) * chunksize + 1_i8
1765 15 : indexMax = minval([i * chunksize, size(mask, kind=i8)])
1766 : ! this is the indexer for the vector
1767 15 : counts = count(mask(indexMin: indexMax), kind=i8)
1768 : ! unpack slices of maximum size
1769 5 : if (counts == (indexMax - indexMin + 1_i8)) then
1770 6 : unpacked(indexMin: indexMax) = vector(currentCounts: currentCounts + counts - 1_i8)
1771 3 : else if (counts == 0_i8) then
1772 6 : unpacked(indexMin: indexMax) = field
1773 : else
1774 0 : unpacked(indexMin: indexMax) = unpack(vector(currentCounts: currentCounts + counts - 1_i8), &
1775 : mask(indexMin: indexMax), &
1776 1 : field)
1777 : end if
1778 : ! advance the counters
1779 5 : currentCounts = currentCounts + counts
1780 5 : i = i + 1_i8
1781 : end do
1782 :
1783 1 : end function unpack_chunkwise_dp
1784 :
1785 0 : function unpack_chunkwise_i1(vector, mask, field, chunksizeArg) result(unpacked)
1786 : !< this is a chunkwise application of the intrinsic unpack function
1787 : !! it has the following restrictions:\n
1788 : !! - vector must be of type i1
1789 : !! - mask must have rank 1
1790 : !! - field must be a scalar
1791 : integer(i1), dimension(:), intent(in) :: vector
1792 : logical, dimension(:), intent(in) :: mask
1793 : integer(i1), intent(in) :: field
1794 : integer(i1), dimension(size(mask, kind=i8)) :: unpacked
1795 : integer(i8), intent(in), optional :: chunksizeArg
1796 :
1797 : integer(i8) :: i, chunksize, indexMin, indexMax, currentCounts, counts
1798 :
1799 0 : if (present(chunksizeArg)) then
1800 0 : chunksize = chunksizeArg
1801 : else
1802 : chunksize = int(huge(0_i4), i8)
1803 : end if
1804 : ! init some values
1805 0 : i = 1_i8
1806 0 : indexMax = i * chunksize
1807 0 : currentCounts = 1_i8
1808 0 : do while (indexMax < size(mask, kind=i8))
1809 : ! get the indices for the mask
1810 0 : indexMin = (i-1) * chunksize + 1_i8
1811 0 : indexMax = minval([i * chunksize, size(mask, kind=i8)])
1812 : ! this is the indexer for the vector
1813 0 : counts = count(mask(indexMin: indexMax), kind=i8)
1814 : ! unpack slices of maximum size
1815 0 : unpacked(indexMin: indexMax) = unpack(vector(currentCounts: currentCounts + counts - 1_i8), &
1816 : mask(indexMin: indexMax), &
1817 0 : field)
1818 : ! advance the counters
1819 0 : currentCounts = currentCounts + counts
1820 0 : i = i + 1_i8
1821 : end do
1822 :
1823 1 : end function unpack_chunkwise_i1
1824 :
1825 :
1826 : END MODULE mo_utils
|
