Initial commit from Polytechnique Montreal

1 files changed, 556 insertions, 0 deletions

diff --git a/Dragon/src/NXTGMD.f b/Dragon/src/NXTGMD.f
new file mode 100644
index 0000000..4713e81
--- /dev/null
+++ b/Dragon/src/NXTGMD.f
@@ -0,0 +1,556 @@
+*DECK NXTGMD
+      SUBROUTINE NXTGMD(IPGEO ,IPTRK ,IPRINT,ITYPBC,ILCELL,NEDIM,
+     >                  NBOCEL,NBUCEL,MAXCEL,NUCELL,IUNFLD,
+     >                  IEDIMG,NAGGEO,ITURN ,MERGE ,IDIRR ,
+     >                  DCMESH,DGMESH)
+*
+*----------
+*
+*Purpose:
+* To obtain general dimension vector for geometry.
+* To evaluate global mesh for assembly.
+*
+*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.
+* ITYPBC  type of boundary conditions where:
+*         ITYPBC=0 for geometry with all Cartesian boundaries;
+*         ITYPBC=1 for geometry with one annular boundary;
+*         ITYPBC=2 for geometry with one hexagonal boundary.
+* ILCELL  cell level.
+* NEDIM   number of elements for general dimension vector.
+* NBOCEL  number of cells in original geometry.
+* NBUCEL  number of cells in unfolded geometry.
+* MAXCEL  maximum number of elements in mesh vector for
+*         each directions.
+* NUCELL  number of cell after unfolding in
+*         $X$, $Y$ and $Z$ directions.
+* IUNFLD  array to identify the generating cell (IUNFLD(1,*))
+*         and the rotation associated with this region in space.
+*
+*Parameters: input/output
+* IEDIMG  general dimension vector for problem where
+* NAGGEO  geometry names.
+* ITURN   geometry turns.
+* MERGE   geometry merge.
+* 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.
+* DCMESH  global mesh for each cell.
+* DGMESH  meshing vector for global geometry.
+*
+*Reference:
+*  G. Marleau,
+*  New Geometries Processing in DRAGON: The NXT: Module,
+*  Report IGE-260, Polytechnique Montreal,
+*  Montreal, 2005.
+*  \\\\
+*  This routine is based on the XELDCL routine written by
+*  R. Roy for the EXCELT: module. It contains an additional
+*  level for cluster subgeometry analysis.
+*
+*----------
+*
+      USE              GANLIB
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      TYPE(C_PTR)      IPGEO,IPTRK
+*      INTEGER          IPGEO,IPTRK
+      INTEGER          IPRINT,ITYPBC,ILCELL,NEDIM,NBOCEL,NBUCEL,MAXCEL,
+     >                 NUCELL(3),IUNFLD(2,NBUCEL)
+      INTEGER          IEDIMG(NEDIM),NAGGEO(3,NBOCEL),
+     >                 ITURN(NBOCEL),MERGE(NBOCEL),IDIRR(NBOCEL)
+      DOUBLE PRECISION DCMESH(3,2,NBOCEL),DGMESH(0:MAXCEL,3)
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTGMD')
+      INTEGER          NSTATE
+      PARAMETER       (NSTATE=40)
+      DOUBLE PRECISION DCUT
+      PARAMETER       (DCUT=1.0D-6)
+*----
+*  Local variables
+*----
+      INTEGER          ISTATG(NSTATE)
+      INTEGER          NDIM,NSDIM,NCDIM,IDIRC
+      INTEGER          MAXMSH,MAXREG,NBTCLS,
+     >                 MAXPIN,MAXMSP,MAXRSP,NBGCLS,
+     >                 NTPIN,NPIN,NBSP,NBRSP,NZU
+      INTEGER          ITC,IDIR,ICEL,ICN,ISPL,ICLS,
+     >                 IX,IY,IZ,IPOS,INTRN,IRTRN
+      INTEGER          ILCMLN,ILCMTY
+      CHARACTER        NAMGG*12,NAMREC*12,NAMCL*12
+      REAL             OFFCEN(3)
+      DOUBLE PRECISION DRW(3),DNW(3)
+      REAL             SIDEH
+      DOUBLE PRECISION ARGS
+      INTEGER          IH,IR,IS,ISS,NCH,NCR,NCS,NSS,ISPLTH
+*----
+*  Allocatable arrays
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IGEN,ISPLT
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NAGCEL,NAGCLS
+      REAL, ALLOCATABLE, DIMENSION(:) :: RMESH,XMESH
+*----
+*  Data
+*----
+      CHARACTER        CDIR(4)*1
+      SAVE             CDIR
+      DATA             CDIR /'X','Y','Z','R'/
+*----
+*  Processing starts:
+*  print routine openning output header if required
+*  and initialize various parameters.
+*----
+      IF(IPRINT .GE. 10) THEN
+        WRITE(IOUT,6000) NAMSBR
+      ENDIF
+      NDIM=IEDIMG(1)
+      MAXMSH=IEDIMG(16)
+      MAXREG=1
+      NBTCLS=0
+      MAXPIN=0
+      MAXMSP=0
+      MAXRSP=0
+      DGMESH(0:MAXCEL,:3)=0.0D0
+      IF(ILCELL .EQ. 1) THEN
+        ALLOCATE(IGEN(NBOCEL),NAGCEL(3,NBOCEL))
+        CALL LCMGET(IPGEO,'CELL',NAGCEL)
+        CALL LCMLEN(IPGEO,'MIX',ILCMLN,ILCMTY)
+        IF(ILCMLN .GT. 0) THEN
+          CALL LCMGET(IPGEO,'MIX',IGEN)
+          DO ICEL=1,NBOCEL
+            ICN=-IGEN(ICEL)
+            DO ITC=1,3
+              NAGGEO(ITC,ICEL)=NAGCEL(ITC,ICN)
+            ENDDO
+          ENDDO
+        ELSE
+          CALL LCMGET(IPGEO,'GENERATING',IGEN)
+          DO ICEL=1,NBOCEL
+            ICN=IGEN(ICEL)
+            DO ITC=1,3
+              NAGGEO(ITC,ICEL)=NAGCEL(ITC,ICN)
+            ENDDO
+          ENDDO
+        ENDIF
+        DEALLOCATE(NAGCEL,IGEN)
+        CALL LCMLEN(IPGEO,'TURN',ILCMLN,ILCMTY)
+        IF(ILCMLN .GT. 0) THEN
+          CALL LCMGET(IPGEO,'TURN',ITURN)
+        ELSE
+          ITURN(:NBOCEL)=1
+        ENDIF
+        CALL LCMLEN(IPGEO,'MERGE',ILCMLN,ILCMTY)
+        IF(ILCMLN .GT. 0) THEN
+*----
+*  Treatment of MERGE not yet implemented
+*----
+          CALL XABORT(NAMSBR//
+     >    ': the module NXT cannot process geometry MERGE option')
+*          CALL LCMGET(IPGEO,'MERGE',MERGE)
+        ELSE
+          DO ICEL=1,NBOCEL
+            MERGE(ICEL)=ICEL
+          ENDDO
+        ENDIF
+      ELSE
+        NAMGG='            '
+        READ(NAMGG,'(3A4)') (NAGGEO(ITC,1),ITC=1,3)
+        ITURN(1)=1
+        MERGE(1)=1
+      ENDIF
+      DO ICEL=1,NBOCEL
+        IF(ILCELL .EQ. 1) THEN
+          WRITE(NAMGG,'(3A4)') (NAGGEO(ITC,ICEL),ITC=1,3)
+          CALL LCMSIX(IPGEO,NAMGG,1)
+        ENDIF
+        ISTATG(:NSTATE)=0
+        CALL LCMGET(IPGEO,'STATE-VECTOR',ISTATG)
+        IRTRN=ITURN(ICEL)
+*----
+*  Find cell dimension and direction
+*  Default is CAR3D
+*----
+        IDIRR(ICEL)=0
+        NSDIM=3
+        IF(ISTATG(1) .EQ. 5  .OR.
+     >     ISTATG(1) .EQ. 8  .OR.
+     >     ISTATG(1) .EQ. 12 .OR.
+     >     ISTATG(1) .EQ. 20 .OR.
+     >     ISTATG(1) .EQ. 26 ) THEN
+          IDIRR(ICEL)=1
+          NSDIM=2
+        ELSE IF(ISTATG(1) .EQ. 21) THEN
+          IDIRR(ICEL)=2
+        ELSE IF(ISTATG(1) .EQ. 22) THEN
+          IDIRR(ICEL)=3
+        ELSE IF(ISTATG(1) .EQ. 23 .OR.
+     >     ISTATG(1) .EQ. 27 ) THEN
+          IDIRR(ICEL)=1
+        ENDIF
+        IF(NDIM .NE. NSDIM) CALL XABORT(NAMSBR//
+     >  ': Geometry and sub-geometry must have the same dimensions')
+        IF(ISTATG(1) .EQ. 12 .OR.
+     >     ISTATG(1) .EQ. 13 ) THEN
+          MAXMSH=MAX(MAXMSH,ISTATG(2),2*ISTATG(3)+1,ISTATG(5))
+          MAXREG=MAX(MAXREG,6*ISTATG(3)*ISTATG(3)*MAX(ISTATG(5),1))
+          ALLOCATE(ISPLT(MAXMSH))
+          MAXMSP=2*ISTATG(3)+1
+          NBRSP=6*ISTATG(3)*ISTATG(3)
+          NAMREC='SPLITH      '
+          CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+          IF(ILCMLN .GT. 0) THEN
+            CALL LCMGET(IPGEO,NAMREC,ISPLTH)
+            ISPLTH=ISPLTH*ISTATG(3)
+            MAXMSP=MAX(2*ISPLTH+1,MAXMSP)
+            NBRSP=MAX(6*ISPLTH*ISPLTH,NBRSP)
+          ENDIF
+          IDIR=3
+          NAMREC='SPLIT'//CDIR(IDIR)//'      '
+          CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+          IF(ILCMLN .GT. 0) THEN
+            CALL LCMGET(IPGEO,NAMREC,ISPLT)
+            NBSP=0
+            DO ISPL=1,ILCMLN
+              NBSP=NBSP+ABS(ISPLT(ISPL))
+            ENDDO
+          ELSE
+            NBSP=MAX(ISTATG(5),1)
+          ENDIF
+          NBRSP=NBRSP*NBSP
+          MAXMSP=MAX(MAXMSP,NBSP+1)
+          MAXRSP=MAX(MAXRSP,NBRSP)
+        ELSE IF(ISTATG(1) .EQ. 26 .OR.
+     >          ISTATG(1) .EQ. 27 ) THEN
+          MAXMSH=MAX(MAXMSH,ISTATG(2)+1,2*ISTATG(3)+1,ISTATG(5))
+          MAXREG=MAX(MAXREG,(ISTATG(2)+1)*(6*ISTATG(3)*ISTATG(3))
+     >                      *MAX(ISTATG(5),1))
+          ALLOCATE(ISPLT(MAXMSH))
+          MAXMSP=2*ISTATG(3)+1
+          NBRSP=6*ISTATG(3)*ISTATG(3)
+          NAMREC='SPLITH      '
+          CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+          IF(ILCMLN .GT. 0) THEN
+            CALL LCMGET(IPGEO,NAMREC,ISPLTH)
+            ISPLTH=ISPLTH*ISTATG(3)
+            MAXMSP=MAX(2*ISPLTH+1,MAXMSP)
+            NBRSP=MAX(6*ISPLTH*ISPLTH,NBRSP)
+          ENDIF
+          IDIR=3
+          NAMREC='SPLIT'//CDIR(IDIR)//'      '
+          CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+          IF(ILCMLN .GT. 0) THEN
+            CALL LCMGET(IPGEO,NAMREC,ISPLT)
+            NBSP=0
+            DO ISPL=1,ILCMLN
+              NBSP=NBSP+ABS(ISPLT(ISPL))
+            ENDDO
+          ELSE
+            NBSP=MAX(ISTATG(5),1)
+          ENDIF
+          NBRSP=(ISTATG(2)+1)*NBRSP*NBSP
+          MAXMSP=MAX(MAXMSP,NBSP+1)
+          MAXRSP=MAX(MAXRSP,NBRSP)
+        ELSE
+          MAXMSH=MAX(MAXMSH,ISTATG(2),ISTATG(3),ISTATG(4),ISTATG(5))
+          MAXREG=MAX(MAXREG,(ISTATG(2)+1)*
+     >           MAX(ISTATG(3),1)*MAX(ISTATG(4),1)*MAX(ISTATG(5),1))
+          ALLOCATE(ISPLT(MAXMSH))
+          NBRSP=1
+          DO IDIR=1,4
+            NAMREC='SPLIT'//CDIR(IDIR)//'      '
+            CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+            IF(ILCMLN .GT. 0) THEN
+              CALL LCMGET(IPGEO,NAMREC,ISPLT)
+              NBSP=0
+              DO ISPL=1,ILCMLN
+                NBSP=NBSP+ABS(ISPLT(ISPL))
+              ENDDO
+*            NBRSP=NBRSP*(NBSP+1)
+              IF(IDIR .EQ. 4) THEN
+                NBRSP=NBRSP*(NBSP+1)
+              ELSE
+                NBRSP=NBRSP*NBSP
+              ENDIF
+              MAXMSP=MAX(MAXMSP,NBSP+1)
+            ENDIF
+          ENDDO
+          MAXRSP=MAX(MAXRSP,NBRSP)
+        ENDIF
+        DEALLOCATE(ISPLT)
+*----
+*  Get off center
+*----
+        INTRN=1
+        OFFCEN(:3)=0.0
+        CALL LCMLEN(IPGEO,'OFFCENTER   ',ILCMLN,ILCMTY)
+        IF(ILCMLN .GT. 0)
+     >  CALL LCMGET(IPGEO,'OFFCENTER   ',OFFCEN)
+        DO IDIR=1,3
+          DRW(IDIR)=OFFCEN(IDIR)
+        ENDDO
+        CALL NXTTRM(IRTRN,INTRN,DRW,DCMESH(1,2,ICEL))
+        IF(ITYPBC .EQ. 2) THEN
+*----
+*  Get H mesh
+*----
+          ALLOCATE(RMESH(MAXMSH+1))
+          CALL LCMGET(IPGEO,'SIDE        ',SIDEH)
+          DRW(1)=DBLE(SIDEH)
+*----
+*  Get z-mesh
+*----
+          IDIR=3
+          NAMREC='MESH'//CDIR(IDIR)//'       '
+          CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+          IF(ILCMLN .GT. 0) THEN
+            CALL LCMGET(IPGEO,NAMREC,RMESH)
+            DRW(IDIR)=DBLE(RMESH(ILCMLN))-DBLE(RMESH(1))
+          ENDIF
+          DEALLOCATE(RMESH)
+          CALL NXTTRM(IRTRN,INTRN,DRW,DCMESH(1,1,ICEL))
+*----
+*  Fill DGMESH according to IUNFLD
+*----
+          IPOS=0
+          NZU=MAX(NUCELL(3),1)
+          NCH=NUCELL(1)
+          ARGS=DBLE(12*NCH-3)
+          NCR=(NINT(SQRT(ARGS))+3)/6
+          DO IZ=1,MAX(NUCELL(3),1)
+            IH=0
+            DO IR=1,NCR
+              NCS=6
+              NSS=MAX(1,IR-1)
+              IF(IR.EQ.1) NCS=1
+              DO IS=1,NCS
+                DO ISS=1,NSS
+                  IH=IH+1
+                  IPOS=IPOS+1
+*----
+*  Locate crown (IR), sector (IS) and hexagon number in sector (ISS)
+*----
+                  IF(IUNFLD(1,IPOS) .EQ. ICEL) THEN
+*----
+*  cell is located at this position
+*  Turn width according to symmetry
+*----
+                    INTRN=IUNFLD(2,IPOS)
+                    CALL NXTTRM(IRTRN,INTRN,DRW,DNW(1))
+                    IF(DGMESH(0,1) .EQ. 0.0D0) THEN
+                      DGMESH(0,1)=DRW(1)
+                      DGMESH(0,2)=DRW(1)
+                    ELSE IF(ABS(DGMESH(0,1)-DRW(1)) .GT. DCUT) THEN
+                      CALL XABORT(NAMSBR//
+     >              ': H- Mesh in assembly is not uniform')
+                    ENDIF
+                    CALL NXTHCL(IPRINT,IR,IS,ISS,DRW(1),
+     >                          DGMESH(IH,1),DGMESH(IH,2))
+                    IF(DGMESH(IZ,3) .EQ. 0.0D0) THEN
+                      DGMESH(IZ,3)=DNW(3)
+                    ELSE IF(ABS(DGMESH(IZ,3)-DNW(3)) .GT. DCUT) THEN
+                      CALL XABORT(NAMSBR//
+     >                ': Z Mesh in assembly is not uniform')
+                    ENDIF
+                  ENDIF
+                ENDDO
+              ENDDO
+            ENDDO
+          ENDDO
+        ELSE
+*----
+*  Get Cartesian mesh
+*----
+          ALLOCATE(XMESH(MAXMSH+1))
+          DO IDIR=1,3
+            NAMREC='MESH'//CDIR(IDIR)//'       '
+            CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+            IF(ILCMLN .GT. 0) THEN
+              CALL LCMGET(IPGEO,NAMREC,XMESH)
+              DRW(IDIR)=DBLE(XMESH(ILCMLN))-DBLE(XMESH(1))
+            ENDIF
+          ENDDO
+          DEALLOCATE(XMESH)
+          CALL NXTTRM(IRTRN,INTRN,DRW,DCMESH(1,1,ICEL))
+*----
+*  Fill DGMESH according to IUNFLD
+*----
+          IPOS=0
+          NZU=MAX(NUCELL(3),1)
+          DO IZ=1,MAX(NUCELL(3),1)
+            DO IY=1,NUCELL(2)
+              DO IX=1,NUCELL(1)
+                IPOS=IPOS+1
+                IF(IUNFLD(1,IPOS) .EQ. ICEL) THEN
+*----
+*  cell is located at this position
+*  Turn width according to symmetry
+*----
+                  INTRN=IUNFLD(2,IPOS)
+                  CALL NXTTRM(IRTRN,INTRN,DRW,DNW(1))
+                  IF(DGMESH(IX,1) .EQ. 0.0D0) THEN
+                    DGMESH(IX,1)=DNW(1)
+                  ELSE IF(ABS(DGMESH(IX,1)-DNW(1)) .GT. DCUT) THEN
+                    CALL XABORT(NAMSBR//
+     >              ': X Mesh in assembly is not uniform')
+                  ENDIF
+                  IF(DGMESH(IY,2) .EQ. 0.0D0) THEN
+                    DGMESH(IY,2)=DNW(2)
+                  ELSE IF(ABS(DGMESH(IY,2)-DNW(2)) .GT. DCUT) THEN
+                    CALL XABORT(NAMSBR//
+     >              ': Y Mesh in assembly is not uniform')
+                  ENDIF
+                  IF(DGMESH(IZ,3) .EQ. 0.0D0) THEN
+                    DGMESH(IZ,3)=DNW(3)
+                  ELSE IF(ABS(DGMESH(IZ,3)-DNW(3)) .GT. DCUT) THEN
+                    CALL XABORT(NAMSBR//
+     >              ': Z Mesh in assembly is not uniform')
+                  ENDIF
+                ENDIF
+              ENDDO
+            ENDDO
+          ENDDO
+        ENDIF
+        IF(ISTATG(13) .GT. 0) THEN
+          NBGCLS=ISTATG(13)
+          NBTCLS=NBTCLS+NBGCLS
+          ALLOCATE(NAGCLS(3,NBGCLS))
+          CALL LCMGET(IPGEO,'CLUSTER',NAGCLS)
+          NTPIN=0
+          DO ICLS=1,NBGCLS
+            WRITE(NAMCL,'(3A4)') (NAGCLS(ITC,ICLS),ITC=1,3)
+            CALL LCMSIX(IPGEO,NAMCL,1)
+            ISTATG(:NSTATE)=0
+            CALL LCMGET(IPGEO,'STATE-VECTOR',ISTATG)
+*----
+*  Find cluster dimension and direction
+*----
+            IDIRC=1
+            NCDIM=3
+            IF(ISTATG(1) .EQ. 3) THEN
+              NCDIM=2
+            ELSE IF(ISTATG(1) .EQ. 10) THEN
+              IDIRC=2
+            ELSE IF(ISTATG(1) .EQ. 11) THEN
+              IDIRC=3
+            ELSE IF(ISTATG(1) .EQ. 6) THEN
+              IDIRC=1
+            ELSE
+              CALL XABORT(NAMSBR//': Invalid pin geometry'//
+     >' -- Only TUBE* permitted')
+            ENDIF
+            IF(NSDIM .NE. NCDIM) CALL XABORT(NAMSBR//
+     >': Sub-geometry and clusters must have the same dimensions')
+            IF(IDIRC .NE. IDIRR(ICEL)) THEN
+              IF(IDIRR(ICEL) .NE. 0) CALL XABORT(NAMSBR//
+     >': Pin and geometry direction are incompatible')
+              IDIRR(ICEL)=IDIRC
+            ENDIF
+            MAXMSH=MAX(MAXMSH,ISTATG(2),ISTATG(3),ISTATG(4),ISTATG(5))
+            MAXREG=MAX(MAXREG,ISTATG(6))
+            CALL LCMGET(IPGEO,'NPIN',NPIN)
+            NTPIN=NTPIN+NPIN
+            ALLOCATE(ISPLT(MAXMSH))
+            NBRSP=1
+            DO IDIR=1,4
+              NAMREC='SPLIT'//CDIR(IDIR)//'      '
+              CALL LCMLEN(IPGEO,NAMREC,ILCMLN,ILCMTY)
+              IF(ILCMLN .GT. 0) THEN
+                CALL LCMGET(IPGEO,NAMREC,ISPLT)
+                NBSP=0
+                DO ISPL=1,ILCMLN
+                  NBSP=NBSP+ABS(ISPLT(ISPL))
+                ENDDO
+                IF(IDIR .EQ. 4) THEN
+                  NBRSP=NBRSP*(NBSP+1)
+                ELSE
+                  NBRSP=NBRSP*NBSP
+                ENDIF
+                MAXMSP=MAX(MAXMSP,NBSP+1)
+              ENDIF
+            ENDDO
+            MAXRSP=MAX(MAXRSP,NBRSP)
+            DEALLOCATE(ISPLT)
+            CALL LCMSIX(IPGEO,NAMCL,2)
+          ENDDO
+          MAXPIN=MAX(MAXPIN,NTPIN)
+*----
+*  Now test if pin overlapp
+*----
+          CALL NXTTPO(IPGEO ,IPRINT,ITYPBC,NBGCLS,NTPIN ,MAXMSH,NCDIM ,
+     >                IDIRR(ICEL) ,DRW   ,OFFCEN,NAGCLS)
+          DEALLOCATE(NAGCLS)
+        ENDIF
+        IF(IDIRR(ICEL) .EQ. 0) IDIRR(ICEL)=1
+        IF(ILCELL .EQ. 1) THEN
+          CALL LCMSIX(IPGEO,NAMGG,2)
+        ENDIF
+      ENDDO
+      IEDIMG(16)=MAXMSH
+      IEDIMG(17)=MAXREG
+      IEDIMG(18)=NBTCLS
+      IEDIMG(19)=MAX(MAXPIN,1)
+      IEDIMG(20)=MAXMSP+1
+      IEDIMG(21)=MAXRSP
+      IF(ITYPBC .EQ. 2) THEN
+        IF(NUCELL(1) .GT. 0) THEN
+          DO IDIR=1,2
+            NAMREC='G00000001SM'//CDIR(IDIR)
+            CALL LCMPUT(IPTRK,NAMREC,(NUCELL(1)+1),4,
+     >                  DGMESH(0,IDIR))
+          ENDDO
+        ENDIF
+        IDIR=3
+        IF(NUCELL(IDIR) .GT. 0) THEN
+          DO IX=1,NUCELL(IDIR)
+            DGMESH(IX,IDIR)=DGMESH(IX-1,IDIR)+DGMESH(IX,IDIR)
+          ENDDO
+          NAMREC='G00000001SM'//CDIR(IDIR)
+          CALL LCMPUT(IPTRK,NAMREC,(NUCELL(IDIR)+1),4,
+     >                DGMESH(0,IDIR))
+        ENDIF
+      ELSE
+        DO IDIR=1,3
+          IF(NUCELL(IDIR) .GT. 0) THEN
+            DO IX=1,NUCELL(IDIR)
+              DGMESH(IX,IDIR)=DGMESH(IX-1,IDIR)+DGMESH(IX,IDIR)
+            ENDDO
+            NAMREC='G00000001SM'//CDIR(IDIR)
+            CALL LCMPUT(IPTRK,NAMREC,(NUCELL(IDIR)+1),4,
+     >                  DGMESH(0,IDIR))
+          ENDIF
+        ENDDO
+      ENDIF
+      CALL LCMPUT(IPTRK,'G00000001DIM',NEDIM,1,IEDIMG)
+*----
+*  Processing finished:
+*  print routine closing output header if required
+*  and return
+*----
+      IF(IPRINT .GE. 10) THEN
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+      END