*DECK NXTMCD
      SUBROUTINE NXTMCD(IPGEO ,IPTRK ,IPRINT,NDIM  ,ILCELL,NBOCEL,
     >                  MAXMSH,MAXREG,MAXPIN,NBTCLS,ITSYM ,IDFEX ,
     >                  DCMESH,NAGGEO,ITURN ,IDIRR ,NFSUR ,NFREG ,
     >                  MXGSUR,MXGREG)
*
*----------
*
*Purpose:
* Create a multicell description for the geometry and verify
* if intrinsic cell symmetry required by boundary conditions
* are satisfied.
*
*Copyright:
* Copyright (C) 2005 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
*Author(s):
* G. Marleau.
*
*Parameters: input
* IPGEO   pointer to the GEOMETRY data structure.
* IPTRK   pointer to the TRACKING data structure.
* IPRINT  print level.
* NDIM    problem dimensions.
* ILCELL  cell level.
* NBOCEL  number of cells in original geometry.
* MAXMSH  maximum number of elements in MESH array.
* MAXREG  maximum number of elements in MIX array.
* MAXPIN  maximum number of pins in clusters.
* NBTCLS  total number of cluster geometry.
* ITSYM   array to identify the symmetry to test for each original
*         cell where:
*         ITSYM(1,*) identify $X$ symmetry;
*         ITSYM(2,*) identify $Y$ symmetry;
*         ITSYM(3,*) identify $Z$ symmetry;
*         ITSYM(4,*) identify $X-Y$ symmetry.
*         A value of 0 indicate that the geometry does not need
*         to be verified while a value of 1 implies a verification
*         of the geometry.
* IDFEX   identify faces associated with external boundary for a
*         generating cell and number of times this cell is used. Here:
*         IDFEX( 1,*)  identify bottom $U$ hexagonal face;
*         IDFEX( 2,*)  identify top $U$ hexagonal face;
*         IDFEX( 3,*)  identify bottom $V$ hexagonal face;
*         IDFEX( 4,*)  identify top $V$ hexagonal face;
*         IDFEX( 5,*)  identify bottom $Z$ face;
*         IDFEX( 6,*)  identify top $Z$ face;
*         IDFEX( 7,*)  not used;
*         IDFEX( 8,*)  not used;
*         IDFEX( 9,*)  identify bottom $W$ hexagonal face;
*         IDFEX(10,*) identify top $W$ hexagonal face.
* DCMESH  global mesh for each cell.
* NAGGEO  names of generating geometries.
* ITURN   geometry turns.
* IDIRR   direction of cell (1 for XYZ, 2 for YZX and 3 for ZXY).
*         Note: for CAR3D without pins IDIRR=1
*         for CAR3D with pins IDIRR specified by pins direction.
*
*Parameters: output
* NFSUR   final number of surfaces.
* NFREG   final number of regions.
* MXGSUR  maximum number of surfaces for any geometry.
* MXGREG  maximum number of region for any geometry.
*
*Reference:
*  G. Marleau,
*  New Geometries Processing in DRAGON: The NXT: Module,
*  Report IGE-260, Polytechnique Montreal,
*  Montreal, 2005.
*
*----------
*
      USE              GANLIB
      IMPLICIT         NONE
*----
*  Subroutine arguments
*----
      TYPE(C_PTR)      IPGEO,IPTRK
*      INTEGER          IPGEO,IPTRK
      INTEGER          IPRINT,NDIM,
     >                 ILCELL,NBOCEL,
     >                 MAXMSH,MAXREG,MAXPIN,NBTCLS,
     >                 ITSYM(4,NBOCEL),IDFEX(0:10,NBOCEL)
      DOUBLE PRECISION DCMESH(3,2,NBOCEL)
      INTEGER          NAGGEO(3,NBOCEL),ITURN(NBOCEL),
     >                 IDIRR(NBOCEL),NFSUR,NFREG,MXGSUR,MXGREG
*----
*  Local parameters
*----
      INTEGER          IOUT
      CHARACTER        NAMSBR*6
      PARAMETER       (IOUT=6,NAMSBR='NXTMCD')
      INTEGER          NSTATE
      PARAMETER       (NSTATE=40)
      DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO
      PARAMETER       (DCUTOF=1.0D-8,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
*----
*  Functions
*----
      DOUBLE PRECISION XDRCST,PI
*----
*  Local variables
*----
      INTEGER          IOFNAC,ICEL,ITC,ITYPG,ITRN,IDIRC,ITST,MAXMSS,
     >                 NBGCLS,ICLS,ITYPC,IGCLS,ILEV,
     >                 NMIX,NMIXS,NMIXC,NMIXCS,
     >                 NREG,NSUR,NREGS,NSURS,NREGC,NSURC,NREGCS,NSURCS,
     >                 NREGN,NSURN,NREGCN,NSURCN
      INTEGER          ISTATG(NSTATE),ISTATC(NSTATE),
     >                 IEDIMX(NSTATE),IEDIMP(NSTATE)
      INTEGER          NM(4),NMC(4),NMS(4),NMCS(4),
     >                 ITSYMC(4),IGSYM,IDSYM
      CHARACTER        NAMGG*12,NAMCL*12,NAMREC*12
      INTEGER          ILCMLN,ILCMTY,ILCMLX,ILCMLY
      INTEGER          NAREG,NASUR,NEREN
      INTEGER          NPIN,IP,IPINT,ICPIN,NPIR,IDGPP,NBPIN
      DOUBLE PRECISION DELTA,TWOPI,DROUT,DHPIN
*----
*  Allocatable arrays
*----
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ITPIN,ISPLT,MIX
      REAL, ALLOCATABLE, DIMENSION(:) :: RMESH,RPIN,APIN,CPINX,CPINY
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DAMESH,DRAPIN
*
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: MIXC,ISPLTR
      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: NAGCLS,INDXSR
      INTEGER, ALLOCATABLE, DIMENSION(:)     :: IDSUR,IDREG
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: DAMESR
*----
*  Data
*----
      CHARACTER        CLEV(2)*1
      SAVE             CLEV
      DATA             CLEV /'C','P'/
*----
*  Scratch storage allocation
*   ITPIN   temporary array for reading pin type.
*   ISPLT   temporary array for reading cell mesh split.
*   MIX     temporary array for reading cell mixture (including HMIX).
*   RMESH   temporary vector for reading cell mesh array.
*   RPIN    temporary vector for reading pin positions.
*   APIN    temporary vector for reading pin angles.
*   DAMESH  temporary vector for storing global mesh array.
*   DRAPIN  temporary vector for storing global pin positions.
*----
      ALLOCATE(ITPIN(3,MAXPIN),ISPLT(MAXMSH,4),MIX(MAXREG,2))
      ALLOCATE(RMESH(0:MAXMSH),RPIN(MAXPIN),APIN(MAXPIN),
     >         CPINX(MAXPIN),CPINY(MAXPIN))
      ALLOCATE(DAMESH(-1:MAXMSH,4),DRAPIN(-1:4,MAXPIN))
*----
*  Processing starts:
*  print routine openning output header if required
*  and initialize various parameters.
*----
      IF(IPRINT .GE. 10) THEN
        WRITE(IOUT,6000) NAMSBR
      ENDIF
      PI=XDRCST('Pi',' ')
      TWOPI=DTWO*PI
      ALLOCATE(MIXC(MAXREG,2,2),NAGCLS(3,NBTCLS),ISPLTR(MAXMSH,4,2))
      MXGSUR=0
      MXGREG=0
      IOFNAC=0
      IGCLS=0
      ICPIN=0
      NAREG=0
      NASUR=0
      DO ICEL=1,NBOCEL
        IF(ILCELL .EQ. 1) THEN
          WRITE(NAMGG,'(3A4)') (NAGGEO(ITC,ICEL),ITC=1,3)
          CALL LCMSIX(IPGEO,NAMGG,1)
        ELSE
          NAMGG='/           '
        ENDIF
        IF(IPRINT .GE. 10) THEN
          WRITE(IOUT,6010) ICEL,NAMGG
        ENDIF
        ISTATG(:NSTATE)=0
        CALL LCMGET(IPGEO,'STATE-VECTOR',ISTATG)
        ITYPG=ISTATG(1)
        IF(ITYPG .EQ. 20) THEN
          IF(ISTATG(2) .EQ. 0) ITYPG=5
        ELSE IF(ITYPG .EQ. 21 .OR. ITYPG .EQ. 22 .OR.
     >          ITYPG .EQ. 23) THEN
          IF(ISTATG(2) .EQ. 0) ITYPG=7
        ELSE IF(ITYPG .EQ. 26) THEN
          IF(ISTATG(2) .EQ. 0) ITYPG=12
        ELSE IF(ITYPG .EQ. 27) THEN
          IF(ISTATG(2) .EQ. 0) ITYPG=13
        ENDIF
        NM(4)=ISTATG(2)
        NM(1)=ISTATG(3)
        NM(2)=ISTATG(4)
        NM(3)=ISTATG(5)
        NMIX=ISTATG(6)
        ITRN=ITURN(ICEL)
        IDIRC=IDIRR(ICEL)
        NBGCLS=ISTATG(13)
        ILEV=1
        ITPIN(:3,:MAXPIN)=0
        ICPIN=0
        IGSYM=0
        DO IDSYM=1,4
          ITSYMC(IDSYM)=0
          IGSYM=IGSYM+ABS(ITSYM(IDSYM,ICEL))
        ENDDO
        IF(NBGCLS .NE. 0) THEN
*----
*  Get pin cluster geometry information
*----
          DRAPIN(-1:4,:MAXPIN)=DZERO
          ILEV=2
          CALL LCMGET(IPGEO,'CLUSTER',NAGCLS(1,IOFNAC+1))
          DO ICLS=1,NBGCLS
            IGCLS=IGCLS+1
            WRITE(NAMCL,'(3A4)') (NAGCLS(ITC,IOFNAC+ICLS),ITC=1,3)
            CALL LCMSIX(IPGEO,NAMCL,1)
            IF(IPRINT .GE. 10) THEN
              WRITE(IOUT,6011) ICLS,NAMCL
            ENDIF
            ISTATC(:NSTATE)=0
            CALL LCMGET(IPGEO,'STATE-VECTOR',ISTATC)
            ITYPC=ISTATC(1)
            IF(ITYPC .EQ. 20) THEN
              IF(ISTATC(2) .EQ. 0) ITYPC=5
            ELSE IF(ITYPC .EQ. 21 .OR. ITYPC .EQ. 22 .OR.
     >              ITYPC .EQ. 23) THEN
              IF(ISTATC(2) .EQ. 0) ITYPC=7
            ENDIF
*----
*  Cartesian or Annular pin
*----
            IPINT=0
            IF(ITYPC .EQ.  5 .OR. ITYPC .EQ.  7 .OR.
     >         ITYPC .EQ. 20 .OR. ITYPC .EQ. 23      ) THEN
              IPINT=-3
            ELSE IF(ITYPC .EQ.  3 .OR. ITYPC .EQ.  6 ) THEN
              IPINT=3
            ELSE IF(ITYPC .EQ. 21) THEN
              IPINT=-1
            ELSE IF(ITYPC .EQ. 10) THEN
              IPINT=1
            ELSE IF(ITYPC .EQ. 22) THEN
              IPINT=-2
            ELSE IF(ITYPC .EQ. 11 ) THEN
              IPINT=2
            ENDIF
            NMC(4)=ISTATC(2)
            NMC(1)=ISTATC(3)
            NMC(2)=ISTATC(4)
            NMC(3)=ISTATC(5)
            NMIXC=ISTATC(6)
            DAMESH(-1:MAXMSH,:4)=DZERO
            ISPLT(:MAXMSH,:4)=1
            MIX(:MAXREG,:2)=0
            MIXC(:MAXREG,:2,:2)=0
            CALL NXTEGI(IPGEO ,IPRINT,ITYPC ,MAXMSH,NMIXC ,NMC   ,
     >                  MAXMSS,NMCS  ,NREGC ,NREGCS,NSURC ,NSURCS,
     >                  MIX   ,ISPLT ,DAMESH,RMESH ,MIXC  )
            ALLOCATE(IDSUR(NSURCS),IDREG(NREGCS))
            IDSUR(:NSURCS)=0
            IDREG(:NREGCS)=0
            NEREN=NREGCS+NSURCS
            ALLOCATE(DAMESR(MAXMSS+2,4,2))
            DAMESR(:MAXMSS+2,:4,:2)=DZERO
            NMIXCS=NREGCS
            DROUT=DAMESH(NMC(4),4)
            IF(ABS(IPINT) .EQ. 3) THEN
*----
*  For pins rotate symmetry by pi/2 around z axis
*----
              IF(ITSYM(1,ICEL) .NE. 0) THEN
                ITSYMC(2)=ITSYM(1,ICEL)
              ENDIF
              IF(ITSYM(2,ICEL) .NE. 0) THEN
                ITSYMC(1)=ITSYM(2,ICEL)
              ENDIF
              ITSYMC(3)=ITSYM(3,ICEL)
              ITSYMC(4)=4*ITSYM(4,ICEL)
            ELSE IF(ABS(IPINT) .EQ. 2) THEN
*----
*  For pins rotate symmetry by pi/2  around y axis
*----
              IF(ITSYM(1,ICEL) .NE. 0) THEN
                ITSYMC(3)=ITSYM(1,ICEL)
              ENDIF
              IF(ITSYM(3,ICEL) .NE. 0) THEN
                ITSYMC(1)=ITSYM(3,ICEL)
              ENDIF
              ITSYMC(2)=ITSYM(2,ICEL)
              ITSYMC(4)=4*ITSYM(4,ICEL)
            ELSE IF(ABS(IPINT) .EQ. 1) THEN
*----
*  For pins rotate symmetry by pi/2  around x axis
*----
              IF(ITSYM(3,ICEL) .NE. 0) THEN
                ITSYMC(2)=ITSYM(3,ICEL)
              ENDIF
              IF(ITSYM(2,ICEL) .NE. 0) THEN
                ITSYMC(3)=ITSYM(2,ICEL)
              ENDIF
              ITSYMC(1)=ITSYM(1,ICEL)
              ITSYMC(4)=4*ITSYM(4,ICEL)
            ENDIF
*----
*  Test if pin cells satisfy symmetry options
*----
            IF(IGSYM .GT. 0) THEN
*----
*  Rotate pin cells and test for intrinsic symmetry
*----
              ITST=1
              CALL NXTRIS(IPRINT,ITYPC ,MAXMSH,NREGC ,ITRN  ,ITST  ,
     >                    ITSYMC(1)    ,NMC   ,MIX   ,ISPLT ,DAMESH,
     >                    NMCS  ,MIXC  ,ISPLTR,
     >                    DAMESR)
            ELSE
*----
*  Rotate pin cells
*----
              ITST=0
              CALL NXTRIS(IPRINT,ITYPC ,MAXMSH,NREGC ,ITRN  ,ITST  ,
     >                    ITSYMC(1)    ,NMC   ,MIX   ,ISPLT ,DAMESH,
     >                    NMCS  ,MIXC  ,ISPLTR,DAMESR)
            ENDIF
*----
*  Final mesh after unfolding
*----
            IF(MOD(ITRN,2) .EQ. 0) THEN
*----
*  First turn geometry type if required
*----
              IF(ITYPC .EQ. 10) THEN
                ITYPC=11
              ELSE IF(ITYPC .EQ.11) THEN
                ITYPC=10
              ELSE IF(ITYPC .EQ.21) THEN
                ITYPC=22
              ELSE IF(ITYPC .EQ.22) THEN
                ITYPC=21
              ENDIF
            ENDIF
*----
*  Then create new mesh
*----
            CALL NXTSGI(IPTRK ,IPRINT,MAXMSH,ITYPC ,IGCLS ,ILEV  ,
     >                  MAXMSS,NMIXC ,NMC   ,MIX   ,DAMESH,ISPLT ,
     >                  NMIXCS,NMCS  ,DAMESR,ITSYMC(1),
     >                  NREGCS,NSURCS,NREGCN,NSURCN,NEREN ,
     >                  IDREG ,IDSUR )
            ALLOCATE(INDXSR(5,NREGCS+NSURCS+1))
            INDXSR(:5,:NREGCS+NSURCS+1)=0
            NBPIN=0
            MXGSUR=MAX(MXGSUR,NSURCS)
            MXGREG=MAX(MXGREG,NREGCS)
            IEDIMP(:NSTATE)=0
            IEDIMP(1)=ITYPC
            IEDIMP(2)=NMCS(4)
            IEDIMP(3)=NMCS(1)
            IEDIMP(4)=NMCS(2)
            IEDIMP(5)=NMCS(3)
            IEDIMP(6)=NMIXCS
            IEDIMP(7)=ILEV
            IEDIMP(8)=NREGCS
            IEDIMP(9)=NSURCS
            IEDIMP(10)=NREGCN
            IEDIMP(11)=NSURCN
            IEDIMP(12)=NAREG+1
            CALL NXTVOL(IPTRK ,IPRINT,MAXMSS,ITYPC ,IDIRC ,IGCLS ,
     >                  ILEV  ,NMCS  ,NREGCS,NSURCS,NREGCN,NSURCN,
     >                  MAXPIN,NBPIN ,ITPIN ,DRAPIN, IDREG,IDSUR ,
     >                  DAMESR,INDXSR,NAREG )
            IEDIMP(13)=NAREG
            DEALLOCATE(DAMESR)
*----
*  Read pin locations
*----
            CALL LCMGET(IPGEO,'NPIN',NPIN)
            CALL LCMLEN(IPGEO,'RPIN',ILCMLN,ILCMTY)
            IF(ILCMLN .GE. 1) THEN
              IF(ILCMLN .EQ. 1) THEN
                CALL LCMGET(IPGEO,'RPIN',RPIN)
                DO IP=1,NPIN
                  DRAPIN(0,ICPIN+IP)=DBLE(RPIN(1))
                ENDDO
              ELSE IF(ILCMLN .EQ. NPIN) THEN
                CALL LCMGET(IPGEO,'RPIN',RPIN)
                DO IP=1,NPIN
                  DRAPIN(0,ICPIN+IP)=DBLE(RPIN(IP))
                ENDDO
              ELSE
                CALL XABORT(NAMSBR//
     >          ': Length of RPIN vector is invalid')
              ENDIF
              CALL LCMLEN(IPGEO,'APIN',ILCMLN,ILCMTY)
              IF(ILCMLN .EQ. 0) THEN
                DRAPIN(-1,ICPIN+1)=DZERO
                DELTA=(DTWO*PI)/DBLE(NPIN)
                DO IP=2,NPIN
                  DRAPIN(-1,ICPIN+IP)=DRAPIN(-1,ICPIN+IP-1)+DELTA
                ENDDO
              ELSE IF(ILCMLN .EQ. 1) THEN
                CALL LCMGET(IPGEO,'APIN',APIN)
                DRAPIN(-1,ICPIN+1)=DBLE(APIN(1))
                DELTA=(DTWO*PI)/DBLE(NPIN)
                DO IP=2,NPIN
                  DRAPIN(-1,ICPIN+IP)=DRAPIN(-1,ICPIN+IP-1)+DELTA
                ENDDO
              ELSE IF(ILCMLN .EQ. NPIN) THEN
                CALL LCMGET(IPGEO,'APIN',APIN)
                DO IP=1,NPIN
                  DRAPIN(-1,ICPIN+IP)=DBLE(APIN(IP))
                ENDDO
              ELSE
                CALL XABORT(NAMSBR//
     >          ': Length of APIN vector is invalid')
              ENDIF
            ELSE
              CALL LCMLEN(IPGEO,'CPINX',ILCMLX,ILCMTY)
              CALL LCMLEN(IPGEO,'CPINY',ILCMLY,ILCMTY)
              IF(ILCMLX .EQ. NPIN .AND. ILCMLY .EQ. NPIN) THEN
                CALL LCMGET(IPGEO,'CPINX',CPINX)
                CALL LCMGET(IPGEO,'CPINY',CPINY)
              ELSE
                CALL XABORT(NAMSBR//
     >': (RPIN,APIN) or (CPINX,CPINY) are absent for pin cluster')
              ENDIF
              DO IP=1,NPIN
                DRAPIN(0,ICPIN+IP)=DBLE(SQRT(CPINX(IP)**2+CPINY(IP)**2))
                DRAPIN(-1,ICPIN+IP)=DBLE(ATAN2(CPINY(IP),CPINX(IP)))
              ENDDO
            ENDIF
*----
*  Find pin height
*----
            IDGPP=MOD(IDIRC+1,3)+1
            IF(NDIM .EQ. 3) THEN
              DHPIN=DAMESH(NMC(IDGPP),IDGPP)-DAMESH(0,IDGPP)
            ELSE
              DHPIN=DONE
            ENDIF
*----
*  Replace angles in range 0 to 2*Pi
*----
            DO IP=1,NPIN
              IF(DRAPIN(-1,ICPIN+IP) .GE. DZERO) THEN
                NPIR=INT((DRAPIN(-1,ICPIN+IP)+DCUTOF)/TWOPI)
                DRAPIN(-1,ICPIN+IP)=DRAPIN(-1,ICPIN+IP)
     >                             -DBLE(NPIR)*TWOPI
              ELSE
                NPIR=INT((DRAPIN(-1,ICPIN+IP)-DCUTOF)/TWOPI)
                DRAPIN(-1,ICPIN+IP)=DRAPIN(-1,ICPIN+IP)
     >                             -DBLE(NPIR-1)*TWOPI
              ENDIF
              DRAPIN(IDGPP,ICPIN+IP)=DHPIN
              DRAPIN(4,ICPIN+IP)=DROUT
              ITPIN(1,ICPIN+IP)=ICEL
              ITPIN(2,ICPIN+IP)=IGCLS
              ITPIN(3,ICPIN+IP)=IPINT
            ENDDO
            IF(IGSYM .GT. 0) THEN
*----
*  Test pin in cluster for symmetry
*----
              CALL NXTTPS(IPRINT,NPIN  ,IDGPP ,ITSYM(1,ICEL),
     >                    DRAPIN(-1,ICPIN+1))
            ENDIF
*----
*  Renumber pin surfaces
*----
            IEDIMP(14)=NASUR+1
            CALL NXTRPS(IPTRK ,IPRINT,ITYPC ,IGCLS ,ILEV  ,
     >                  NREGCS,NSURCS,NSURCN,IDFEX(0,ICEL),
     >                  INDXSR,DHPIN ,DCMESH(1,1,ICEL),
     >                  NASUR ,IDSUR )
            DEALLOCATE(INDXSR)
            DEALLOCATE(IDSUR,IDREG)
            IEDIMP(15)=NASUR
            IEDIMP(16)=ICPIN+1
            IEDIMP(17)=NPIN
            IEDIMP(20)=IDIRC
            WRITE(NAMREC,'(A1,I8.8,A3)') CLEV(ILEV),IGCLS,'DIM'
            CALL LCMPUT(IPTRK,NAMREC,NSTATE,1,IEDIMP)
            ICPIN=ICPIN+NPIN
            CALL LCMSIX(IPGEO,NAMCL,2)
          ENDDO
*----
*  Save cluster information for this cell
*----
          ILEV=1
          IOFNAC=IOFNAC+NBGCLS
          WRITE(NAMREC,'(A1,I8.8,A3)') CLEV(ILEV),ICEL,'PIN'
          CALL LCMPUT(IPTRK,NAMREC,6*ICPIN,4,DRAPIN)
          WRITE(NAMREC,'(A1,I8.8,A3)') CLEV(ILEV),ICEL,'PNT'
          CALL LCMPUT(IPTRK,NAMREC,3*ICPIN,1,ITPIN)
        ENDIF
        DAMESH(-1:MAXMSH,:4)=DZERO
        ISPLT(:MAXMSH,:4)=1
        MIX(:MAXREG,:2)=0
        MIXC(:MAXREG,:2,:2)=0
        IF(IPRINT .GE. 10) THEN
          WRITE(IOUT,6012)
        ENDIF
*----
*  Get main geometry information
*----
        CALL NXTEGI(IPGEO ,IPRINT,ITYPG ,MAXMSH,NMIX  ,NM    ,
     >              MAXMSS,NMS   ,NREG  ,NREGS ,NSUR  ,NSURS ,
     >              MIX   ,ISPLT ,DAMESH,
     >              RMESH ,MIXC  )
        NMIXS=NREGS
        ALLOCATE(IDREG(NREGS),IDSUR(NSURS))
        IDREG(:NREGS)=0
        IDSUR(:NSURS)=0
        NEREN=NREGS+NSURS
        ALLOCATE(DAMESR(MAXMSS+2,4,2))
        DAMESR(:MAXMSS+2,:4,:2)=DZERO
        IF(ITYPG .EQ. 12 .OR. ITYPG .EQ. 13 .OR. 
     >     ITYPG .EQ. 26 .OR. ITYPG .EQ. 27) THEN
          IF(IGSYM .GT. 0) THEN
*----
*  Rotate global geometry and test for intrinsic symmetry
*----
            ITST=1
            CALL NXTRTS(IPRINT,ITYPG ,MAXMSH,NREG  ,ITRN  ,ITST,
     >                  ITSYM(1,ICEL),NM    ,MIX   ,ISPLT ,DAMESH,
     >                  NMS   ,MIXC  ,ISPLTR,DAMESR)
          ELSE
*----
*  Rotate global geometry
*----
            ITST=0
            CALL NXTRTS(IPRINT,ITYPG ,MAXMSH,NREG  ,ITRN  ,ITST,
     >                  ITSYM(1,ICEL),NM    ,MIX   ,ISPLT ,DAMESH,
     >                  NMS   ,MIXC  ,ISPLTR,DAMESR)
          ENDIF
        ELSE
          IF(IGSYM .GT. 0) THEN
*----
*  Rotate global geometry and test for intrinsic symmetry
*----
            ITST=1
            CALL NXTRIS(IPRINT,ITYPG ,MAXMSH,NREG  ,ITRN  ,ITST,
     >                  ITSYM(1,ICEL),NM    ,MIX   ,ISPLT ,DAMESH,
     >                  NMS   ,MIXC  ,ISPLTR,DAMESR)
          ELSE
*----
*  Rotate global geometry
*----
            ITST=0
            CALL NXTRIS(IPRINT,ITYPG ,MAXMSH,NREG  ,ITRN  ,ITST,
     >                  ITSYM(1,ICEL),NM    ,MIX   ,ISPLT ,DAMESH,
     >                  NMS   ,MIXC  ,ISPLTR,DAMESR)
          ENDIF
        ENDIF
*----
*  Final mesh after unfolding
*----
        IF(MOD(ITRN,2) .EQ. 0) THEN
*----
*  First turn geometry type if required
*----
          IF(ITYPG .EQ. 10) THEN
            ITYPG=11
          ELSE IF(ITYPG .EQ.11) THEN
            ITYPG=10
          ELSE IF(ITYPG .EQ.21) THEN
            ITYPG=22
          ELSE IF(ITYPG .EQ.22) THEN
            ITYPG=21
          ENDIF
        ENDIF
*----
*  Then create new mesh
*----
        IF(ITYPG .EQ. 12 .OR. ITYPG .EQ. 13 .OR. 
     >     ITYPG .EQ. 26 .OR. ITYPG .EQ. 27) THEN
          CALL NXTSGT(IPTRK ,IPRINT,MAXMSH,ITYPG ,ICEL  ,ILEV  ,
     >                MAXMSS,NMIX  ,NM    ,MIX   ,DAMESH,ISPLT ,
     >                NMIXS ,NMS   ,DAMESR,
     >                ITSYM(1,ICEL),
     >                NREGS ,NSURS ,NREGN ,NSURN ,NEREN ,
     >                IDREG ,IDSUR )
        ELSE
          CALL NXTSGI(IPTRK ,IPRINT,MAXMSH,ITYPG ,ICEL  ,ILEV  ,
     >                MAXMSS,NMIX  ,NM    ,MIX   ,DAMESH,ISPLT ,
     >                NMIXS ,NMS   ,DAMESR,
     >                ITSYM(1,ICEL),
     >                NREGS ,NSURS ,NREGN ,NSURN ,NEREN ,
     >                IDREG ,IDSUR )
        ENDIF
        ALLOCATE(INDXSR(5,NREGS+NSURS+1))
        INDXSR(:5,:NREGS+NSURS+1)=0
        NBPIN=ICPIN
        MXGSUR=MAX(MXGSUR,NSURS)
        MXGREG=MAX(MXGREG,NREGS)
        IEDIMX(:NSTATE)=0
        IEDIMX(1)=ITYPG
        IEDIMX(2)=NMS(4)
        IEDIMX(3)=NMS(1)
        IEDIMX(4)=NMS(2)
        IEDIMX(5)=NMS(3)
        IEDIMX(6)=NMIXS
        IEDIMX(7)=ILEV
        IEDIMX(8)=NREGS
        IEDIMX(9)=NSURS
        IEDIMX(10)=NREGN
        IEDIMX(11)=NSURN
        IEDIMX(12)=NAREG+1
        CALL NXTVOL(IPTRK ,IPRINT,MAXMSS,ITYPG ,IDIRC ,ICEL  ,
     >              ILEV  ,NMS   ,NREGS ,NSURS ,NREGN ,NSURN ,
     >              MAXPIN,NBPIN ,ITPIN ,DRAPIN,IDREG ,IDSUR ,
     >              DAMESR,INDXSR,NAREG )
*----
*  Renumber cell surfaces
*----
        IEDIMX(13)=NAREG
        IEDIMX(14)=NASUR+1
        CALL NXTRCS(IPTRK ,IPRINT,ICEL ,ILEV  ,
     >              NREGS ,NSURS ,NSURN ,IDFEX(0,ICEL),
     >              INDXSR,NASUR ,IDSUR )
        IEDIMX(15)=NASUR
        IEDIMX(16)=NBGCLS
        IEDIMX(17)=IGCLS-NBGCLS+1
        IEDIMX(18)=NBPIN
        IEDIMX(19)=IDFEX(0,ICEL)
        IEDIMX(20)=IDIRC
        DEALLOCATE(INDXSR)
        WRITE(NAMREC,'(A1,I8.8,A3)') CLEV(ILEV),ICEL,'DIM'
        CALL LCMPUT(IPTRK,NAMREC,NSTATE,1,IEDIMX)
        DEALLOCATE(DAMESR)
        DEALLOCATE(IDSUR,IDREG)
        IF(ILCELL .EQ. 1) THEN
          CALL LCMSIX(IPGEO,NAMGG,2)
        ENDIF
      ENDDO
*----
*  Processing finished:
*  print routine closing output header if required
*  and return
*----
      DEALLOCATE(ISPLTR,NAGCLS,MIXC)
      IF(IPRINT .GE. 10) THEN
        WRITE(IOUT,6001) NAMSBR
      ENDIF
      NFSUR=NASUR
      NFREG=NAREG
*----
*  Scratch storage deallocation
*----
      DEALLOCATE(DRAPIN,DAMESH)
      DEALLOCATE(CPINY,CPINX,APIN,RPIN,RMESH)
      DEALLOCATE(MIX,ISPLT)
      RETURN
*----
*  Output formats
*----
 6000 FORMAT('(* Output from --',A6,'-- follows ')
 6001 FORMAT('   Output from --',A6,'-- completed *)')
 6010 FORMAT(' Processing CELL    = ',I5.4,5X,A12)
 6011 FORMAT(' Processing CLUSTER = ',I5.4,5X,A12)
 6012 FORMAT(' Global cell processing ')
      END