LCOV - code coverage report
Current view: | top level - hemco/HEMCO/src/Core - messy_ncregrid_base.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 1113 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 33 | 0.0 %

              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             : ! ******************************************************************

Generated by: LCOV version 1.14