Line data Source code
1 : ! ******************************************************************
2 : ! ------------------------------------------------------------------
3 : MODULE MESSY_NCREGRID_BASE
4 : ! ------------------------------------------------------------------
5 : ! Author: Patrick Joeckel, MPICH, Mainz, June 2002
6 : ! ******************************************************************
7 :
8 : !#if defined(MESSY)
9 : ! USE messy_main_constants_mem, ONLY: SP, DP, I4, I8
10 : !#else
11 : ! USE typeSizes, ONLY: SP => FourByteReal &
12 : ! , DP => EightByteReal &
13 : ! , I4 => FourByteInt &
14 : ! , I8 => EightByteInt
15 : !#endif
16 : USE HCO_ERROR_MOD, ONLY : SP, DP, I4, I8
17 :
18 : IMPLICIT NONE
19 :
20 : INTRINSIC :: ASSOCIATED, PRESENT, PRODUCT, REAL, SIZE &
21 : , INT, ABS, DBLE, MAX, MIN, SIGN, TRIM &
22 : , MAXVAL, MINVAL, IAND, NULL
23 :
24 : !GanLuo+ PRIVATE :: ASSOCIATED, PRESENT, PRODUCT, REAL, SIZE &
25 : !GanLuo+ , INT, ABS, DBLE, MAX, MIN, SIGN, TRIM &
26 : !GanLuo+ , MAXVAL, MINVAL, IAND, NULL
27 :
28 : ! REGRIDDING VARIABLE TYPES
29 : INTEGER, PARAMETER :: RG_INT = 1
30 : INTEGER, PARAMETER :: RG_EXT = 2
31 : INTEGER, PARAMETER :: RG_IDX = 3
32 : INTEGER, PARAMETER :: RG_IXF = 4
33 :
34 : ! REGRIDDING DATA TYPES
35 : INTEGER, PARAMETER :: VTYPE_UNDEF = 0
36 : INTEGER, PARAMETER :: VTYPE_INT = 1
37 : INTEGER, PARAMETER :: VTYPE_REAL = 2
38 : INTEGER, PARAMETER :: VTYPE_DOUBLE = 3
39 : INTEGER, PARAMETER :: VTYPE_BYTE = 4
40 : INTEGER, PARAMETER :: VTYPE_CHAR = 5
41 :
42 : ! MESSAGE LEVEL TYPES
43 : INTEGER, PARAMETER :: RGMLE = 0 ! ERROR
44 : INTEGER, PARAMETER :: RGMLEC = 1 ! ERROR CONTINUED
45 : INTEGER, PARAMETER :: RGMLVL = 2 ! LITTLE VERBOSE
46 : INTEGER, PARAMETER :: RGMLVLC = 3 ! LITTLE VERBOSE CONTINUED
47 : INTEGER, PARAMETER :: RGMLW = 4 ! WARNING
48 : INTEGER, PARAMETER :: RGMLWC = 5 ! WARNING CONTINUED
49 : INTEGER, PARAMETER :: RGMLVM = 6 ! MEDIUM VERBOSE
50 : INTEGER, PARAMETER :: RGMLVMC = 7 ! MEDIUM VERBOSE CONTINUED
51 : INTEGER, PARAMETER :: RGMLI = 8 ! INFO
52 : INTEGER, PARAMETER :: RGMLIC = 9 ! INFO CONTINUED
53 :
54 : ! MESSAGE OUTPUT LEVEL
55 : INTEGER, PARAMETER :: MSGMODE_S = 0 ! SILENT
56 : INTEGER, PARAMETER :: MSGMODE_E = 1 ! ERROR MESSAGES
57 : INTEGER, PARAMETER :: MSGMODE_VL = 2 ! LITTLE VERBOSE
58 : INTEGER, PARAMETER :: MSGMODE_W = 4 ! WARNING MESSAGES
59 : INTEGER, PARAMETER :: MSGMODE_VM = 8 ! MEDIUM VERBOSE
60 : INTEGER, PARAMETER :: MSGMODE_I = 16 ! INFO MESSAGES
61 : ! INTEGER, SAVE :: MSGMODE = MSGMODE_S + MSGMODE_E + MSGMODE_VL &
62 : ! + MSGMODE_W + MSGMODE_VM + MSGMODE_I
63 : INTEGER, SAVE :: MSGMODE = 1
64 :
65 : TYPE narray
66 : ! n-dimenional array as 1D (LINEAR) array (REAL)
67 : INTEGER (I4) :: n = 0 ! number of dimensions
68 : INTEGER (I4), DIMENSION(:), POINTER :: dim => NULL()! dim. vector
69 : REAL (SP) , DIMENSION(:), POINTER :: vr => NULL() ! real values
70 : REAL (DP) , DIMENSION(:), POINTER :: vd => NULL() ! double values
71 : INTEGER (I8), DIMENSION(:), POINTER :: vi => NULL() ! integer values
72 : INTEGER (I4), DIMENSION(:), POINTER :: vb => NULL() ! byte values
73 : CHARACTER, DIMENSION(:), POINTER :: vc => NULL() ! char. values
74 : END TYPE narray
75 :
76 : TYPE axis
77 : ! hyper-axis (for curvilinear coordinates)
78 : TYPE (narray) :: dat ! interface bounds
79 : LOGICAL :: lm = .false. ! modulo axis ?
80 : INTEGER (I4) :: ndp = 0 ! number of dependencies
81 : ! FIRST IN LIST MUST BE INDEPENDENT !!!
82 : ! E.G. IF DIM 3 DEPENDS ON DIM 1,2, and 5
83 : ! dep = (/ 3,1,2,5 /)
84 : INTEGER (I4), DIMENSION(:), POINTER :: dep => NULL() ! dependencies
85 : END TYPE axis
86 :
87 : !! NOTE: DOES NOT WORK PROPERLY FOR SOME COMPILERS ...
88 : ! INTERFACE ASSIGNMENT (=)
89 : ! MODULE PROCEDURE COPY_NARRAY
90 : ! MODULE PROCEDURE COPY_AXIS
91 : ! END INTERFACE
92 :
93 : INTERFACE QSORT
94 : ! THE GOOD OLD (RECURSIVE) QUICKSORT ALGORITHM FOR LINEAR ARRAYS
95 : MODULE PROCEDURE QSORT_I ! INTEGER
96 : #if !(defined(__SX__))
97 : MODULE PROCEDURE QSORT_B ! BYTE
98 : MODULE PROCEDURE QSORT_R ! REAL
99 : #endif
100 : MODULE PROCEDURE QSORT_D ! DOUBLE PRECISION
101 : END INTERFACE
102 :
103 : INTERFACE OVL
104 : ! CALCULATES THE OVERLAP BETWEEN TWO INTERVALS
105 : ! (AS FRACTION OF THE TWO INTERVALS)
106 : ! > ALSO APPLICABLE TO 'MODULO' INTERVALS
107 : #if !(defined(__SX__))
108 : MODULE PROCEDURE OVL_RR ! REAL - REAL
109 : MODULE PROCEDURE OVL_RD ! REAL - DOUBLE PRECISION
110 : MODULE PROCEDURE OVL_DR ! DOUBLE PRECISION - REAL
111 : #endif
112 : MODULE PROCEDURE OVL_DD ! DOUBLE PRECISION - DOUBLE PRECISION
113 : END INTERFACE
114 :
115 : INTERFACE OVL_1D
116 : ! CALCULATES THE OVERLAP MATRICES BETWEEN TWO CONTINUOUS
117 : ! INTERVAL SEQUENCES (IN UNITS OF INTERVAL FRACTIONS)
118 : ! > ALSO APPLICABLE TO 'MODULO' SEQUENCES
119 : #if !(defined(__SX__))
120 : MODULE PROCEDURE OVL_1D_RR ! REAL - REAL
121 : MODULE PROCEDURE OVL_1D_RD ! REAL - DOUBLE PRECISION
122 : MODULE PROCEDURE OVL_1D_DR ! DOUBLE PRECISION - REAL
123 : #endif
124 : MODULE PROCEDURE OVL_1D_DD ! DOUBLE PRECISION - DOUBLE PRECISION
125 : END INTERFACE
126 :
127 : INTERFACE RGMSG ! MESSAGE OUTPUT
128 : MODULE PROCEDURE RGMSG_C
129 : MODULE PROCEDURE RGMSG_I
130 : MODULE PROCEDURE RGMSG_IA
131 : MODULE PROCEDURE RGMSG_R
132 : END INTERFACE
133 :
134 : CONTAINS
135 :
136 : ! ------------------------------------------------------------------
137 0 : INTEGER FUNCTION QTYPE_NARRAY(na)
138 :
139 : IMPLICIT NONE
140 :
141 : ! I/O
142 : TYPE (narray), INTENT(IN) :: na
143 :
144 0 : QTYPE_NARRAY = VTYPE_UNDEF
145 0 : IF (ASSOCIATED(na%vi)) QTYPE_NARRAY = VTYPE_INT
146 0 : IF (ASSOCIATED(na%vb)) QTYPE_NARRAY = VTYPE_BYTE
147 0 : IF (ASSOCIATED(na%vc)) QTYPE_NARRAY = VTYPE_CHAR
148 0 : IF (ASSOCIATED(na%vr)) QTYPE_NARRAY = VTYPE_REAL
149 0 : IF (ASSOCIATED(na%vd)) QTYPE_NARRAY = VTYPE_DOUBLE
150 :
151 0 : END FUNCTION QTYPE_NARRAY
152 : ! ------------------------------------------------------------------
153 :
154 : ! ------------------------------------------------------------------
155 0 : SUBROUTINE INIT_NARRAY(na, n, dim, qtype)
156 :
157 : IMPLICIT NONE
158 :
159 : ! I/O
160 : TYPE (narray), INTENT(INOUT) :: na
161 : INTEGER, INTENT(IN), OPTIONAL :: n
162 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: dim
163 : INTEGER, OPTIONAL :: qtype
164 :
165 : ! LOCAL
166 : CHARACTER(LEN=*), PARAMETER :: substr = 'INIT_NARRAY'
167 : INTEGER :: status
168 : INTEGER :: len
169 :
170 0 : na%n = 0
171 0 : IF (ASSOCIATED(na%dim)) THEN
172 0 : IF (SIZE(na%dim) > 0) THEN
173 0 : DEALLOCATE(na%dim, STAT=status)
174 0 : CALL ERRMSG(substr,status,1)
175 : END IF
176 0 : NULLIFY(na%dim)
177 : END IF
178 :
179 0 : IF (ASSOCIATED(na%vi)) THEN
180 0 : IF (SIZE(na%vi) > 0) THEN
181 0 : DEALLOCATE(na%vi, STAT=status)
182 0 : CALL ERRMSG(substr,status,2)
183 : END IF
184 0 : NULLIFY(na%vi)
185 : END IF
186 :
187 0 : IF (ASSOCIATED(na%vr)) THEN
188 0 : IF (SIZE(na%vr) > 0) THEN
189 0 : DEALLOCATE(na%vr, STAT=status)
190 0 : CALL ERRMSG(substr,status,3)
191 : END IF
192 0 : NULLIFY(na%vr)
193 : END IF
194 :
195 0 : IF (ASSOCIATED(na%vd)) THEN
196 0 : IF (SIZE(na%vd) > 0) THEN
197 0 : DEALLOCATE(na%vd, STAT=status)
198 0 : CALL ERRMSG(substr,status,4)
199 : END IF
200 0 : NULLIFY(na%vd)
201 : END IF
202 :
203 0 : IF (ASSOCIATED(na%vc)) THEN
204 0 : IF (SIZE(na%vc) > 0) THEN
205 0 : DEALLOCATE(na%vc, STAT=status)
206 0 : CALL ERRMSG(substr,status,5)
207 : END IF
208 0 : NULLIFY(na%vc)
209 : END IF
210 :
211 0 : IF (ASSOCIATED(na%vb)) THEN
212 0 : IF (SIZE(na%vb) > 0) THEN
213 0 : DEALLOCATE(na%vb, STAT=status)
214 0 : CALL ERRMSG(substr,status,6)
215 : END IF
216 0 : NULLIFY(na%vb)
217 : END IF
218 :
219 0 : IF (PRESENT(dim) .AND. PRESENT(n)) THEN
220 0 : IF (SIZE(dim) /= n) THEN
221 : CALL RGMSG(substr, RGMLE, &
222 0 : 'NUMBER OF DIMENSIONS MISMATCH !' )
223 : END IF
224 : END IF
225 :
226 0 : len = 0
227 0 : IF (PRESENT(n)) THEN
228 0 : na%n = n
229 0 : ALLOCATE(na%dim(n), STAT=status)
230 0 : CALL ERRMSG(substr,status,7)
231 0 : na%dim(:) = 0
232 : END IF
233 :
234 0 : IF (PRESENT(dim)) THEN
235 0 : IF (.NOT. ASSOCIATED(na%dim)) THEN
236 0 : na%n = SIZE(dim)
237 0 : ALLOCATE(na%dim(na%n), STAT=status)
238 0 : CALL ERRMSG(substr,status,8)
239 : END IF
240 0 : na%dim(:) = dim(:)
241 0 : len = PRODUCT(na%dim(:))
242 : END IF
243 :
244 0 : IF (PRESENT(qtype).AND.(len > 0)) THEN
245 0 : SELECT CASE(qtype)
246 : CASE(VTYPE_INT)
247 0 : ALLOCATE(na%vi(len), STAT=status)
248 0 : CALL ERRMSG(substr,status,9)
249 0 : na%vi(:) = 0
250 : CASE(VTYPE_REAL)
251 0 : ALLOCATE(na%vr(len), STAT=status)
252 0 : CALL ERRMSG(substr,status,10)
253 0 : na%vr(:) = 0.0
254 : CASE(VTYPE_DOUBLE)
255 0 : ALLOCATE(na%vd(len), STAT=status)
256 0 : CALL ERRMSG(substr,status,11)
257 0 : na%vd(:) = REAL(0.0, DP)
258 : CASE(VTYPE_CHAR)
259 0 : ALLOCATE(na%vc(len), STAT=status)
260 0 : CALL ERRMSG(substr,status,12)
261 0 : na%vc(:) = ' '
262 : CASE(VTYPE_BYTE)
263 0 : ALLOCATE(na%vb(len), STAT=status)
264 0 : CALL ERRMSG(substr,status,13)
265 0 : na%vi(:) = 0
266 : CASE(VTYPE_UNDEF)
267 : ! DO NOTHING, KEEP UNDEFINED
268 : CASE DEFAULT
269 : CALL RGMSG(substr, RGMLE, &
270 0 : 'REQUESTED TYPE FOR N-ARRAY IS UNRECOGNIZED !' )
271 : END SELECT
272 : END IF
273 :
274 0 : END SUBROUTINE INIT_NARRAY
275 : ! ------------------------------------------------------------------
276 :
277 : ! ------------------------------------------------------------------
278 0 : SUBROUTINE INIT_AXIS(a)
279 :
280 : IMPLICIT NONE
281 :
282 : ! I/O
283 : TYPE (axis), INTENT(INOUT) :: a
284 :
285 : ! LOCAL
286 : CHARACTER(LEN=*), PARAMETER :: substr = 'INIT_AXIS'
287 : INTEGER :: status
288 :
289 :
290 0 : CALL INIT_NARRAY(a%dat)
291 0 : a%lm = .false.
292 0 : a%ndp = 0
293 0 : IF (ASSOCIATED(a%dep)) THEN
294 0 : IF (SIZE(a%dep) > 0) THEN
295 0 : DEALLOCATE(a%dep, STAT=status)
296 0 : CALL ERRMSG(substr,status,1)
297 : END IF
298 0 : NULLIFY(a%dep)
299 : END IF
300 :
301 0 : END SUBROUTINE INIT_AXIS
302 : ! ------------------------------------------------------------------
303 :
304 : ! ------------------------------------------------------------------
305 0 : SUBROUTINE COPY_NARRAY(d, s)
306 :
307 : IMPLICIT NONE
308 :
309 : ! I/O
310 : TYPE (narray), INTENT(OUT) :: d
311 : TYPE (narray), INTENT(IN) :: s
312 :
313 : ! LOCAL
314 : CHARACTER(LEN=*), PARAMETER :: substr = 'COPY_NARRAY'
315 : INTEGER :: n
316 : INTEGER :: vtype
317 : INTEGER :: status
318 :
319 : ! INIT
320 0 : n = 0
321 :
322 0 : d%n = s%n
323 0 : IF (ASSOCIATED(s%dim)) THEN
324 0 : ALLOCATE(d%dim(d%n),STAT=status)
325 0 : CALL ERRMSG(substr,status,1)
326 0 : d%dim(:) = s%dim(:)
327 : END IF
328 :
329 0 : vtype = QTYPE_NARRAY(s)
330 0 : SELECT CASE(vtype)
331 : CASE(VTYPE_INT)
332 0 : n = SIZE(s%vi)
333 0 : IF (n > 0) THEN
334 0 : ALLOCATE(d%vi(n),STAT=status)
335 0 : CALL ERRMSG(substr,status,2)
336 0 : d%vi(:) = s%vi(:)
337 : END IF
338 : CASE(VTYPE_REAL)
339 0 : n = SIZE(s%vr)
340 0 : IF (n > 0) THEN
341 0 : ALLOCATE(d%vr(n),STAT=status)
342 0 : CALL ERRMSG(substr,status,3)
343 0 : d%vr(:) = s%vr(:)
344 : END IF
345 : CASE(VTYPE_DOUBLE)
346 0 : n = SIZE(s%vd)
347 0 : IF (n > 0) THEN
348 0 : ALLOCATE(d%vd(n),STAT=status)
349 0 : CALL ERRMSG(substr,status,4)
350 0 : d%vd(:) = s%vd(:)
351 : END IF
352 : CASE(VTYPE_CHAR)
353 0 : n = SIZE(s%vc)
354 0 : IF (n > 0) THEN
355 0 : ALLOCATE(d%vc(n),STAT=status)
356 0 : CALL ERRMSG(substr,status,5)
357 0 : d%vc(:) = s%vc(:)
358 : END IF
359 : CASE(VTYPE_BYTE)
360 0 : n = SIZE(s%vb)
361 0 : IF (n > 0) THEN
362 0 : ALLOCATE(d%vb(n),STAT=status)
363 0 : CALL ERRMSG(substr,status,6)
364 0 : d%vb(:) = s%vb(:)
365 : END IF
366 : CASE(VTYPE_UNDEF)
367 : ! DO NOTHING, EMPTY N-ARRAY IS COPIED
368 : CASE DEFAULT
369 : CALL RGMSG(substr, RGMLE, &
370 0 : 'N-ARRAY OF UNRECOGNIZED TYPE CANNOT BE COPIED !' )
371 : END SELECT
372 :
373 0 : END SUBROUTINE COPY_NARRAY
374 : ! ------------------------------------------------------------------
375 :
376 : ! ------------------------------------------------------------------
377 0 : SUBROUTINE COPY_AXIS(d, s)
378 :
379 : IMPLICIT NONE
380 :
381 : ! I/O
382 : TYPE (axis), INTENT(OUT) :: d
383 : TYPE (axis), INTENT(IN) :: s
384 :
385 : ! LOCAL
386 : CHARACTER(LEN=*), PARAMETER :: substr = 'COPY_AXIS'
387 : INTEGER :: status
388 :
389 0 : CALL COPY_NARRAY(d%dat, s%dat)
390 0 : d%lm = s%lm
391 0 : d%ndp = s%ndp
392 0 : IF (ASSOCIATED(s%dep)) THEN
393 0 : ALLOCATE(d%dep(d%ndp),STAT=status)
394 0 : CALL ERRMSG(substr,status,1)
395 0 : d%dep(:) = s%dep(:)
396 : END IF
397 :
398 0 : END SUBROUTINE COPY_AXIS
399 : ! ------------------------------------------------------------------
400 :
401 : ! ------------------------------------------------------------------
402 0 : SUBROUTINE SORT_NARRAY(na, nx, reverse)
403 :
404 : IMPLICIT NONE
405 :
406 : ! I/O
407 : TYPE (narray), INTENT(INOUT) :: na
408 : TYPE (narray), INTENT(INOUT) :: nx
409 : LOGICAL , INTENT(IN) ,OPTIONAL :: reverse
410 :
411 : ! LOCAL
412 : CHARACTER(LEN=*), PARAMETER :: substr = 'SORT_NARRAY'
413 : INTEGER :: vtype
414 : LOGICAL :: lrev
415 : INTEGER :: n, i
416 : TYPE (narray) :: nh
417 :
418 0 : IF (PRESENT(reverse)) THEN
419 0 : lrev = reverse
420 : ELSE
421 : lrev = .false. ! DEFAULT
422 : END IF
423 :
424 0 : IF (na%n == 0) THEN
425 0 : CALL RGMSG(substr, RGMLW, 'EMPTY ARRAY ! NOTHING TO DO !')
426 0 : RETURN
427 : END IF
428 :
429 0 : IF (na%n > 1) THEN
430 : CALL RGMSG(substr, RGMLW, &
431 0 : 'SORTING A ',na%n,'-DIMENSIONAL N-ARRAY AS LINEAR ARRAY !')
432 : END IF
433 :
434 0 : IF (lrev) THEN
435 :
436 0 : vtype = QTYPE_NARRAY(nx)
437 0 : IF (vtype /= VTYPE_INT) THEN
438 0 : CALL RGMSG(substr, RGMLE, 'INDEX N-ARRAY MUST BE OF TYPE INTEGER !')
439 : END IF
440 0 : n = SIZE(nx%vi)
441 0 : CALL COPY_NARRAY(nh, na) ! COPY TO BE RE-ORDERED
442 0 : vtype = QTYPE_NARRAY(na)
443 : SELECT CASE(vtype)
444 : CASE(VTYPE_REAL)
445 0 : DO i=1,n
446 0 : na%vr(nx%vi(i)) = nh%vr(i)
447 : END DO
448 : CASE(VTYPE_DOUBLE)
449 0 : DO i=1,n
450 0 : na%vd(nx%vi(i)) = nh%vd(i)
451 : END DO
452 : CASE(VTYPE_INT)
453 0 : DO i=1,n
454 0 : na%vi(nx%vi(i)) = nh%vi(i)
455 : END DO
456 : CASE(VTYPE_BYTE)
457 0 : DO i=1,n
458 0 : na%vb(nx%vi(i)) = nh%vb(i)
459 : END DO
460 : CASE(VTYPE_CHAR)
461 0 : CALL RGMSG(substr,RGMLE,'UN-SORTING OF TYPE CHAR IS NOT IMPLEMENTED !')
462 : CASE(VTYPE_UNDEF)
463 0 : CALL RGMSG(substr,RGMLE,'ARRAY OF UNDEFINED TYPE CANNOT BE UN-SORTED !')
464 : CASE DEFAULT
465 0 : CALL RGMSG(substr,RGMLE,'ARRAY OF UNRECOGNIZED TYPE CANNOT BE UN-SORTED !')
466 : END SELECT
467 : ! CLEAN UP
468 0 : CALL INIT_NARRAY(nh)
469 :
470 : ELSE
471 :
472 0 : CALL INIT_NARRAY(nx, na%n, na%dim)
473 : !
474 0 : vtype = QTYPE_NARRAY(na)
475 0 : SELECT CASE(vtype)
476 : CASE(VTYPE_REAL)
477 0 : CALL QSORT_R(na%vr, nx%vi)
478 : CASE(VTYPE_DOUBLE)
479 0 : CALL QSORT_D(na%vd, nx%vi)
480 : CASE(VTYPE_INT)
481 0 : CALL QSORT_I(na%vi, nx%vi)
482 : CASE(VTYPE_BYTE)
483 0 : CALL QSORT_B(na%vb, nx%vi)
484 : CASE(VTYPE_CHAR)
485 0 : CALL RGMSG(substr,RGMLE,'SORTING OF TYPE CHAR IS NOT IMPLEMENTED !')
486 : CASE(VTYPE_UNDEF)
487 0 : CALL RGMSG(substr,RGMLE,'ARRAY OF UNDEFINED TYPE CANNOT BE SORTED !')
488 : CASE DEFAULT
489 0 : CALL RGMSG(substr,RGMLE,'ARRAY OF UNRECOGNIZED TYPE CANNOT BE SORTED !')
490 : END SELECT
491 :
492 : END IF ! (reverse ?)
493 :
494 0 : END SUBROUTINE SORT_NARRAY
495 : ! ------------------------------------------------------------------
496 :
497 : ! ------------------------------------------------------------------
498 0 : SUBROUTINE REORDER_NARRAY(na, nx)
499 :
500 : IMPLICIT NONE
501 :
502 : ! I/O
503 : TYPE (narray), INTENT(INOUT) :: na ! n-array to reorder
504 : TYPE (narray), INTENT(IN) :: nx ! index n-array
505 :
506 : ! LOCAL
507 : CHARACTER(LEN=*), PARAMETER :: substr = 'REORDER_NARRAY'
508 : TYPE (narray) :: nh ! copy of na
509 : INTEGER :: vtype
510 : INTEGER :: i, n
511 :
512 0 : IF (na%n == 0) THEN
513 0 : CALL RGMSG(substr, RGMLW, 'EMPTY ARRAY ! NOTHING TO DO !')
514 0 : RETURN
515 : END IF
516 :
517 0 : vtype = QTYPE_NARRAY(nx)
518 0 : IF (vtype /= VTYPE_INT) THEN
519 0 : CALL RGMSG(substr, RGMLE, 'INDEX N-ARRAY MUST BE OF TYPE INT !')
520 : END IF
521 :
522 0 : IF (na%n > 1) THEN
523 : CALL RGMSG(substr, RGMLW, 'REORDERING A ',na%n, &
524 0 : '-DIMENSIONAL N-ARRAY AS LINEAR ARRAY !')
525 : END IF
526 :
527 0 : IF (na%n /= nx%n) THEN
528 : CALL RGMSG(substr, RGMLE, 'DIMENSION MISMATCH BETWEEN N-ARRAY (',&
529 0 : na%n,')' , .false.)
530 0 : CALL RGMSG(substr, RGMLEC, 'AND INDEX N-ARRAY (',nx%n,') !')
531 : END IF
532 :
533 0 : IF (PRODUCT(na%dim) /= PRODUCT(nx%dim)) THEN
534 0 : CALL RGMSG(substr, RGMLE, 'LENGTH OF N-ARRAY MISMATCH BETWEEN', .false.)
535 0 : CALL RGMSG(substr, RGMLEC, 'N-ARRAY (',PRODUCT(na%dim),') AND', .false.)
536 0 : CALL RGMSG(substr, RGMLEC, 'INDEX N-ARRAY (',PRODUCT(nx%dim),') !')
537 : END IF
538 :
539 0 : vtype = QTYPE_NARRAY(na)
540 0 : CALL COPY_NARRAY(nh, na)
541 0 : n = PRODUCT(na%dim)
542 :
543 : SELECT CASE(vtype)
544 : CASE(VTYPE_REAL)
545 0 : DO i=1,n
546 0 : na%vr(i) = nh%vr(nx%vi(i))
547 : END DO
548 : CASE(VTYPE_DOUBLE)
549 0 : DO i=1,n
550 0 : na%vd(i) = nh%vd(nx%vi(i))
551 : END DO
552 : CASE(VTYPE_INT)
553 0 : DO i=1,n
554 0 : na%vi(i) = nh%vi(nx%vi(i))
555 : END DO
556 : CASE(VTYPE_BYTE)
557 0 : DO i=1,n
558 0 : na%vb(i) = nh%vb(nx%vi(i))
559 : END DO
560 : CASE(VTYPE_CHAR)
561 0 : DO i=1,n
562 0 : na%vc(i) = nh%vc(nx%vi(i))
563 : END DO
564 : CASE(VTYPE_UNDEF)
565 0 : CALL RGMSG(substr, RGMLE, 'ARRAY OF UNDEFINED TYPE CANNOT BE UN-SORTED !')
566 : CASE DEFAULT
567 0 : CALL RGMSG(substr, RGMLE, 'ARRAY OF UNRECOGNIZED TYPE CANNOT BE UN-SORTED !')
568 : END SELECT
569 :
570 : ! CLEAN UP
571 0 : CALL INIT_NARRAY(nh)
572 :
573 0 : END SUBROUTINE REORDER_NARRAY
574 : ! ------------------------------------------------------------------
575 :
576 : ! ------------------------------------------------------------------
577 0 : SUBROUTINE DOUBLE_NARRAY(na)
578 :
579 : IMPLICIT NONE
580 :
581 : ! I/O
582 : TYPE (narray), INTENT(INOUT) :: na
583 :
584 : ! LOCAL
585 : CHARACTER(LEN=*), PARAMETER :: substr = 'DOUBLE_NARRAY'
586 : INTEGER :: vtype
587 : INTEGER :: status
588 :
589 0 : vtype = QTYPE_NARRAY(na)
590 0 : SELECT CASE(vtype)
591 : CASE(VTYPE_REAL)
592 0 : ALLOCATE(na%vd(SIZE(na%vr)), STAT=status)
593 0 : CALL ERRMSG(substr,status,1)
594 0 : na%vd(:) = REAL(na%vr(:), DP)
595 0 : DEALLOCATE(na%vr, STAT=status)
596 0 : CALL ERRMSG(substr,status,2)
597 0 : NULLIFY(na%vr)
598 : CASE(VTYPE_DOUBLE)
599 : ! NOTHING TO DO
600 : CASE(VTYPE_INT)
601 0 : ALLOCATE(na%vd(SIZE(na%vi)), STAT=status)
602 0 : CALL ERRMSG(substr,status,3)
603 0 : na%vd(:) = REAL(na%vi(:), DP)
604 0 : DEALLOCATE(na%vi, STAT=status)
605 0 : CALL ERRMSG(substr,status,4)
606 0 : NULLIFY(na%vi)
607 : CASE(VTYPE_BYTE)
608 0 : ALLOCATE(na%vd(SIZE(na%vi)), STAT=status)
609 0 : CALL ERRMSG(substr,status,5)
610 0 : na%vd(:) = REAL(na%vb(:), DP)
611 0 : DEALLOCATE(na%vb, STAT=status)
612 0 : CALL ERRMSG(substr,status,6)
613 0 : NULLIFY(na%vb)
614 : CASE(VTYPE_CHAR)
615 0 : CALL RGMSG(substr, RGMLE, 'CHAR CANNOT BE CONVERTED TO DOUBLEPRECISION !')
616 : CASE(VTYPE_UNDEF)
617 0 : CALL RGMSG(substr, RGMLE, 'UNDEFINED N-ARRAY TYPE !')
618 : CASE DEFAULT
619 0 : CALL RGMSG(substr, RGMLE, 'UNRECOGNIZED N-ARRAY TYPE !')
620 : END SELECT
621 :
622 0 : END SUBROUTINE DOUBLE_NARRAY
623 : ! ------------------------------------------------------------------
624 :
625 : ! ------------------------------------------------------------------
626 0 : SUBROUTINE SCALE_NARRAY(na, sc)
627 :
628 : IMPLICIT NONE
629 :
630 : ! I/O
631 : TYPE (narray), INTENT(INOUT) :: na ! N-array
632 : REAL , INTENT(IN) :: sc ! scaling factor
633 :
634 : ! LOCAL
635 : CHARACTER(LEN=*), PARAMETER :: substr = 'SCALE_NARRAY'
636 : INTEGER :: vtype
637 : INTEGER :: status
638 : INTEGER :: i
639 :
640 0 : vtype = QTYPE_NARRAY(na)
641 :
642 0 : SELECT CASE(vtype)
643 : CASE(VTYPE_REAL)
644 0 : na%vr(:) = na%vr(:) * REAL(sc, SP)
645 : CASE(VTYPE_DOUBLE)
646 0 : na%vd(:) = na%vd(:) * REAL(sc, DP)
647 : CASE(VTYPE_INT)
648 0 : CALL RGMSG(substr, RGMLI, 'N-ARRAY OF TYPE INTEGER CONVERTED TO REAL !')
649 0 : ALLOCATE(na%vr(SIZE(na%vi)), STAT=status)
650 0 : CALL ERRMSG(substr,status,1)
651 0 : DO i=1, SIZE(na%vi)
652 0 : na%vr(i) = REAL(na%vi(i), SP) * REAL(sc, SP)
653 : END DO
654 0 : DEALLOCATE(na%vi, STAT=status)
655 0 : CALL ERRMSG(substr,status,2)
656 0 : NULLIFY(na%vi)
657 : CASE(VTYPE_BYTE)
658 0 : CALL RGMSG(substr, RGMLI, 'N-ARRAY OF TYPE BYTE CONVERTED TO REAL !')
659 0 : ALLOCATE(na%vr(SIZE(na%vb)), STAT=status)
660 0 : CALL ERRMSG(substr,status,3)
661 0 : DO i=1, SIZE(na%vb)
662 0 : na%vr(i) = REAL(INT(na%vb(i), I8), SP) * REAL(sc, SP)
663 : END DO
664 0 : DEALLOCATE(na%vb, STAT=status)
665 0 : CALL ERRMSG(substr,status,4)
666 0 : NULLIFY(na%vb)
667 : CASE(VTYPE_CHAR)
668 0 : CALL RGMSG(substr, RGMLE, 'N-ARRAY OF TYPE CHAR CANNOT BE SCALED !')
669 : CASE(VTYPE_UNDEF)
670 0 : CALL RGMSG(substr, RGMLE, 'UNDEFINED N-ARRAY TYPE !')
671 : CASE DEFAULT
672 0 : CALL RGMSG(substr, RGMLE, 'UNRECOGNIZED N-ARRAY TYPE !')
673 : END SELECT
674 :
675 0 : END SUBROUTINE SCALE_NARRAY
676 : ! ------------------------------------------------------------------
677 :
678 : ! ------------------------------------------------------------------
679 0 : SUBROUTINE CAT_NARRAY(na, nb)
680 :
681 : IMPLICIT NONE
682 :
683 : ! I/O
684 : TYPE(narray), INTENT(INOUT) :: na
685 : TYPE(narray), INTENT(in) :: nb
686 :
687 : ! LOCAL
688 : CHARACTER(LEN=*), PARAMETER :: substr = 'CAT_NARRAY'
689 : TYPE(narray) :: nh
690 : INTEGER :: vtype1, vtype2
691 : INTEGER :: n,m,i
692 :
693 0 : vtype1 = QTYPE_NARRAY(na)
694 0 : vtype2 = QTYPE_NARRAY(nb)
695 :
696 0 : IF (vtype2 == VTYPE_UNDEF) THEN
697 0 : CALL RGMSG(substr, RGMLE, 'N-ARRAY TO BE APPENDED MUST BE DEFINED !')
698 : ELSE
699 0 : IF (vtype1 == VTYPE_UNDEF) THEN
700 0 : CALL COPY_NARRAY(na, nb)
701 0 : RETURN
702 : ELSE
703 0 : IF (vtype1 /= vtype2) THEN
704 0 : CALL RGMSG(substr, RGMLE, 'N-ARRAYS MUST BE OF SAME TYPE !')
705 : END IF
706 : END IF
707 0 : IF ((na%n /= 1).OR.(nb%n /= 1)) THEN
708 0 : CALL RGMSG(substr, RGMLE, 'N-ARRAYS MUST BE 1-DIMENSIONAL !')
709 : END IF
710 : END IF
711 :
712 0 : n = na%dim(1)
713 0 : m = nb%dim(1)
714 0 : CALL INIT_NARRAY(nh, 1, (/ n+m /), vtype2)
715 :
716 : SELECT CASE(vtype2)
717 : CASE(VTYPE_REAL)
718 0 : DO i=1, n
719 0 : nh%vr(i) = na%vr(i)
720 : END DO
721 0 : DO i=1, m
722 0 : nh%vr(n+i) = nb%vr(i)
723 : END DO
724 : CASE(VTYPE_DOUBLE)
725 0 : DO i=1, n
726 0 : nh%vd(i) = na%vd(i)
727 : END DO
728 0 : DO i=1, m
729 0 : nh%vd(n+i) = nb%vd(i)
730 : END DO
731 : CASE(VTYPE_INT)
732 0 : DO i=1, n
733 0 : nh%vi(i) = na%vi(i)
734 : END DO
735 0 : DO i=1, m
736 0 : nh%vi(n+i) = nb%vi(i)
737 : END DO
738 : CASE(VTYPE_BYTE)
739 0 : DO i=1, n
740 0 : nh%vb(i) = na%vb(i)
741 : END DO
742 0 : DO i=1, m
743 0 : nh%vb(n+i) = nb%vb(i)
744 : END DO
745 : CASE(VTYPE_CHAR)
746 0 : DO i=1, n
747 0 : nh%vc(i) = na%vc(i)
748 : END DO
749 0 : DO i=1, m
750 0 : nh%vc(n+i) = nb%vc(i)
751 : END DO
752 : CASE(VTYPE_UNDEF)
753 : ! NOTHING TO DO FOR UNDEFINED ARRAYS
754 : CASE DEFAULT
755 : ! NOTHING TO DO FOR UNRECOGNIZED ARRAYS
756 : END SELECT
757 :
758 0 : CALL COPY_NARRAY(na, nh)
759 : ! CLEAN UP
760 0 : CALL INIT_NARRAY(nh)
761 :
762 0 : END SUBROUTINE CAT_NARRAY
763 : ! ------------------------------------------------------------------
764 :
765 : ! ------------------------------------------------------------------
766 0 : RECURSIVE SUBROUTINE QSORT_I(data,index,ileft,iright)
767 :
768 : IMPLICIT NONE
769 :
770 : ! I/O
771 : INTEGER (I8), DIMENSION(:), INTENT(INOUT) :: data ! data to sort
772 : INTEGER (I8), DIMENSION(:), POINTER :: index ! index list
773 : INTEGER (I8), INTENT(IN), OPTIONAL :: ileft, iright
774 :
775 : ! LOCAL
776 : CHARACTER(LEN=*), PARAMETER :: substr = 'QSORT_I'
777 : INTEGER (I8) :: temp ! temporal data
778 : INTEGER :: n ! LENGTH OF LIST
779 : INTEGER (I8) :: left, right
780 : INTEGER (I8) :: i, last, apu
781 : INTEGER (I8) :: ti ! temporal index
782 : INTEGER :: status
783 :
784 : ! INIT
785 0 : n = SIZE(data)
786 0 : IF (.NOT.ASSOCIATED(index)) THEN
787 0 : ALLOCATE(index(n), STAT=status)
788 0 : CALL ERRMSG(substr,status,1)
789 0 : DO i=1, n
790 0 : index(i) = i
791 : END DO
792 : END IF
793 :
794 0 : IF (PRESENT(ileft)) THEN
795 0 : left = ileft
796 : ELSE
797 0 : left = 1
798 : END IF
799 0 : IF (PRESENT(iright)) THEN
800 0 : right = iright
801 : ELSE
802 0 : right = n
803 : END IF
804 :
805 0 : IF(left >= right) RETURN
806 :
807 0 : apu = (left+right)/2
808 :
809 0 : temp = data(left)
810 0 : data(left) = data(apu)
811 0 : data(apu) = temp
812 :
813 0 : ti = index(left)
814 0 : index(left) = index(apu)
815 0 : index(apu) = ti
816 :
817 0 : last = left
818 :
819 0 : DO i=left+1,right
820 0 : IF(data(i) < data(left)) THEN
821 0 : last = last+1
822 0 : temp = data(last)
823 0 : data(last) = data(i)
824 0 : data(i) = temp
825 :
826 0 : ti = index(last)
827 0 : index(last) = index(i)
828 0 : index(i) = ti
829 : ENDIF
830 : END DO
831 :
832 0 : temp = data(left)
833 0 : data(left) = data(last)
834 0 : data(last) = temp
835 :
836 0 : ti = index(left)
837 0 : index(left) = index(last)
838 0 : index(last) = ti
839 :
840 0 : CALL QSORT_I(data,index,left,last-1)
841 0 : CALL QSORT_I(data,index,last+1,right)
842 :
843 0 : RETURN
844 :
845 : END SUBROUTINE QSORT_I
846 : ! ------------------------------------------------------------------
847 :
848 : ! ------------------------------------------------------------------
849 0 : RECURSIVE SUBROUTINE QSORT_B(data,index,ileft,iright)
850 :
851 : IMPLICIT NONE
852 :
853 : ! I/O
854 : INTEGER (I4), DIMENSION(:), INTENT(INOUT) :: data ! data to sort
855 : INTEGER (I8), DIMENSION(:), POINTER :: index ! index list
856 : INTEGER (I8), INTENT(IN), OPTIONAL :: ileft, iright
857 :
858 : ! LOCAL
859 : CHARACTER(LEN=*), PARAMETER :: substr = 'QSORT_B'
860 : INTEGER (I4) :: temp ! temporal data
861 : INTEGER :: n ! LENGTH OF LIST
862 : INTEGER (I8) :: left, right
863 : INTEGER (I8) :: i, last, apu
864 : INTEGER (I8) :: ti ! temporal index
865 : INTEGER :: status
866 :
867 : ! INIT
868 0 : n = SIZE(data)
869 0 : IF (.NOT.ASSOCIATED(index)) THEN
870 0 : ALLOCATE(index(n), STAT=status)
871 0 : CALL ERRMSG(substr,status,1)
872 0 : DO i=1, n
873 0 : index(i) = i
874 : END DO
875 : END IF
876 :
877 0 : IF (PRESENT(ileft)) THEN
878 0 : left = ileft
879 : ELSE
880 0 : left = 1
881 : END IF
882 0 : IF (PRESENT(iright)) THEN
883 0 : right = iright
884 : ELSE
885 0 : right = n
886 : END IF
887 :
888 0 : IF(left >= right) RETURN
889 :
890 0 : apu = (left+right)/2
891 :
892 0 : temp = data(left)
893 0 : data(left) = data(apu)
894 0 : data(apu) = temp
895 :
896 0 : ti = index(left)
897 0 : index(left) = index(apu)
898 0 : index(apu) = ti
899 :
900 0 : last = left
901 :
902 0 : DO i=left+1,right
903 0 : IF(data(i) < data(left)) THEN
904 0 : last = last+1
905 0 : temp = data(last)
906 0 : data(last) = data(i)
907 0 : data(i) = temp
908 :
909 0 : ti = index(last)
910 0 : index(last) = index(i)
911 0 : index(i) = ti
912 : ENDIF
913 : END DO
914 :
915 0 : temp = data(left)
916 0 : data(left) = data(last)
917 0 : data(last) = temp
918 :
919 0 : ti = index(left)
920 0 : index(left) = index(last)
921 0 : index(last) = ti
922 :
923 0 : CALL QSORT_B(data,index,left,last-1)
924 0 : CALL QSORT_B(data,index,last+1,right)
925 :
926 0 : RETURN
927 :
928 : END SUBROUTINE QSORT_B
929 : ! ------------------------------------------------------------------
930 :
931 : ! ------------------------------------------------------------------
932 0 : RECURSIVE SUBROUTINE QSORT_R(data,index,ileft,iright)
933 :
934 : IMPLICIT NONE
935 :
936 : ! I/O
937 : REAL (SP), DIMENSION(:), INTENT(INOUT) :: data ! data to sort
938 : INTEGER (I8), DIMENSION(:), POINTER :: index ! index list
939 : INTEGER (I8), INTENT(IN), OPTIONAL :: ileft, iright
940 :
941 : ! LOCAL
942 : CHARACTER(LEN=*), PARAMETER :: substr = 'QSORT_R'
943 : REAL (SP) :: temp ! temporal data
944 : INTEGER :: n ! LENGTH OF LIST
945 : INTEGER (I8) :: left, right
946 : INTEGER (I8) :: i, last, apu
947 : INTEGER (I8) :: ti ! temporal index
948 : INTEGER :: status
949 :
950 : ! INIT
951 0 : n = SIZE(data)
952 0 : IF (.NOT.ASSOCIATED(index)) THEN
953 0 : ALLOCATE(index(n),STAT=status)
954 0 : CALL ERRMSG(substr,status,1)
955 0 : DO i=1, n
956 0 : index(i) = i
957 : END DO
958 : END IF
959 :
960 0 : IF (PRESENT(ileft)) THEN
961 0 : left = ileft
962 : ELSE
963 0 : left = 1
964 : END IF
965 0 : IF (PRESENT(iright)) THEN
966 0 : right = iright
967 : ELSE
968 0 : right = n
969 : END IF
970 :
971 0 : IF(left >= right) RETURN
972 :
973 0 : apu = (left+right)/2
974 :
975 0 : temp = data(left)
976 0 : data(left) = data(apu)
977 0 : data(apu) = temp
978 :
979 0 : ti = index(left)
980 0 : index(left) = index(apu)
981 0 : index(apu) = ti
982 :
983 0 : last = left
984 :
985 0 : DO i=left+1,right
986 0 : IF(data(i) < data(left)) THEN
987 0 : last = last+1
988 0 : temp = data(last)
989 0 : data(last) = data(i)
990 0 : data(i) = temp
991 :
992 0 : ti = index(last)
993 0 : index(last) = index(i)
994 0 : index(i) = ti
995 : ENDIF
996 : END DO
997 :
998 0 : temp = data(left)
999 0 : data(left) = data(last)
1000 0 : data(last) = temp
1001 :
1002 0 : ti = index(left)
1003 0 : index(left) = index(last)
1004 0 : index(last) = ti
1005 :
1006 0 : CALL QSORT_R(data,index,left,last-1)
1007 0 : CALL QSORT_R(data,index,last+1,right)
1008 :
1009 0 : RETURN
1010 :
1011 : END SUBROUTINE QSORT_R
1012 : ! ------------------------------------------------------------------
1013 :
1014 : ! ------------------------------------------------------------------
1015 0 : RECURSIVE SUBROUTINE QSORT_D(data,index,ileft,iright)
1016 :
1017 : IMPLICIT NONE
1018 :
1019 : ! I/O
1020 : REAL (DP), DIMENSION(:), INTENT(INOUT) :: data ! data to sort
1021 : INTEGER (I8), DIMENSION(:), POINTER :: index ! index list
1022 : INTEGER (I8), INTENT(IN), OPTIONAL :: ileft, iright
1023 :
1024 : ! LOCAL
1025 : CHARACTER(LEN=*), PARAMETER :: substr = 'QSORT_D'
1026 : REAL (DP) :: temp ! temporal data
1027 : INTEGER :: n ! LENGTH OF LIST
1028 : INTEGER (I8) :: left, right
1029 : INTEGER (I8) :: i, last, apu
1030 : INTEGER (I8) :: ti ! temporal index
1031 : INTEGER :: status
1032 :
1033 : ! INIT
1034 0 : n = SIZE(data)
1035 0 : IF (.NOT.ASSOCIATED(index)) THEN
1036 0 : ALLOCATE(index(n), STAT=status)
1037 0 : CALL ERRMSG(substr,status,1)
1038 0 : DO i=1, n
1039 0 : index(i) = i
1040 : END DO
1041 : END IF
1042 :
1043 0 : IF (PRESENT(ileft)) THEN
1044 0 : left = ileft
1045 : ELSE
1046 0 : left = 1
1047 : END IF
1048 0 : IF (PRESENT(iright)) THEN
1049 0 : right = iright
1050 : ELSE
1051 0 : right = n
1052 : END IF
1053 :
1054 0 : IF(left >= right) RETURN
1055 :
1056 0 : apu = (left+right)/2
1057 :
1058 0 : temp = data(left)
1059 0 : data(left) = data(apu)
1060 0 : data(apu) = temp
1061 :
1062 0 : ti = index(left)
1063 0 : index(left) = index(apu)
1064 0 : index(apu) = ti
1065 :
1066 0 : last = left
1067 :
1068 0 : DO i=left+1,right
1069 0 : IF(data(i) < data(left)) THEN
1070 0 : last = last+1
1071 0 : temp = data(last)
1072 0 : data(last) = data(i)
1073 0 : data(i) = temp
1074 :
1075 0 : ti = index(last)
1076 0 : index(last) = index(i)
1077 0 : index(i) = ti
1078 : ENDIF
1079 : END DO
1080 :
1081 0 : temp = data(left)
1082 0 : data(left) = data(last)
1083 0 : data(last) = temp
1084 :
1085 0 : ti = index(left)
1086 0 : index(left) = index(last)
1087 0 : index(last) = ti
1088 :
1089 0 : CALL QSORT_D(data,index,left,last-1)
1090 0 : CALL QSORT_D(data,index,last+1,right)
1091 :
1092 0 : RETURN
1093 :
1094 : END SUBROUTINE QSORT_D
1095 : ! ------------------------------------------------------------------
1096 :
1097 : ! ------------------------------------------------------------------
1098 0 : SUBROUTINE OVL_RR (sl,sr,dl,dr,fs,fd,modulo)
1099 :
1100 : IMPLICIT NONE
1101 :
1102 : ! I/O
1103 : REAL (SP), INTENT(IN) :: sl, sr ! 'box' edges
1104 : REAL (SP), INTENT(IN) :: dl, dr ! 'box' edges
1105 : REAL (DP), INTENT(IN), OPTIONAL :: modulo ! shift for 'modulo' axis
1106 : REAL (DP), INTENT(OUT) :: fs, fd ! overlap fractions
1107 :
1108 : ! LOCAL
1109 : REAL (DP) :: x1, x2
1110 : REAL (DP) :: shift
1111 : INTEGER :: fx
1112 :
1113 0 : x1 = MAX(MIN(dl,dr),MIN(sl,sr))
1114 0 : x2 = MIN(MAX(dl,dr),MAX(sl,sr))
1115 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1116 0 : fs = (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1117 0 : fd = (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1118 :
1119 0 : IF (PRESENT(modulo)) THEN
1120 :
1121 0 : shift = REAL(INT((dl/modulo) + REAL(1.5, DP)), DP) * modulo
1122 0 : x1 = MAX(DBLE(MIN(dl,dr))+shift,DBLE(MIN(sl,sr)))
1123 0 : x2 = MIN(DBLE(MAX(dl,dr))+shift,DBLE(MAX(sl,sr)))
1124 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1125 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1126 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1127 :
1128 0 : shift = REAL(INT((sl/modulo) + REAL(1.5, DP)), DP) * modulo
1129 0 : x1 = MAX(DBLE(MIN(dl,dr)),DBLE(MIN(sl,sr))+shift)
1130 0 : x2 = MIN(DBLE(MAX(dl,dr)),DBLE(MAX(sl,sr))+shift)
1131 0 : fx = INT(SIGN(REAl(-0.5, DP),x2-x1)+REAL(1., DP))
1132 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1133 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1134 : END IF
1135 :
1136 0 : END SUBROUTINE OVL_RR
1137 : ! ------------------------------------------------------------------
1138 :
1139 : ! ------------------------------------------------------------------
1140 0 : SUBROUTINE OVL_RD (sl,sr,dl,dr,fs,fd,modulo)
1141 :
1142 : IMPLICIT NONE
1143 :
1144 : ! I/O
1145 : REAL (SP), INTENT(IN) :: sl, sr ! 'box' edges
1146 : REAL (DP), INTENT(IN) :: dl, dr ! 'box' edges
1147 : REAL (DP), INTENT(IN), OPTIONAL :: modulo ! shift for 'modulo' axis
1148 : REAL (DP), INTENT(OUT) :: fs, fd ! overlap fractions
1149 :
1150 : ! LOCAL
1151 : REAL (DP) :: x1, x2
1152 : REAL (DP) :: shift
1153 : INTEGER :: fx
1154 :
1155 0 : x1 = MAX(MIN(dl,dr),DBLE(MIN(sl,sr)))
1156 0 : x2 = MIN(MAX(dl,dr),DBLE(MAX(sl,sr)))
1157 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAl(1., DP))
1158 0 : fs = (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1159 0 : fd = (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1160 :
1161 0 : IF (PRESENT(modulo)) THEN
1162 :
1163 0 : shift = REAL(INT((dl/modulo) + REAL(1.5, DP)), DP) * modulo
1164 0 : x1 = MAX(MIN(dl,dr)+shift,DBLE(MIN(sl,sr)))
1165 0 : x2 = MIN(MAX(dl,dr)+shift,DBLE(MAX(sl,sr)))
1166 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1167 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1168 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1169 :
1170 0 : shift = REAL(INT((sl/modulo) + REAL(1.5, DP)), DP) * modulo
1171 0 : x1 = MAX(MIN(dl,dr),DBLE(MIN(sl,sr))+shift)
1172 0 : x2 = MIN(MAX(dl,dr),DBLE(MAX(sl,sr))+shift)
1173 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1174 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1175 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1176 : END IF
1177 :
1178 0 : END SUBROUTINE OVL_RD
1179 : ! ------------------------------------------------------------------
1180 :
1181 : ! ------------------------------------------------------------------
1182 0 : SUBROUTINE OVL_DR (sl,sr,dl,dr,fs,fd,modulo)
1183 :
1184 : IMPLICIT NONE
1185 :
1186 : ! I/O
1187 : REAL (DP), INTENT(IN) :: sl, sr ! 'box' edges
1188 : REAL (SP), INTENT(IN) :: dl, dr ! 'box' edges
1189 : REAL (DP), INTENT(IN), OPTIONAL :: modulo ! shift for 'modulo' axis
1190 : REAL (DP), INTENT(OUT) :: fs, fd ! overlap fractions
1191 :
1192 : ! LOCAL
1193 : REAL (DP) :: x1, x2
1194 : REAL (DP) :: shift
1195 : INTEGER :: fx
1196 :
1197 0 : x1 = MAX(DBLE(MIN(dl,dr)),MIN(sl,sr))
1198 0 : x2 = MIN(DBLE(MAX(dl,dr)),MAX(sl,sr))
1199 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1200 0 : fs = (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1201 0 : fd = (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1202 :
1203 0 : IF (PRESENT(modulo)) THEN
1204 :
1205 0 : shift = REAL(INT((dl/modulo) + REAL(1.5, DP)), DP) * modulo
1206 0 : x1 = MAX(DBLE(MIN(dl,dr))+shift,MIN(sl,sr))
1207 0 : x2 = MIN(DBLE(MAX(dl,dr))+shift,MAX(sl,sr))
1208 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1209 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1210 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1211 :
1212 0 : shift = REAL(INT((sl/modulo) + REAL(1.5, DP)), DP) * modulo
1213 0 : x1 = MAX(DBLE(MIN(dl,dr)),MIN(sl,sr)+shift)
1214 0 : x2 = MIN(DBLE(MAX(dl,dr)),MAX(sl,sr)+shift)
1215 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1216 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1217 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1218 : END IF
1219 :
1220 0 : END SUBROUTINE OVL_DR
1221 : ! ------------------------------------------------------------------
1222 :
1223 : ! ------------------------------------------------------------------
1224 0 : SUBROUTINE OVL_DD (sl,sr,dl,dr,fs,fd,modulo)
1225 :
1226 : IMPLICIT NONE
1227 :
1228 : ! I/O
1229 : REAL (DP), INTENT(IN) :: sl, sr ! 'box' edges
1230 : REAL (DP), INTENT(IN) :: dl, dr ! 'box' edges
1231 : REAL (DP), INTENT(IN), OPTIONAL :: modulo ! shift for 'modulo' axis
1232 : REAL (DP), INTENT(OUT) :: fs, fd ! overlap fractions
1233 :
1234 : ! LOCAL
1235 : REAL (DP) :: x1, x2
1236 : REAL (DP) :: shift
1237 : INTEGER :: fx
1238 :
1239 0 : x1 = MAX(MIN(dl,dr),MIN(sl,sr))
1240 0 : x2 = MIN(MAX(dl,dr),MAX(sl,sr))
1241 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1242 0 : fs = (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1243 0 : fd = (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1244 :
1245 0 : IF (PRESENT(modulo)) THEN
1246 :
1247 0 : shift = REAL(INT((dl/modulo) + REAL(1.5, DP)), DP) * modulo
1248 0 : x1 = MAX(MIN(dl,dr)+shift,MIN(sl,sr))
1249 0 : x2 = MIN(MAX(dl,dr)+shift,MAX(sl,sr))
1250 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1251 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1252 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1253 :
1254 0 : shift = REAL(INT((sl/modulo) + REAL(1.5, DP)), DP) * modulo
1255 0 : x1 = MAX(MIN(dl,dr),MIN(sl,sr)+shift)
1256 0 : x2 = MIN(MAX(dl,dr),MAX(sl,sr)+shift)
1257 0 : fx = INT(SIGN(REAL(-0.5, DP),x2-x1)+REAL(1., DP))
1258 0 : fs = fs + (x2-x1)*REAL(fx, DP)/ABS(sr-sl)
1259 0 : fd = fd + (x2-x1)*REAL(fx, DP)/ABS(dr-dl)
1260 : END IF
1261 :
1262 0 : END SUBROUTINE OVL_DD
1263 : ! ------------------------------------------------------------------
1264 :
1265 : ! ------------------------------------------------------------------
1266 0 : SUBROUTINE OVL_1D_RR(s, d, mfs, mfd, lmod)
1267 :
1268 : IMPLICIT NONE
1269 :
1270 : ! I/O
1271 : REAL (SP), DIMENSION(:), INTENT(IN) :: s ! source bounds
1272 : REAL (SP), DIMENSION(:), INTENT(IN) :: d ! dest. bounds
1273 : ! REAL (DP), DIMENSION(:,:), POINTER :: mfs, mfd ! overlap matrices
1274 : REAL (DP), DIMENSION(:,:) :: mfs, mfd ! overlap matrices
1275 : LOGICAL, INTENT(IN) :: lmod ! modulo axis ?
1276 :
1277 : ! LOCAL
1278 : CHARACTER(LEN=*), PARAMETER :: substr = 'OVL_1D_RR'
1279 : INTEGER :: n,m ! dimensions
1280 : INTEGER :: i,j ! counter
1281 : ! INTEGER :: status
1282 : REAL (DP) :: vmod ! modulo value
1283 :
1284 : ! INIT
1285 0 : n = SIZE(s)-1
1286 0 : m = SIZE(d)-1
1287 :
1288 : ! ALLOCATE MEMORY
1289 : ! ALLOCATE(mfs(n,m),STAT=status)
1290 : ! CALL ERRMSG('OVL_1D_RR',status,1)
1291 : ! ALLOCATE(mfd(m,n),STAT=status)
1292 : ! CALL ERRMSG('OVL_1D_RR',status,2)
1293 :
1294 : ! CHECK CONSISTENCY
1295 0 : IF ((SIZE(mfs,1) /= n).OR.(SIZE(mfs,2) /= m)) THEN
1296 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1297 : END IF
1298 0 : IF ((SIZE(mfd,2) /= n).OR.(SIZE(mfd,1) /= m)) THEN
1299 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1300 : END IF
1301 :
1302 0 : IF (lmod) THEN
1303 0 : vmod = ABS(s(n+1))+ABS(s(1))
1304 0 : DO i=1, n
1305 0 : DO j=1, m
1306 0 : CALL OVL_RR(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i),vmod)
1307 : END DO
1308 : END DO
1309 : ELSE
1310 0 : DO i=1, n
1311 0 : DO j=1, m
1312 0 : CALL OVL_RR(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i))
1313 : END DO
1314 : END DO
1315 : END IF
1316 :
1317 0 : END SUBROUTINE OVL_1D_RR
1318 : ! ------------------------------------------------------------------
1319 :
1320 : ! ------------------------------------------------------------------
1321 0 : SUBROUTINE OVL_1D_RD(s, d, mfs, mfd, lmod)
1322 :
1323 : IMPLICIT NONE
1324 :
1325 : ! I/O
1326 : REAL (SP), DIMENSION(:), INTENT(IN) :: s ! source bounds
1327 : REAL (DP), DIMENSION(:), INTENT(IN) :: d ! dest. bounds
1328 : ! REAL (DP), DIMENSION(:,:), POINTER :: mfs, mfd ! overlap matrices
1329 : REAL (DP), DIMENSION(:,:) :: mfs, mfd ! overlap matrices
1330 : LOGICAL, INTENT(IN) :: lmod ! modulo axis ?
1331 :
1332 : ! LOCAL
1333 : CHARACTER(LEN=*), PARAMETER :: substr = 'OVL_1D_RD'
1334 : INTEGER :: n,m ! dimensions
1335 : INTEGER :: i,j ! counter
1336 : ! INTEGER :: status
1337 : REAL (DP) :: vmod ! modulo value
1338 :
1339 : ! INIT
1340 0 : n = SIZE(s)-1
1341 0 : m = SIZE(d)-1
1342 :
1343 : ! ALLOCATE MEMORY
1344 : ! ALLOCATE(mfs(n,m),STAT=status)
1345 : ! CALL ERRMSG('OVL_1D_RD',status,1)
1346 : ! ALLOCATE(mfd(m,n),STAT=status)
1347 : ! CALL ERRMSG('OVL_1D_RD',status,2)
1348 :
1349 : ! CHECK CONSISTENCY
1350 0 : IF ((SIZE(mfs,1) /= n).OR.(SIZE(mfs,2) /= m)) THEN
1351 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1352 : END IF
1353 0 : IF ((SIZE(mfd,2) /= n).OR.(SIZE(mfd,1) /= m)) THEN
1354 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1355 : END IF
1356 :
1357 0 : IF (lmod) THEN
1358 0 : vmod = ABS(s(n+1))+ABS(s(1))
1359 0 : DO i=1, n
1360 0 : DO j=1, m
1361 0 : CALL OVL_RD(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i),vmod)
1362 : END DO
1363 : END DO
1364 : ELSE
1365 0 : DO i=1, n
1366 0 : DO j=1, m
1367 0 : CALL OVL_RD(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i))
1368 : END DO
1369 : END DO
1370 : END IF
1371 :
1372 0 : END SUBROUTINE OVL_1D_RD
1373 : ! ------------------------------------------------------------------
1374 :
1375 : ! ------------------------------------------------------------------
1376 0 : SUBROUTINE OVL_1D_DR(s, d, mfs, mfd, lmod)
1377 :
1378 : IMPLICIT NONE
1379 :
1380 : ! I/O
1381 : REAL (DP), DIMENSION(:), INTENT(IN) :: s ! source bounds
1382 : REAL (SP), DIMENSION(:), INTENT(IN) :: d ! dest. bounds
1383 : ! REAL (DP), DIMENSION(:,:), POINTER :: mfs, mfd ! overlap matrices
1384 : REAL (DP), DIMENSION(:,:) :: mfs, mfd ! overlap matrices
1385 : LOGICAL, INTENT(IN) :: lmod ! modulo axis ?
1386 :
1387 : ! LOCAL
1388 : CHARACTER(LEN=*), PARAMETER :: substr = 'OVL_1D_DR'
1389 : INTEGER :: n,m ! dimensions
1390 : INTEGER :: i,j ! counter
1391 : ! INTEGER :: status
1392 : REAL (DP) :: vmod ! modulo value
1393 :
1394 : ! INIT
1395 0 : n = SIZE(s)-1
1396 0 : m = SIZE(d)-1
1397 :
1398 : ! ALLOCATE MEMORY
1399 : ! ALLOCATE(mfs(n,m),STAT=status)
1400 : ! CALL ERRMSG('OVL_1D_DR',status,1)
1401 : ! ALLOCATE(mfd(m,n),STAT=status)
1402 : ! CALL ERRMSG('OVL_1D_DR',status,2)
1403 :
1404 : ! CHECK CONSISTENCY
1405 0 : IF ((SIZE(mfs,1) /= n).OR.(SIZE(mfs,2) /= m)) THEN
1406 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1407 : END IF
1408 0 : IF ((SIZE(mfd,2) /= n).OR.(SIZE(mfd,1) /= m)) THEN
1409 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1410 : END IF
1411 :
1412 0 : IF (lmod) THEN
1413 0 : vmod = ABS(s(n+1))+ABS(s(1))
1414 0 : DO i=1, n
1415 0 : DO j=1, m
1416 0 : CALL OVL_DR(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i),vmod)
1417 : END DO
1418 : END DO
1419 : ELSE
1420 0 : DO i=1, n
1421 0 : DO j=1, m
1422 0 : CALL OVL_DR(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i))
1423 : END DO
1424 : END DO
1425 : END IF
1426 :
1427 0 : END SUBROUTINE OVL_1D_DR
1428 : ! ------------------------------------------------------------------
1429 :
1430 : ! ------------------------------------------------------------------
1431 0 : SUBROUTINE OVL_1D_DD(s, d, mfs, mfd, lmod)
1432 :
1433 : IMPLICIT NONE
1434 :
1435 : ! I/O
1436 : REAL (DP), DIMENSION(:), INTENT(IN) :: s ! source bounds
1437 : REAL (DP), DIMENSION(:), INTENT(IN) :: d ! dest. bounds
1438 : ! REAL (DP), DIMENSION(:,:), POINTER :: mfs, mfd ! overlap matrices
1439 : REAL (DP), DIMENSION(:,:) :: mfs, mfd ! overlap matrices
1440 : LOGICAL, INTENT(IN) :: lmod ! modulo axis ?
1441 :
1442 : ! LOCAL
1443 : CHARACTER(LEN=*), PARAMETER :: substr = 'OVL_1D_DD'
1444 : INTEGER :: n,m ! dimensions
1445 : INTEGER :: i,j ! counter
1446 : ! INTEGER :: status
1447 : REAL (DP) :: vmod ! modulo value
1448 :
1449 : ! INIT
1450 0 : n = SIZE(s)-1
1451 0 : m = SIZE(d)-1
1452 :
1453 : ! ALLOCATE MEMORY
1454 : ! ALLOCATE(mfs(n,m),STAT=status)
1455 : ! CALL ERRMSG('OVL_1D_DD',status,1)
1456 : ! ALLOCATE(mfd(m,n),STAT=status)
1457 : ! CALL ERRMSG('OVL_1D_DD',status,2)
1458 :
1459 : ! CHECK CONSISTENCY
1460 0 : IF ((SIZE(mfs,1) /= n).OR.(SIZE(mfs,2) /= m)) THEN
1461 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1462 : END IF
1463 0 : IF ((SIZE(mfd,2) /= n).OR.(SIZE(mfd,1) /= m)) THEN
1464 0 : CALL RGMSG(substr, RGMLE, 'INVALID SIZE OF MATRIX !')
1465 : END IF
1466 :
1467 0 : IF (lmod) THEN
1468 0 : vmod = ABS(s(n+1))+ABS(s(1))
1469 0 : DO i=1, n
1470 0 : DO j=1, m
1471 0 : CALL OVL_DD(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i),vmod)
1472 : END DO
1473 : END DO
1474 : ELSE
1475 0 : DO i=1, n
1476 0 : DO j=1, m
1477 0 : CALL OVL_DD(s(i),s(i+1),d(j),d(j+1),mfs(i,j),mfd(j,i))
1478 : END DO
1479 : END DO
1480 : END IF
1481 :
1482 0 : END SUBROUTINE OVL_1D_DD
1483 : ! ------------------------------------------------------------------
1484 :
1485 : ! ------------------------------------------------------------------
1486 0 : INTEGER FUNCTION POSITION(dim, vec)
1487 :
1488 : ! THIS FUNCTION CALCULATES THE POSITION NUMBER IN A
1489 : ! 1-D (LINEAR) ARRAY, GIVEN THAT THE ARRAY SHOULD BE INTERPRETED
1490 : ! AS N-D ARRAY WITH DIMENSIONS
1491 : ! dim = (d1, d2, d3, ..., dN)
1492 : ! OF THE ELEMENT
1493 : ! vec = (v1, v2, v3, ..., vN)
1494 : !
1495 :
1496 : IMPLICIT NONE
1497 :
1498 : ! I/O
1499 : INTEGER, DIMENSION(:), INTENT(IN) :: dim
1500 : INTEGER, DIMENSION(:), INTENT(IN) :: vec
1501 :
1502 : ! LOCAL
1503 : INTEGER :: i
1504 : INTEGER :: n
1505 : INTEGER :: dacc
1506 :
1507 0 : IF (SIZE(dim) /= SIZE(vec)) THEN
1508 0 : POSITION = 0
1509 0 : RETURN
1510 : END IF
1511 :
1512 0 : DO i=1, SIZE(dim)
1513 0 : IF (vec(i) > dim(i)) THEN
1514 0 : POSITION = 0
1515 0 : RETURN
1516 : END IF
1517 : END DO
1518 :
1519 0 : n = vec(1)
1520 0 : dacc = 1
1521 0 : DO i=2,SIZE(dim)
1522 0 : dacc = dacc*dim(i-1)
1523 0 : n = n + dacc*(vec(i)-1)
1524 : END DO
1525 :
1526 0 : POSITION = n
1527 :
1528 0 : END FUNCTION POSITION
1529 : ! ------------------------------------------------------------------
1530 :
1531 : ! ------------------------------------------------------------------
1532 0 : SUBROUTINE ELEMENT(dim, n, vec)
1533 :
1534 : ! THIS SUBROUTINE CALCULATES THE ELEMENT VECTOR
1535 : ! vec = (v1, v2, v3, ..., vN)
1536 : ! OF THE ELEMENT WITH POSITION n IN A
1537 : ! 1-D (LINEAR) ARRAY, GIVEN THAT THE ARRAY SHOULD BE INTERPRETED
1538 : ! AS N-D ARRAY WITH DIMENSIONS
1539 : ! dim = (d1, d2, d3, ..., dN)
1540 : !
1541 :
1542 : IMPLICIT NONE
1543 :
1544 : ! I/O
1545 : INTEGER, DIMENSION(:), INTENT(IN) :: dim ! dimension vector
1546 : INTEGER, INTENT(IN) :: n ! element in linear array
1547 : INTEGER, DIMENSION(:), POINTER :: vec ! element vector
1548 :
1549 : ! LOCAL
1550 : CHARACTER(LEN=*), PARAMETER :: substr = 'ELEMENT'
1551 : INTEGER :: m ! COPY of n
1552 : INTEGER :: i ! counter
1553 0 : INTEGER , DIMENSION(:), ALLOCATABLE :: dacc
1554 : INTEGER :: l ! length of dim
1555 : INTEGER :: status
1556 :
1557 0 : m = n
1558 0 : l = SIZE(dim)
1559 0 : IF (ASSOCIATED(vec)) THEN
1560 0 : IF (SIZE(vec) > 0) THEN
1561 0 : DEALLOCATE(vec, STAT=status)
1562 0 : CALL ERRMSG(substr,status,1)
1563 : END IF
1564 0 : NULLIFY(vec)
1565 : END IF
1566 0 : ALLOCATE(vec(l), STAT=status)
1567 0 : CALL ERRMSG(substr,status,2)
1568 0 : vec(:) = 0
1569 :
1570 0 : ALLOCATE(dacc(l))
1571 0 : dacc(1) = 1
1572 0 : DO i=2, l
1573 0 : dacc(i) = dacc(i-1)*dim(i-1)
1574 : END DO
1575 :
1576 0 : IF (m > dacc(l)*dim(l)) RETURN
1577 :
1578 0 : DO i=l, 2, -1
1579 0 : vec(i) = (m-1)/dacc(i)+1
1580 0 : m = m - (vec(i)-1)*dacc(i)
1581 : END DO
1582 0 : vec(1) = m
1583 :
1584 0 : DEALLOCATE(dacc, stat=STATUS)
1585 0 : CALL ERRMSG(substr,status,3)
1586 :
1587 0 : END SUBROUTINE ELEMENT
1588 : ! ------------------------------------------------------------------
1589 :
1590 : ! ------------------------------------------------------------------
1591 0 : RECURSIVE SUBROUTINE NREGRID(s, sax, dax, d, RG_TYPE &
1592 : ,sovl, dovl, rcnt &
1593 0 : ,gm, gnr, gsvec, gdvec, gdio &
1594 : ,gsf, gdf &
1595 : )
1596 :
1597 : IMPLICIT NONE
1598 :
1599 : ! I/O
1600 : TYPE (narray), DIMENSION(:), INTENT(IN) :: s ! source data fields
1601 : TYPE (narray), DIMENSION(:), POINTER :: d ! dest. data fields
1602 : TYPE (axis) , DIMENSION(:), INTENT(IN) :: sax ! axes (source)
1603 : TYPE (axis) , DIMENSION(:), INTENT(INOUT) :: dax ! axes (dest)
1604 : INTEGER , DIMENSION(:), INTENT(IN) :: RG_TYPE ! regridding type
1605 : REAL (DP), DIMENSION(:), POINTER :: sovl ! glob. overlap fraction
1606 : REAL (DP), DIMENSION(:), POINTER :: dovl ! glob. overlap fraction
1607 : INTEGER, DIMENSION(:), POINTER :: rcnt ! recursion level counter
1608 :
1609 : ! FOR RECURSION
1610 : INTEGER, OPTIONAL, INTENT(IN) :: gm ! current dimension
1611 : INTEGER, OPTIONAL, INTENT(IN) :: gnr ! no of dims to regrid
1612 : INTEGER, OPTIONAL, DIMENSION(:) :: gsvec ! cnt. vector (source)
1613 : INTEGER, OPTIONAL, DIMENSION(:) :: gdvec ! cnt vector (dest.)
1614 : INTEGER, OPTIONAL, DIMENSION(:) :: gdio ! dimension order
1615 : REAL (DP), OPTIONAL, INTENT(IN) :: gsf ! current overlap fraction
1616 : REAL (DP), OPTIONAL, INTENT(IN) :: gdf ! current overlap fraction
1617 :
1618 : ! FOR 1st STEP -> RECURSION
1619 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: dio ! dimension order
1620 : INTEGER :: m ! current dimenion
1621 : INTEGER :: nr ! no of dims to regrid
1622 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: svec ! cnt. vector (source)
1623 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: dvec ! cnt. vector (dest.)
1624 : REAL (DP) :: sf ! current overlap fraction
1625 : REAL (DP) :: df ! current overlap fraction
1626 : REAL (DP) :: sf0 ! current overlap fraction
1627 : REAL (DP) :: df0 ! current overlap fraction
1628 :
1629 : ! FOR INVARIANT DEPENDENT AXIS
1630 0 : TYPE (axis) , DIMENSION(:), ALLOCATABLE :: psax ! axes (source)
1631 0 : TYPE (axis) , DIMENSION(:), ALLOCATABLE :: pdax ! axes (dest)
1632 0 : TYPE (narray), DIMENSION(:), POINTER :: pd ! dest. data
1633 0 : REAL (DP), DIMENSION(:), POINTER :: psovl ! glob. overlap fraction
1634 0 : REAL (DP), DIMENSION(:), POINTER :: pdovl ! glob. overlap fraction
1635 0 : INTEGER, DIMENSION(:), POINTER :: prcnt ! recursion level counter
1636 :
1637 : ! LOCAL
1638 : CHARACTER(LEN=*), PARAMETER :: substr = 'NREGRID'
1639 : INTEGER :: i,j, k ! counter
1640 : INTEGER :: nvar ! number of fields
1641 : INTEGER :: id ! dimenions loop ounter
1642 : INTEGER :: nid ! count no. of dimensions
1643 : INTEGER :: ndim ! number of dimensions
1644 0 : REAL (DP) , DIMENSION(:), ALLOCATABLE :: sb ! source bounds
1645 0 : REAL (DP) , DIMENSION(:), ALLOCATABLE :: db ! dest. bounds
1646 : INTEGER :: vo ! offset for bound vec.
1647 : INTEGER :: vl ! length of bound vec.
1648 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: vec ! element vector
1649 0 : REAL (DP), DIMENSION(:,:), ALLOCATABLE :: ms ! overlap matrix
1650 0 : REAL (DP), DIMENSION(:,:), ALLOCATABLE :: md ! overlap matrix
1651 : INTEGER :: ns, nd ! dim. of ovl-matrices
1652 : INTEGER :: status ! error status
1653 : INTEGER :: vtype, svtype ! variable type
1654 : LOGICAL :: lflag ! switch
1655 : LOGICAL :: lsdef ! flag for defined dim.
1656 : LOGICAL :: lddef ! flag for defined dim.
1657 : LOGICAL :: lsdep ! flag for dependent dim.
1658 : INTEGER :: novl ! number of overlaps
1659 : REAL (DP) :: zso, zdo ! local overlap
1660 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: hio ! for output
1661 : CHARACTER(LEN=1000) :: hstr ! for output
1662 :
1663 : ! INIT
1664 : ! NUMBER OF DIMENSIONS
1665 0 : ndim = SIZE(sax)
1666 :
1667 0 : IF (PRESENT(gm)) THEN
1668 0 : m = gm
1669 : ELSE
1670 : m = 1
1671 : END IF
1672 :
1673 0 : IF (PRESENT(gnr)) THEN
1674 0 : nr = gnr
1675 : ELSE
1676 0 : nr = 0
1677 : END IF
1678 :
1679 0 : IF (PRESENT(gsvec)) THEN
1680 0 : ALLOCATE(svec(SIZE(gsvec)), STAT=status)
1681 0 : CALL ERRMSG(substr,status,1)
1682 0 : svec(:) = gsvec(:)
1683 : ELSE
1684 0 : ALLOCATE(svec(s(1)%n),STAT=status)
1685 0 : CALL ERRMSG(substr,status,2)
1686 0 : svec(:) = 1
1687 : END IF
1688 :
1689 0 : IF (PRESENT(gdvec)) THEN
1690 0 : ALLOCATE(dvec(SIZE(gdvec)), STAT=status)
1691 0 : CALL ERRMSG(substr,status,3)
1692 0 : dvec(:) = gdvec(:)
1693 : ELSE
1694 0 : ALLOCATE(dvec(s(1)%n),STAT=status)
1695 0 : CALL ERRMSG(substr,status,4)
1696 0 : dvec(:) = 1
1697 : END IF
1698 :
1699 0 : IF (PRESENT(gdio)) THEN
1700 0 : ALLOCATE(dio(SIZE(gdio)), STAT=status)
1701 0 : CALL ERRMSG(substr,status,5)
1702 0 : dio(:) = gdio(:)
1703 : ELSE
1704 0 : ALLOCATE(dio(ndim),STAT=status)
1705 0 : CALL ERRMSG(substr,status,6)
1706 0 : dio(:) = 0
1707 : END IF
1708 :
1709 0 : IF (PRESENT(gsf)) THEN
1710 0 : sf0 = gsf
1711 : ELSE
1712 0 : sf0 = REAL(1.0, DP)
1713 : END IF
1714 :
1715 0 : IF (PRESENT(gdf)) THEN
1716 0 : df0 = gdf
1717 : ELSE
1718 0 : df0 = REAL(1.0, DP)
1719 : END IF
1720 :
1721 0 : nvar = SIZE(s)
1722 :
1723 : ! ....................................................................
1724 : ! (A) ONLY TO BE DONE ONCE (FIRST STEP)
1725 0 : IF (m == 1) THEN
1726 : ! 1. CHECK SOME BASIC REQUIREMENTS
1727 0 : IF (SIZE(RG_TYPE) /= nvar) THEN
1728 : CALL RGMSG(substr,RGMLE, &
1729 0 : 'NUMBER OF REGRIDDING TYPE MISMATCH IN SOURCE !')
1730 : END IF
1731 : !
1732 0 : DO k=1, nvar
1733 0 : IF (s(k)%n /= SIZE(sax)) THEN
1734 0 : CALL RGMSG(substr,RGMLE,'NUMBER OF DIMENSIONS MISMATCH', .false.)
1735 0 : CALL RGMSG(substr,RGMLEC,'FOR VARIABLE ',k,':', .false.)
1736 0 : CALL RGMSG(substr,RGMLEC,'VAR_DIMENSIONS: ',s(k)%n,' ', .false.)
1737 0 : CALL RGMSG(substr,RGMLEC,'AXES : ',SIZE(sax),' ')
1738 : END IF
1739 : END DO
1740 : !
1741 0 : DO k=1, nvar
1742 0 : DO i=1, s(k)%n
1743 0 : vtype = QTYPE_NARRAY(sax(i)%dat)
1744 0 : IF (vtype /= VTYPE_UNDEF) THEN
1745 0 : IF (s(k)%dim(i) /= sax(i)%dat%dim(1)-1) THEN
1746 : CALL RGMSG(substr,RGMLE,'SOURCE DIMENSION MISMATCH', &
1747 0 : .false.)
1748 0 : CALL RGMSG(substr,RGMLEC,'FOR DIMENSION ',i,' ', .false.)
1749 0 : CALL RGMSG(substr,RGMLEC,'OF VARIABLE ',k,' ', .false.)
1750 : CALL RGMSG(substr,RGMLEC,'BETWEEN DATA WITH LENGTH ', &
1751 0 : s(k)%dim(i),' ', .false.)
1752 : CALL RGMSG(substr,RGMLEC,'AND AXIS WITH LENGTH ', &
1753 0 : sax(i)%dat%dim(1),' ')
1754 : END IF
1755 : END IF
1756 : END DO
1757 : END DO
1758 : !
1759 0 : IF (SIZE(dax) /= SIZE(sax)) THEN
1760 : CALL RGMSG(substr,RGMLE,'NUMBER OF DIMENSIONS OF SOURCE: ', &
1761 0 : SIZE(sax),' ',.false.)
1762 : CALL RGMSG(substr,RGMLEC,'NUMBER OF DIMENSIONS OF DEST. : ', &
1763 0 : SIZE(dax),' ')
1764 : END IF
1765 : !
1766 :
1767 : ! 2. CHECK DIMENSIONS
1768 0 : nid = 0
1769 0 : CALL RGMSG(substr, RGMLI, 'SCANNING ',ndim, ' DIMENSIONS ...')
1770 : ! 2.1: SOURCE DIMESNION TYPE
1771 0 : DO id=1,ndim
1772 : ! 2.1.1.: SOURCE DIM MUST BE REAL, DOUBLE, OR UNDEFINED
1773 0 : svtype = QTYPE_NARRAY(sax(id)%dat)
1774 0 : SELECT CASE(svtype)
1775 : CASE(VTYPE_REAL)
1776 : CALL RGMSG(substr, RGMLIC, '... SOURCE DIMENSION',id, &
1777 0 : ' IS OF TYPE REAL' )
1778 0 : lsdef = .true.
1779 : CASE(VTYPE_DOUBLE)
1780 : CALL RGMSG(substr, RGMLIC, '... SOURCE DIMENSION',id, &
1781 0 : ' IS OF TYPE DOUBLE PRECISION' )
1782 0 : lsdef = .true.
1783 : CASE(VTYPE_UNDEF)
1784 : CALL RGMSG(substr, RGMLIC, '... SOURCE DIMENSION',id, &
1785 0 : ' IS UNDEFINED' )
1786 0 : lsdef = .false.
1787 : CASE DEFAULT
1788 : !CASE(VTYPE_INT)
1789 : !CASE(VTYPE_CHAR)
1790 : !CASE(VTYPE_BYTE)
1791 : CALL RGMSG(substr, RGMLE, 'REGRIDDING IS ONLY POSSIBLE FOR', &
1792 0 : .false.)
1793 : CALL RGMSG(substr, RGMLEC, &
1794 0 : 'DIMENSIONS OF TYPE REAL OR DOUBLE PRECISION !')
1795 : END SELECT
1796 : ! 2.1.2 CHECK DEPENDENCIES
1797 0 : IF (sax(id)%ndp > 1) THEN
1798 0 : lsdep = .true.
1799 0 : DO i=1, sax(id)%ndp
1800 0 : IF (QTYPE_NARRAY(sax(sax(id)%dep(i))%dat) == VTYPE_UNDEF) THEN
1801 : CALL RGMSG(substr, RGMLE, 'SOURCE DIMENSION ', id, &
1802 0 : ' IS DEPENDENT ON', .false.)
1803 0 : CALL RGMSG(substr, RGMLEC, 'DIMENSION ',sax(id)%dep(i), &
1804 0 : ' WHICH IS UNDEFINED !')
1805 : END IF
1806 : END DO
1807 : ELSE
1808 : lsdep = .false.
1809 : END IF
1810 : ! 2.2.1.: DEST. DIM MUST BE REAL, DOUBLE, OR UNDEFINED
1811 0 : vtype = QTYPE_NARRAY(dax(id)%dat)
1812 0 : SELECT CASE(vtype)
1813 : CASE(VTYPE_REAL)
1814 : CALL RGMSG(substr, RGMLIC, '... DEST. DIMENSION',id, &
1815 0 : ' IS OF TYPE REAL' )
1816 0 : lddef = .true.
1817 : CASE(VTYPE_DOUBLE)
1818 : CALL RGMSG(substr, RGMLIC, '... DEST. DIMENSION',id, &
1819 0 : ' IS OF TYPE DOUBLE PRECISION')
1820 0 : lddef = .true.
1821 : CASE(VTYPE_UNDEF)
1822 : CALL RGMSG(substr, RGMLIC, '... DEST. DIMENSION',id, &
1823 0 : ' IS UNDEFINED')
1824 0 : lddef = .false.
1825 : CASE DEFAULT
1826 : !CASE(VTYPE_INT)
1827 : !CASE(VTYPE_CHAR)
1828 : !CASE(VTYPE_BYTE)
1829 : CALL RGMSG(substr, RGMLE, 'REGRIDDING IS ONLY POSSIBLE FOR', &
1830 0 : .false.)
1831 : CALL RGMSG(substr, RGMLEC, &
1832 0 : 'DIMENSIONS OF TYPE REAL OR DOUBLE PRECISION !')
1833 : END SELECT
1834 : ! 2.2.2 CHECK DEPENDENCIES
1835 0 : IF (dax(id)%ndp > 1) THEN
1836 0 : DO i=1, dax(id)%ndp
1837 0 : IF (QTYPE_NARRAY(dax(dax(id)%dep(i))%dat) == VTYPE_UNDEF) THEN
1838 : CALL RGMSG(substr, RGMLE, 'DEST. DIMENSION ', id, &
1839 0 : ' IS DEPENDENT ON', .false.)
1840 0 : CALL RGMSG(substr, RGMLEC, 'DIMENSION ',dax(id)%dep(i), &
1841 0 : ' WHICH IS UNDEFINED !')
1842 : END IF
1843 : END DO
1844 : END IF
1845 : !
1846 : ! 2.3 CHECK CONSISTENCY BETWEEN SOURCE AND DEST.
1847 : ! 2.3.1 UNDEF dimensions cannot be regridded
1848 0 : IF ((.NOT.lsdef).AND.(lddef)) THEN
1849 : CALL RGMSG(substr, RGMLE, &
1850 0 : 'DIMENSION OF TYPE UNDEFINED', .false.)
1851 : CALL RGMSG(substr, RGMLEC, &
1852 0 : 'CANNOT BE REGRIDDED ! FOR INVARIANT', .false.)
1853 : CALL RGMSG(substr, RGMLEC, &
1854 0 : 'DIMENSION, DEST. DIMENSION MUST ALSO', .false.)
1855 0 : CALL RGMSG(substr, RGMLEC, 'BE UNDEFINED !')
1856 : END IF
1857 : ! 2.3.2 INVARIANT, DEPENDENT DIMENSIONS HAVE TO BE PRE-REGRIDDED
1858 0 : IF (lsdep.AND.(.NOT.lddef)) THEN
1859 0 : CALL RGMSG(substr, RGMLI, 'FOUND INVRIANT DEPENDENT DIMENSION')
1860 : CALL RGMSG(substr, RGMLIC, &
1861 0 : ' >>> START REGRIDDING OF DEPENDENT DIMENSION ...')
1862 : ! REGRIDDING AXES
1863 0 : dax(id)%lm = sax(id)%lm
1864 0 : dax(id)%ndp = sax(id)%ndp
1865 0 : ALLOCATE(dax(id)%dep(dax(id)%ndp),STAT=status)
1866 0 : CALL ERRMSG(substr,status,7)
1867 0 : dax(id)%dep(:) = sax(id)%dep(:)
1868 0 : ALLOCATE(psax(sax(id)%ndp),STAT=status)
1869 0 : CALL ERRMSG(substr,status,8)
1870 0 : ALLOCATE(pdax(dax(id)%ndp),STAT=status)
1871 0 : CALL ERRMSG(substr,status,9)
1872 0 : CALL INIT_AXIS(psax(1)) ! FIRST DIMENSION IS INVARIANT ...
1873 0 : CALL INIT_AXIS(pdax(1)) ! ...
1874 0 : DO i=2, sax(id)%ndp
1875 0 : CALL COPY_AXIS(psax(i), sax(sax(id)%dep(i)))
1876 0 : CALL COPY_AXIS(pdax(i), dax(dax(id)%dep(i)))
1877 : END DO
1878 : !
1879 0 : CALL NREGRID( (/sax(id)%dat/) ,psax, pdax, pd &
1880 0 : ,(/ RG_INT/) , psovl, pdovl, prcnt)
1881 : CALL NREGRID_STAT(psax, pdax, psovl, pdovl &
1882 0 : ,(/sax(id)%dat/), pd, prcnt)
1883 : !
1884 0 : CALL COPY_NARRAY(dax(id)%dat, pd(1))
1885 0 : DEALLOCATE(psax, pdax, psovl, pdovl, STAT=status)
1886 0 : CALL ERRMSG(substr,status,10)
1887 0 : NULLIFY(psovl, pdovl)
1888 0 : DEALLOCATE(pd, STAT=status)
1889 0 : CALL ERRMSG(substr,status,11)
1890 0 : NULLIFY(pd)
1891 0 : DEALLOCATE(prcnt, STAT=status)
1892 0 : CALL ERRMSG(substr,status,12)
1893 0 : NULLIFY(prcnt)
1894 : CALL RGMSG(substr, RGMLIC, &
1895 0 : ' <<< ... END REGRIDDING OF DEPENDENT DIMENSION !')
1896 : END IF
1897 :
1898 : END DO
1899 0 : CALL RGMSG(substr, RGMLIC, '... END SCANNING ',ndim, ' DIMENSIONS !')
1900 :
1901 : ! 3. CALCULATE REGRIDDING ORDER OF DIMENSIONS
1902 0 : CALL RGMSG(substr, RGMLI, 'PROCESSING ORDER OF ',ndim, ' DIMENSIONS ...')
1903 0 : DO id=1,s(1)%n
1904 : ! 1st: ALL INDEPENDENT + TO BE REGRIDDED
1905 0 : vtype = QTYPE_NARRAY(dax(id)%dat)
1906 0 : lflag = (sax(id)%ndp == 1).AND. & ! source dim. independent
1907 : ! destination dimension available
1908 : ((vtype == VTYPE_REAL).OR.(vtype == VTYPE_DOUBLE)).AND. &
1909 : ! destination dimenions is independent
1910 0 : (dax(id)%ndp == 1)
1911 0 : IF (lflag) THEN
1912 0 : nid = nid + 1
1913 0 : dio(nid) = id
1914 0 : CALL RGMSG(substr, RGMLIC, ' ... DIMENSION ',id, ' IS INDEPENDENT')
1915 : END IF
1916 : END DO
1917 0 : DO id=1,s(1)%n
1918 : ! 2nd: ALL DEPENDENT + TO BE REGRIDDED
1919 0 : vtype = QTYPE_NARRAY(dax(id)%dat)
1920 0 : lflag = ((sax(id)%ndp > 1).OR. & ! source dim. indepndent
1921 : (dax(id)%ndp > 1)).AND. &
1922 : ! destination dimension available
1923 0 : ((vtype == VTYPE_REAL).OR.(vtype == VTYPE_DOUBLE))
1924 0 : IF (lflag) THEN
1925 0 : nid = nid + 1
1926 0 : dio(nid) = id
1927 0 : CALL RGMSG(substr, RGMLIC, ' ... DIMENSION ',id, ' IS DEPENDENT')
1928 : END IF
1929 : END DO
1930 : !
1931 0 : nr = nid ! THIS IS THE NUMBER OF DIMs TO REGRID
1932 : !
1933 0 : DO id=1,s(1)%n
1934 : ! 3rd: REST IS INVARIANT
1935 0 : vtype = QTYPE_NARRAY(dax(id)%dat)
1936 0 : lflag = (vtype /= VTYPE_REAL).AND.(vtype /= VTYPE_DOUBLE)
1937 0 : IF (lflag) THEN
1938 0 : nid = nid + 1
1939 0 : dio(nid) = id
1940 0 : CALL RGMSG(substr, RGMLIC, ' ... DIMENSION ',id, ' IS INVARIANT')
1941 : END IF
1942 : END DO
1943 : CALL RGMSG(substr, RGMLIC, &
1944 0 : '... END PROCESSING ORDER OF ',ndim, ' DIMENSIONS !')
1945 0 : IF (nid < s(1)%n) THEN
1946 0 : CALL RGMSG(substr, RGMLE, 'UNRECOGNIZED DIMENSION !', .false.)
1947 0 : CALL RGMSG(substr, RGMLEC, 'CHECK TYPE OF BOUNDS', .false.)
1948 0 : CALL RGMSG(substr, RGMLEC, '(MUST BE REAL OR DOUBLE PRECISION) !')
1949 : END IF
1950 0 : IF (nid > s(1)%n) THEN
1951 0 : CALL RGMSG(substr, RGMLE, 'AMBIGIOUS DIMENSION !')
1952 : END IF
1953 :
1954 : ! ! testing only
1955 : ! write(*,*) 'number of dims to regrid: ', nr
1956 : ! write(*,*) 'nid: ', nid
1957 : ! write(*,*) 'dio: ', dio
1958 :
1959 : ! 4. ALLOCATE SPACE FOR REGRIDDED DATA
1960 :
1961 : ! ! testing only
1962 : ! write(*,*) 'nvar: ', nvar
1963 :
1964 0 : ALLOCATE(d(nvar))
1965 0 : DO k=1, nvar ! LOOP OVER VARIABLES
1966 0 : nid = 1
1967 0 : d(k)%n = s(k)%n ! COPY NUMBER OF DIMENSIONS
1968 0 : ALLOCATE(d(k)%dim(d(k)%n),STAT=status)
1969 0 : CALL ERRMSG(substr,status,13)
1970 :
1971 : ! ! testing only
1972 : ! write(*,*) 'k: ', k
1973 : ! write(*,*) 'd(k)%n: ', d(k)%n
1974 :
1975 0 : DO id=1, s(k)%n
1976 0 : vtype = QTYPE_NARRAY(dax(id)%dat)
1977 0 : IF ((vtype == VTYPE_REAL).OR.(vtype == VTYPE_DOUBLE)) THEN
1978 : ! new destination dimension
1979 : ! INTERFACES HAVE +1 ENTRY
1980 : ! 1st DIMENSION MUST BE THE INDEPENDENT !!!
1981 0 : d(k)%dim(id) = (dax(id)%dat%dim(1) -1)
1982 0 : nid = nid * d(k)%dim(id)
1983 : ! COPY MODULO INFO
1984 0 : dax(id)%lm = sax(id)%lm
1985 :
1986 : ! ! testing only
1987 : ! write(*,*) 'id: ', id
1988 : ! write(*,*) 'd(k)%dim(id): ', d(k)%dim(id)
1989 :
1990 : ELSE
1991 : ! keep source dimension
1992 0 : d(k)%dim(id) = s(k)%dim(id)
1993 0 : nid = nid * d(k)%dim(id)
1994 0 : CALL COPY_AXIS(dax(id), sax(id))
1995 : END IF ! dest. dim or source dim.
1996 : END DO
1997 :
1998 0 : vtype = QTYPE_NARRAY(s(k))
1999 0 : SELECT CASE(vtype)
2000 : CASE(VTYPE_REAL)
2001 0 : ALLOCATE(d(k)%vr(nid),STAT=status)
2002 0 : CALL ERRMSG(substr,status,14)
2003 0 : d(k)%vr(:) = REAL(0.0, SP)
2004 : CASE(VTYPE_DOUBLE)
2005 0 : ALLOCATE(d(k)%vd(nid),STAT=status)
2006 0 : CALL ERRMSG(substr,status,15)
2007 0 : d(k)%vd(:) = REAL(0.0, DP)
2008 :
2009 : ! ! testing only
2010 : ! write(*,*) 'allocate double type:'
2011 : ! write(*,*) 'nid: ', nid
2012 : ! write(*,*) 'size d(k)%vd: ', size(d(k)%vd)
2013 :
2014 : CASE DEFAULT
2015 : !CASE(VTYPE_INT)
2016 : !CASE(VTYPE_CHAR)
2017 : !CASE(VTYPE_BYTE)
2018 : !CASE(VTYPE_UNDEF)
2019 : CALL RGMSG(substr, RGMLE, 'REGRIDDING IS ONLY POSSIBLE FOR', &
2020 0 : .false.)
2021 : CALL RGMSG(substr, RGMLEC, &
2022 0 : 'DATA OF TYPE REAL OR DOUBLE PRECISION !')
2023 : END SELECT
2024 :
2025 : END DO ! LOOP OVER VARIABLES
2026 :
2027 : ! 5. SOME DIAGNOSTIC OUTPUT
2028 : CALL RGMSG(substr, RGMLVM, ' DIMENSIONS TO REGRID : ', &
2029 0 : dio(1:nr),' ')
2030 : !
2031 0 : ALLOCATE(hio(nr),STAT=status)
2032 0 : CALL ERRMSG(substr,status,16)
2033 0 : DO i=1, nr
2034 0 : hio(i) = s(1)%dim(dio(i))
2035 : END DO
2036 : CALL RGMSG(substr, RGMLVM, ' - LENGTH(S) IN SOURCE : ', &
2037 0 : hio,' ')
2038 : !
2039 0 : DO i=1, nr
2040 0 : hio(i) = d(1)%dim(dio(i))
2041 : END DO
2042 : CALL RGMSG(substr, RGMLVM, ' - LENGTH(S) IN DEST. : ', &
2043 0 : hio,' ')
2044 0 : DEALLOCATE(hio,STAT=status)
2045 0 : CALL ERRMSG(substr,status,17)
2046 : !
2047 : CALL RGMSG(substr, RGMLVM, ' INVARIANT DIMENSIONS : ', &
2048 0 : dio(nr+1:ndim), ' ')
2049 : CALL RGMSG(substr, RGMLVM, ' NUMBER OF FIELDS : ', &
2050 0 : nvar,' ')
2051 : !
2052 0 : CALL RGMSG(substr, RGMLVM, ' INVARIANT DIMENSION LENGTH(S) : ')
2053 0 : ALLOCATE(hio(ndim-nr),STAT=status)
2054 0 : CALL ERRMSG(substr,status,18)
2055 0 : DO k=1,nvar
2056 0 : WRITE(hstr, *) ' FIELD ',k,': '
2057 0 : DO i=nr+1,ndim
2058 0 : hio(i-nr) = s(k)%dim(dio(i))
2059 : END DO
2060 0 : CALL RGMSG(substr, RGMLVM, TRIM(hstr), hio,' ')
2061 : END DO
2062 0 : DEALLOCATE(hio,STAT=status)
2063 0 : CALL ERRMSG(substr,status,18)
2064 :
2065 : ! 6. ALLOCATE SPACE FOR DIAGNOSTIC OUTPUT AND INITIALIZE
2066 0 : ALLOCATE(sovl(s(1)%n),STAT=status)
2067 0 : CALL ERRMSG(substr,status,19)
2068 0 : sovl(:) = REAL(0.0, DP)
2069 0 : ALLOCATE(dovl(d(1)%n),STAT=status)
2070 0 : CALL ERRMSG(substr,status,20)
2071 0 : dovl(:) = REAL(0.0, DP)
2072 0 : ALLOCATE(rcnt(ndim),STAT=status)
2073 0 : CALL ERRMSG(substr,status,21)
2074 0 : rcnt(:) = 0
2075 :
2076 0 : CALL RGMSG(substr, RGMLVL, ' STARTING NREGRID ... ')
2077 :
2078 : END IF ! (A) ONLY TO BE DONE ONCE AT FIRST STEP
2079 : ! ....................................................................
2080 :
2081 : ! SET CURRENT DIMENSION
2082 0 : nid = dio(m)
2083 0 : rcnt(m) = rcnt(m) + 1
2084 :
2085 : ! ! testing only
2086 : ! write(*,*) 'm: ', m
2087 : ! write(*,*) 'nid: ', nid
2088 : ! write(*,*) 'rcnt(m): ', rcnt(m)
2089 : ! write(*,*) 'dio: ', dio
2090 :
2091 : ! ....................................................................
2092 : ! (B) CONDITION FOR END OF RECURSION
2093 0 : IF (m == SIZE(dio)) THEN ! LAST (INVARIANT !) DIMENSION REACHED
2094 0 : zso = REAL(0.0, DP)
2095 0 : zdo = REAL(0.0, DP)
2096 0 : DO k=1, nvar ! LOOP OVER VARIABLES
2097 0 : SELECT CASE(RG_TYPE(k))
2098 : CASE(RG_EXT)
2099 0 : DO j=1, s(k)%dim(nid)
2100 0 : dvec(nid) = j
2101 0 : svec(nid) = j ! INVARIANT !!!
2102 : ! CALCULATE DEST. DATA POINT
2103 0 : vo = POSITION(d(k)%dim,dvec) ! POSITION IN DEST.
2104 0 : vl = POSITION(s(k)%dim,svec) ! POSITION IN SOURCE
2105 0 : vtype = QTYPE_NARRAY(s(k))
2106 : SELECT CASE(vtype)
2107 : CASE(VTYPE_REAL)
2108 0 : d(k)%vr(vo) = d(k)%vr(vo) + s(k)%vr(vl) * REAL(sf0, SP)
2109 : CASE(VTYPE_DOUBLE)
2110 0 : d(k)%vd(vo) = d(k)%vd(vo) + s(k)%vd(vl) * sf0
2111 : CASE DEFAULT
2112 : !CASE(VTYPE_INT)
2113 : !CASE(VTYPE_CHAR)
2114 : !CASE(VTYPE_BYTE)
2115 : !CASE(VTYPE_UNDEF)
2116 : CALL RGMSG(substr, RGMLE, &
2117 0 : 'REGRIDDING IS ONLY POSSIBLE FOR DATA', .false.)
2118 : CALL RGMSG(substr, RGMLEC, &
2119 0 : 'OF TYPE REAL OR DOUBLE PRECISION !')
2120 : END SELECT
2121 : ! OVERLAP IS 1 FOR INVARIANT DIMENSIONS
2122 : ! INTERVAL LENGTH, NUMBER OF VARIABLES
2123 : zso = zso + REAL(1.0, DP)/ &
2124 0 : (REAL(s(k)%dim(nid), DP)*REAL(nvar, DP))
2125 : zdo = zdo + REAL(1.0, DP)/ &
2126 0 : (REAL(s(k)%dim(nid), DP)*REAL(nvar, DP))
2127 : END DO
2128 : ! RESET FOR NEXT/PREVIOUS RECURSION STEP
2129 0 : dvec(nid) = 1
2130 0 : svec(nid) = 1
2131 : CASE(RG_INT, RG_IDX, RG_IXF)
2132 : !CASE(RG_IDX, RG_IXF)
2133 : ! NOTE: IDX VAR WAS CONVERTET TO INDEX-FRACTION
2134 0 : DO j=1, s(k)%dim(nid)
2135 0 : dvec(nid) = j
2136 0 : svec(nid) = j ! INVARIANT !!!
2137 : ! CALCULATE DEST. DATA POINT
2138 0 : vo = POSITION(d(k)%dim,dvec) ! POSITION IN DEST.
2139 0 : vl = POSITION(s(k)%dim,svec) ! POSITION IN SOURCE
2140 0 : vtype = QTYPE_NARRAY(s(k))
2141 : SELECT CASE(vtype)
2142 : CASE(VTYPE_REAL)
2143 0 : d(k)%vr(vo) = d(k)%vr(vo) + s(k)%vr(vl) * REAL(df0, SP)
2144 : CASE(VTYPE_DOUBLE)
2145 0 : d(k)%vd(vo) = d(k)%vd(vo) + s(k)%vd(vl) * df0
2146 : CASE DEFAULT
2147 : !CASE(VTYPE_INT)
2148 : !CASE(VTYPE_CHAR)
2149 : !CASE(VTYPE_BYTE)
2150 : !CASE(VTYPE_UNDEF)
2151 : CALL RGMSG(substr, RGMLE, &
2152 0 : 'REGRIDDING IS ONLY POSSIBLE FOR DATA', .false.)
2153 : CALL RGMSG(substr, RGMLEC, &
2154 0 : 'OF TYPE REAL OR DOUBLE PRECISION !')
2155 : END SELECT
2156 : ! OVERLAP IS 1 FOR INVARIANT DIMENSIONS:
2157 : ! INTERVAL LENGTH, NUMBER OF VARIABLES
2158 : zso = zso + REAL(1.0, DP)/ &
2159 0 : (REAL(s(k)%dim(nid), DP)*REAL(nvar, DP))
2160 : zdo = zdo + REAL(1.0, DP)/ &
2161 0 : (REAL(s(k)%dim(nid), DP)*REAL(nvar, DP))
2162 : END DO
2163 : ! RESET FOR NEXT/PREVIOUS RECURSION STEP
2164 0 : dvec(nid) = 1
2165 0 : svec(nid) = 1
2166 : !CASE DEFAULT
2167 : END SELECT
2168 :
2169 : END DO ! LOOP OVER VARIABLES
2170 :
2171 : ! SUM OVERLAP
2172 0 : sovl(nid) = sovl(nid) + zso
2173 0 : dovl(nid) = dovl(nid) + zdo
2174 :
2175 : ! CLEAN UP
2176 0 : DEALLOCATE(dio, STAT=status)
2177 0 : CALL ERRMSG(substr,status,22)
2178 0 : DEALLOCATE(svec, dvec, STAT=status)
2179 0 : CALL ERRMSG(substr,status,23)
2180 :
2181 0 : RETURN ! DONE ... ... END RECURSION
2182 :
2183 : END IF ! (B) CONDITION FOR END OF RECURSION
2184 : ! ....................................................................
2185 :
2186 : ! ....................................................................
2187 : ! (C) RECURSION STEP
2188 0 : IF (m <= nr) THEN ! (C1): DIMENSION TO REGRID
2189 :
2190 : ! ! testing only
2191 : ! write(*,*) 'regrid along dimension ', m
2192 :
2193 : ! GET SOURCE BOUNDS
2194 : ! 1st dimension is independent -> length of bounds vector
2195 0 : vl = sax(nid)%dat%dim(1) ! length is length of 1st (indep.) dim.
2196 0 : ALLOCATE(sb(vl),STAT=status)
2197 0 : CALL ERRMSG(substr,status,24)
2198 :
2199 0 : ALLOCATE(vec(sax(nid)%ndp),STAT=status) ! dependency element vector
2200 0 : CALL ERRMSG(substr,status,25)
2201 0 : vec(1) = 1 ! start at 1st position of independent dim.
2202 0 : DO i=2,sax(nid)%ndp ! loop over dependent dims
2203 0 : vec(i) = svec(sax(nid)%dep(i)) ! ... actual counter of dep. dims
2204 : END DO
2205 : ! CALCULATE OFFSET ...
2206 0 : vo = POSITION(sax(nid)%dat%dim, vec)
2207 : ! ... AND GET BOUNDS
2208 0 : vtype = QTYPE_NARRAY(sax(nid)%dat)
2209 : SELECT CASE(vtype)
2210 : CASE(VTYPE_REAL)
2211 0 : sb(:) = sax(nid)%dat%vr(vo:(vo+vl-1))
2212 : CASE(VTYPE_DOUBLE)
2213 0 : sb(:) = sax(nid)%dat%vd(vo:(vo+vl-1))
2214 : CASE DEFAULT
2215 : !CASE(VTYPE_INT)
2216 : !CASE(VTYPE_CHAR)
2217 : !CASE(VTYPE_BYTE)
2218 : !CASE(VTYPE_UNDEF)
2219 0 : CALL RGMSG(substr, RGMLE, 'SOURCE BOUNDS MUST BE OF TYPE', .false.)
2220 0 : CALL RGMSG(substr, RGMLEC, 'REAL OR DOUBLE PRECISION !')
2221 : END SELECT
2222 : !
2223 0 : DEALLOCATE(vec,STAT=status)
2224 0 : CALL ERRMSG(substr,status,26)
2225 :
2226 : ! ! testing only
2227 : ! write(*,*) 'source bounds: vl, vo, sb ', vl, vo, sb
2228 :
2229 : ! GET DEST. BOUNDS
2230 : ! 1st dimension is independent -> length of bounds vector
2231 0 : vl = dax(nid)%dat%dim(1) ! length is length of 1st (indep.) dim.
2232 0 : ALLOCATE(db(vl),STAT=status)
2233 0 : CALL ERRMSG(substr,status,27)
2234 :
2235 0 : ALLOCATE(vec(dax(nid)%ndp),STAT=status) ! dependency element vector
2236 0 : CALL ERRMSG(substr,status,28)
2237 0 : vec(1) = 1 ! start at 1st position of independent dim.
2238 0 : DO i=2,dax(nid)%ndp ! loop over dependent dims
2239 0 : vec(i) = dvec(dax(nid)%dep(i)) ! ... actual counter of dep. dims
2240 : END DO
2241 : ! CALCULATE OFFSET ...
2242 0 : vo = POSITION(dax(nid)%dat%dim, vec)
2243 : ! ... AND GET BOUNDS
2244 0 : vtype = QTYPE_NARRAY(dax(nid)%dat)
2245 : SELECT CASE(vtype)
2246 : CASE(VTYPE_REAL)
2247 0 : db(:) = dax(nid)%dat%vr(vo:(vo+vl-1))
2248 : CASE(VTYPE_DOUBLE)
2249 0 : db(:) = dax(nid)%dat%vd(vo:(vo+vl-1))
2250 : CASE DEFAULT
2251 : !CASE(VTYPE_INT)
2252 : !CASE(VTYPE_CHAR)
2253 : !CASE(VTYPE_BYTE)
2254 : !CASE(VTYPE_UNDEF)
2255 0 : CALL RGMSG(substr, RGMLE, 'DEST. BOUNDS MUST BE OF TYPE', .false.)
2256 0 : CALL RGMSG(substr, RGMLEC, 'REAL OR DOUBLE PRECISION !')
2257 : END SELECT
2258 : !
2259 0 : DEALLOCATE(vec,STAT=status)
2260 0 : CALL ERRMSG(substr,status,29)
2261 :
2262 : ! ! testing only
2263 : ! write(*,*) 'destination bounds: vl, vo, db ', vl, vo, db
2264 :
2265 : ! REGRID ALONG THIS DIMENSION
2266 : ! ALLOCATE SPACE FOR MATRICES
2267 0 : ns = SIZE(sb)-1
2268 0 : nd = SIZE(db)-1
2269 0 : ALLOCATE(ms(ns,nd), STAT=status)
2270 0 : CALL ERRMSG(substr,status,30)
2271 0 : ALLOCATE(md(nd,ns), STAT=status)
2272 0 : CALL ERRMSG(substr,status,31)
2273 : CALL OVL_1D(sb, db, ms, md, sax(nid)%lm)
2274 :
2275 : ! SAVE OVERLAP OF NEXT DIMENSION CALCULATED SO FAR ...
2276 0 : zso = sovl(dio(m+1))
2277 0 : zdo = dovl(dio(m+1))
2278 : ! ... AND RESET VECTOR ELEMENT TO ZERO
2279 0 : sovl(dio(m+1)) = REAL(0.0, DP)
2280 0 : dovl(dio(m+1)) = REAL(0.0, DP)
2281 0 : novl = 0
2282 : ! LOOP OVER ALL MATRIX ELEMENTS AND REGRID ALONG NEXT DIMENSION ...
2283 : ! ... IF MATRIX ELEMENT IS NON-ZERO
2284 : !
2285 0 : DO i=1, ns
2286 0 : svec(nid) = i
2287 0 : DO j=1, nd
2288 0 : dvec(nid) = j
2289 0 : sf = sf0*ms(svec(nid),dvec(nid))
2290 0 : df = df0*md(dvec(nid),svec(nid))
2291 0 : IF (sf > 0.0) THEN ! NON-ZERO
2292 0 : novl = novl + 1
2293 : ! CALCULATE TOTAL OVERLAP FRACTION FOR THIS DIMENSION
2294 0 : sovl(nid) = sovl(nid) + &
2295 0 : ms(svec(nid),dvec(nid)) &
2296 0 : *(sb(svec(nid)+1)-sb(svec(nid))) &
2297 0 : /(sb(ns+1)-sb(1))
2298 0 : dovl(nid) = dovl(nid) + &
2299 0 : md(dvec(nid),svec(nid)) &
2300 0 : *(db(dvec(nid)+1)-db(dvec(nid))) &
2301 0 : /(db(nd+1)-db(1))
2302 : ! REGRID NEXT DIMENSION
2303 : CALL NREGRID(s, sax, dax, d, RG_TYPE, sovl, dovl, rcnt &
2304 : ,m+1, nr, svec, dvec &
2305 0 : ,dio, sf, df)
2306 : END IF
2307 : END DO
2308 : END DO
2309 :
2310 : ! AVERAGE OVERLAP FRACTIONS:
2311 : ! RESTORE OLD OVERLAP FRACTION AND ADD NEW
2312 0 : IF (novl /= 0) THEN
2313 0 : sovl(dio(m+1)) = zso + sovl(dio(m+1))/REAL(novl, DP)
2314 0 : dovl(dio(m+1)) = zdo + dovl(dio(m+1))/REAL(novl, DP)
2315 : ELSE
2316 0 : sovl(dio(m+1)) = zso
2317 0 : dovl(dio(m+1)) = zdo
2318 : END IF
2319 :
2320 : ! RESET FOR NEXT/PREVIOUS RECURSION STEP
2321 0 : svec(nid) = 1
2322 0 : dvec(nid) = 1
2323 :
2324 : ! CLEAN
2325 0 : DEALLOCATE(ms, md, stat=status)
2326 0 : CALL ERRMSG(substr,status,32)
2327 0 : DEALLOCATE(sb, db, stat=status)
2328 0 : CALL ERRMSG(substr,status,33)
2329 :
2330 : ELSE ! (C2): INVARIANT DIMENSION: ( m > nr ) .......................
2331 :
2332 : ! SAVE OVERLAP OF NEXT DIMENSION CALCULATED SO FAR ...
2333 0 : zso = sovl(dio(m+1))
2334 0 : zdo = dovl(dio(m+1))
2335 : ! ... AND RESET VECTOR ELEMENT TO ZERO
2336 0 : sovl(dio(m+1)) = REAL(0.0, DP)
2337 0 : dovl(dio(m+1)) = REAL(0.0, DP)
2338 : ! LOOP OVER ALL DIAGONAL MATRIX ELEMENTS ...
2339 : ! ... AND REGRID NEXT DIMENSION
2340 0 : DO j=1, s(1)%dim(nid)
2341 0 : dvec(nid) = j
2342 0 : svec(nid) = j ! INVARIANT !!!
2343 : ! REGRID NEXT DIMENSION
2344 : CALL NREGRID(s, sax, dax, d, RG_TYPE, sovl, dovl, rcnt &
2345 : ,m+1, nr, svec, dvec &
2346 0 : ,dio, sf0, df0 )
2347 : ! OVERLAP IS 1 FOR INVARIANT DIMENSIONS
2348 0 : sovl(nid) = sovl(nid)+REAL(1.0, DP)/REAL(s(1)%dim(nid), DP)
2349 0 : dovl(nid) = dovl(nid)+REAL(1.0, DP)/REAL(s(1)%dim(nid), DP)
2350 : END DO
2351 :
2352 : ! AVERAGE OVERLAP FRACTIONS:
2353 : ! RESTORE OLD OVERLAP FRACTION AND ADD NEW
2354 0 : sovl(dio(m+1)) = zso + sovl(dio(m+1))/REAL(s(1)%dim(nid), DP)
2355 0 : dovl(dio(m+1)) = zdo + dovl(dio(m+1))/REAL(s(1)%dim(nid), DP)
2356 :
2357 : ! RESET FOR NEXT/PREVIOUS RECURSION STEP
2358 0 : svec(nid) = 1
2359 0 : dvec(nid) = 1
2360 :
2361 : END IF ! (C) RECURSION STEP
2362 : ! ....................................................................
2363 :
2364 : ! CLEAN UP
2365 0 : DEALLOCATE(dio, STAT=status)
2366 0 : CALL ERRMSG(substr,status,34)
2367 0 : DEALLOCATE(svec, dvec, STAT=status)
2368 0 : CALL ERRMSG(substr,status,35)
2369 :
2370 0 : END SUBROUTINE NREGRID
2371 : ! ------------------------------------------------------------------
2372 :
2373 : ! ------------------------------------------------------------------
2374 0 : SUBROUTINE NREGRID_STAT(sax, dax, sovl, dovl, nai, nao, rcnt)
2375 :
2376 : IMPLICIT NONE
2377 :
2378 : ! I/O
2379 : TYPE (axis), DIMENSION(:), INTENT(IN) :: sax, dax
2380 : REAL (DP), DIMENSION(:), INTENT(IN) :: sovl, dovl
2381 : TYPE (narray), DIMENSION(:), INTENT(IN) :: nai, nao
2382 : INTEGER, DIMENSION(:), INTENT(IN) :: rcnt
2383 :
2384 : ! LOCAL
2385 : INTEGER :: i
2386 : INTEGER :: vtype
2387 : REAL (DP) :: div
2388 :
2389 0 : WRITE(*,*) ' NREGRID STATISTICS:'
2390 0 : WRITE(*,*) ' .......................................................'
2391 0 : WRITE(*,*) ' NO. OF RECURSION LEVELS : ',SIZE(rcnt)
2392 0 : WRITE(*,*) ' RECURSION LEVELS PROCESSED: ',rcnt
2393 0 : WRITE(*,*) ' '
2394 0 : DO i=1, SIZE(sovl)
2395 : !
2396 0 : IF (i > 1) THEN
2397 0 : IF (rcnt(i-1) /= 0) THEN
2398 0 : div = REAL(rcnt(i-1), DP)
2399 : ELSE
2400 : div = REAL(1.0, DP)
2401 : END IF
2402 : ELSE
2403 : div = REAL(1.0, DP)
2404 : END IF
2405 : !
2406 0 : IF (sax(i)%lm.OR.dax(i)%lm) THEN
2407 0 : WRITE(*,*) ' DIMENSION ',i,': (MODULO)'
2408 : ELSE
2409 0 : WRITE(*,*) ' DIMENSION ',i,':'
2410 : END IF
2411 0 : vtype = QTYPE_NARRAY(sax(i)%dat)
2412 0 : SELECT CASE(vtype)
2413 : CASE(VTYPE_REAL)
2414 0 : WRITE(*,*) ' SOURCE: ',MINVAL(sax(i)%dat%vr),&
2415 0 : MAXVAL(sax(i)%dat%vr)
2416 : CASE(VTYPE_DOUBLE)
2417 0 : WRITE(*,*) ' SOURCE: ',MINVAL(sax(i)%dat%vd),&
2418 0 : MAXVAL(sax(i)%dat%vd)
2419 : CASE(VTYPE_INT)
2420 0 : WRITE(*,*) ' SOURCE: ',MINVAL(sax(i)%dat%vi),&
2421 0 : MAXVAL(sax(i)%dat%vi)
2422 : CASE(VTYPE_BYTE)
2423 0 : WRITE(*,*) ' SOURCE: ',MINVAL(sax(i)%dat%vb),&
2424 0 : MAXVAL(sax(i)%dat%vb)
2425 : CASE(VTYPE_CHAR)
2426 0 : WRITE(*,*) ' SOURCE: CHAR NOT SUPPORTED'
2427 : CASE DEFAULT
2428 0 : WRITE(*,*) ' SOURCE: <UNDEFINED>'
2429 : END SELECT
2430 0 : WRITE(*,*) ' SOURCE REGION COVERED ON AVERAGE: ',sovl(i)/div
2431 : !
2432 0 : vtype = QTYPE_NARRAY(dax(i)%dat)
2433 0 : SELECT CASE(vtype)
2434 : CASE(VTYPE_REAL)
2435 0 : WRITE(*,*) ' DEST. : ',MINVAL(dax(i)%dat%vr),&
2436 0 : MAXVAL(dax(i)%dat%vr)
2437 : CASE(VTYPE_DOUBLE)
2438 0 : WRITE(*,*) ' DEST. : ',MINVAL(dax(i)%dat%vd),&
2439 0 : MAXVAL(dax(i)%dat%vd)
2440 : CASE(VTYPE_INT)
2441 0 : WRITE(*,*) ' DEST. : ',MINVAL(dax(i)%dat%vi),&
2442 0 : MAXVAL(dax(i)%dat%vi)
2443 : CASE(VTYPE_BYTE)
2444 0 : WRITE(*,*) ' DEST. : ',MINVAL(dax(i)%dat%vb),&
2445 0 : MAXVAL(dax(i)%dat%vb)
2446 : CASE(VTYPE_CHAR)
2447 0 : WRITE(*,*) ' DEST. : CHAR NOT SUPPORTED'
2448 : CASE DEFAULT
2449 0 : WRITE(*,*) ' DEST. : <UNDEFINED>'
2450 : END SELECT
2451 : !
2452 0 : WRITE(*,*) ' DEST. REGION COVERED ON AVERAGE: ',dovl(i)/div
2453 : END DO
2454 : !
2455 0 : WRITE(*,*) ' '
2456 0 : WRITE(*,*) ' VARIABLE RANGES: '
2457 0 : DO i=1, SIZE(nai)
2458 0 : vtype = QTYPE_NARRAY(nai(i))
2459 0 : SELECT CASE(vtype)
2460 : CASE(VTYPE_REAL)
2461 0 : WRITE(*,*) ' (',i,'): <REAL>'
2462 0 : WRITE(*,*) ' SOURCE: ',MINVAL(nai(i)%vr),MAXVAL(nai(i)%vr)
2463 0 : WRITE(*,*) ' DEST. : ',MINVAL(nao(i)%vr),MAXVAL(nao(i)%vr)
2464 : CASE(VTYPE_DOUBLE)
2465 0 : WRITE(*,*) ' (',i,'): <DOUBLE PRECISION>'
2466 0 : WRITE(*,*) ' SOURCE: ',MINVAL(nai(i)%vd),MAXVAL(nai(i)%vd)
2467 0 : WRITE(*,*) ' DEST. : ',MINVAL(nao(i)%vd),MAXVAL(nao(i)%vd)
2468 : CASE(VTYPE_INT)
2469 0 : WRITE(*,*) ' (',i,'): <INTEGER>'
2470 0 : WRITE(*,*) ' SOURCE: ',MINVAL(nai(i)%vi),MAXVAL(nai(i)%vi)
2471 0 : WRITE(*,*) ' DEST. : ',MINVAL(nao(i)%vi),MAXVAL(nao(i)%vi)
2472 : CASE(VTYPE_BYTE)
2473 0 : WRITE(*,*) ' (',i,'): <BYTE>'
2474 0 : WRITE(*,*) ' SOURCE: ',MINVAL(nai(i)%vb),MAXVAL(nai(i)%vb)
2475 0 : WRITE(*,*) ' DEST. : ',MINVAL(nao(i)%vb),MAXVAL(nao(i)%vb)
2476 : CASE(VTYPE_CHAR)
2477 0 : WRITE(*,*) ' (',i,'): <CHAR>'
2478 0 : WRITE(*,*) ' SOURCE: CHAR NOT SUPPORTED'
2479 0 : WRITE(*,*) ' DEST. : CHAR NOT SUPPORTED'
2480 : CASE DEFAULT
2481 0 : WRITE(*,*) ' (',i,'): <UNDEFINED>'
2482 0 : WRITE(*,*) ' SOURCE: <UNDEFINED>'
2483 0 : WRITE(*,*) ' DEST. : <UNDEFINED>'
2484 : END SELECT
2485 : END DO
2486 : !
2487 0 : WRITE(*,*) ' .......................................................'
2488 :
2489 0 : END SUBROUTINE NREGRID_STAT
2490 : ! ------------------------------------------------------------------
2491 :
2492 : ! ------------------------------------------------------------------
2493 0 : SUBROUTINE ERRMSG(routine, status, pos)
2494 :
2495 : IMPLICIT NONE
2496 :
2497 : ! I/O
2498 : CHARACTER(LEN=*), INTENT(IN) :: routine
2499 : INTEGER, INTENT(IN) :: status
2500 : INTEGER, INTENT(IN) :: pos
2501 :
2502 0 : IF (status == 0) THEN
2503 : RETURN
2504 : ELSE
2505 0 : CALL RGMSG(routine, RGMLE, 'ERROR STATUS ',status,' ', .false.)
2506 0 : CALL RGMSG(routine, RGMLEC, 'AT POSITION ',pos,' !')
2507 : END IF
2508 :
2509 : END SUBROUTINE ERRMSG
2510 : ! ------------------------------------------------------------------
2511 :
2512 : ! ------------------------------------------------------------------
2513 0 : SUBROUTINE RGMSG_C(routine, level, c, lstop)
2514 :
2515 : #if defined(MPI)
2516 : USE messy_ncregrid_mpi, ONLY: ncregrid_abort
2517 : #endif
2518 :
2519 : IMPLICIT NONE
2520 :
2521 : ! I/O
2522 : CHARACTER(LEN=*), INTENT(IN) :: routine
2523 : INTEGER, INTENT(IN) :: level
2524 : CHARACTER(LEN=*), INTENT(IN) :: c
2525 : LOGICAL, INTENT(IN), OPTIONAL :: lstop
2526 :
2527 : ! LOCAL
2528 : LOGICAL :: llstop
2529 :
2530 : ! INIT
2531 0 : IF (PRESENT(lstop)) THEN
2532 0 : llstop = lstop
2533 : ELSE
2534 0 : IF ((level == RGMLE).OR.(level == RGMLEC)) THEN
2535 : llstop = .true. ! STOP ON ERROR
2536 : ELSE
2537 0 : llstop = .false.
2538 : END IF
2539 : END IF
2540 :
2541 0 : SELECT CASE(level)
2542 : CASE(RGMLE) ! ERROR MESSAGE
2543 0 : IF (IAND(MSGMODE, MSGMODE_E) == MSGMODE_E) THEN
2544 0 : WRITE(*,*) '*** ',TRIM(routine),' ERROR: '
2545 0 : WRITE(*,*) ' ',TRIM(c)
2546 : END IF
2547 : CASE(RGMLEC) ! ERROR MESSAGE CONTINUED
2548 0 : IF (IAND(MSGMODE, MSGMODE_E) == MSGMODE_E) THEN
2549 0 : WRITE(*,*) ' ',TRIM(c)
2550 : END IF
2551 : CASE(RGMLVL) ! LITTLE VERBOSE
2552 0 : IF (IAND(MSGMODE, MSGMODE_VL) == MSGMODE_VL) THEN
2553 0 : WRITE(*,*) TRIM(c)
2554 : END IF
2555 : CASE(RGMLVLC) ! LITTLE VERBOSE CONTINUED
2556 0 : IF (IAND(MSGMODE, MSGMODE_VL) == MSGMODE_VL) THEN
2557 0 : WRITE(*,*) ' ',TRIM(c)
2558 : END IF
2559 : CASE(RGMLW) ! WARNING MESSAGE
2560 0 : IF (IAND(MSGMODE, MSGMODE_W) == MSGMODE_W) THEN
2561 0 : WRITE(*,*) '+++ ',TRIM(routine),' WARNING: '
2562 0 : WRITE(*,*) ' ',TRIM(c)
2563 : END IF
2564 : CASE(RGMLWC) ! WARNING MESSAGE CONTINUED
2565 0 : IF (IAND(MSGMODE, MSGMODE_W) == MSGMODE_W) THEN
2566 0 : WRITE(*,*) ' ',TRIM(c)
2567 : END IF
2568 : CASE(RGMLVM) ! MEDIUM VERBOSE
2569 0 : IF (IAND(MSGMODE, MSGMODE_VM) == MSGMODE_VM) THEN
2570 0 : WRITE(*,*) TRIM(c)
2571 : END IF
2572 : CASE(RGMLVMC) ! MEDIUM VERBOSE CONTINUED
2573 0 : IF (IAND(MSGMODE, MSGMODE_VM) == MSGMODE_VM) THEN
2574 0 : WRITE(*,*) ' ',TRIM(c)
2575 : END IF
2576 : CASE(RGMLI) ! INFO MESSAGE
2577 0 : IF (IAND(MSGMODE, MSGMODE_I) == MSGMODE_I) THEN
2578 0 : WRITE(*,*) '=== ',TRIM(routine),' INFO: '
2579 0 : WRITE(*,*) ' ',TRIM(c)
2580 : END IF
2581 : CASE(RGMLIC) ! INFO MESSAGE CONTINUED
2582 0 : IF (IAND(MSGMODE, MSGMODE_I) == MSGMODE_I) THEN
2583 0 : WRITE(*,*) ' ',TRIM(c)
2584 : END IF
2585 : CASE DEFAULT
2586 : END SELECT
2587 :
2588 : #if defined(MPI)
2589 : IF (llstop) CALL ncregrid_abort('ncregrid')
2590 : #else
2591 0 : IF (llstop) STOP
2592 : #endif
2593 :
2594 0 : END SUBROUTINE RGMSG_C
2595 : ! ------------------------------------------------------------------
2596 :
2597 : ! ------------------------------------------------------------------
2598 0 : SUBROUTINE RGMSG_I(routine, level, c1, i, c2, lstop)
2599 :
2600 : IMPLICIT NONE
2601 :
2602 : ! I/O
2603 : CHARACTER(LEN=*), INTENT(IN) :: routine
2604 : INTEGER, INTENT(IN) :: level
2605 : CHARACTER(LEN=*), INTENT(IN) :: c1
2606 : INTEGER, INTENT(IN) :: i
2607 : CHARACTER(LEN=*), INTENT(IN) :: c2
2608 : LOGICAL, INTENT(IN), OPTIONAL :: lstop
2609 :
2610 : ! LOCAL
2611 : CHARACTER(LEN=1000) :: istr = ''
2612 :
2613 0 : WRITE(istr,*) TRIM(c1),i,TRIM(c2)
2614 0 : CALL RGMSG_C(routine, level, TRIM(istr), lstop)
2615 :
2616 0 : END SUBROUTINE RGMSG_I
2617 : ! ------------------------------------------------------------------
2618 :
2619 : ! ------------------------------------------------------------------
2620 0 : SUBROUTINE RGMSG_IA(routine, level, c1, i, c2, lstop)
2621 :
2622 : IMPLICIT NONE
2623 :
2624 : ! I/O
2625 : CHARACTER(LEN=*), INTENT(IN) :: routine
2626 : INTEGER, INTENT(IN) :: level
2627 : CHARACTER(LEN=*), INTENT(IN) :: c1
2628 : INTEGER, DIMENSION(:), INTENT(IN) :: i
2629 : CHARACTER(LEN=*), INTENT(IN) :: c2
2630 : LOGICAL, INTENT(IN), OPTIONAL :: lstop
2631 :
2632 : ! LOCAL
2633 : CHARACTER(LEN=1000) :: istr = ''
2634 :
2635 0 : WRITE(istr,*) TRIM(c1),i,TRIM(c2)
2636 0 : CALL RGMSG_C(routine, level, TRIM(istr), lstop)
2637 :
2638 0 : END SUBROUTINE RGMSG_IA
2639 : ! ------------------------------------------------------------------
2640 :
2641 : ! ------------------------------------------------------------------
2642 0 : SUBROUTINE RGMSG_R(routine, level, c1, r, c2, lstop)
2643 :
2644 : IMPLICIT NONE
2645 :
2646 : ! I/O
2647 : CHARACTER(LEN=*), INTENT(IN) :: routine
2648 : INTEGER, INTENT(IN) :: level
2649 : CHARACTER(LEN=*), INTENT(IN) :: c1
2650 : REAL, INTENT(IN) :: r
2651 : CHARACTER(LEN=*), INTENT(IN) :: c2
2652 : LOGICAL, INTENT(IN), OPTIONAL :: lstop
2653 :
2654 : ! LOCAL
2655 : CHARACTER(LEN=1000) :: rstr = ''
2656 :
2657 0 : WRITE(rstr,*) TRIM(c1),r,TRIM(c2)
2658 0 : CALL RGMSG_C(routine, level, TRIM(rstr), lstop)
2659 :
2660 0 : END SUBROUTINE RGMSG_R
2661 : ! ------------------------------------------------------------------
2662 :
2663 : ! ******************************************************************
2664 0 : END MODULE MESSY_NCREGRID_BASE
2665 : ! ******************************************************************
|
