C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LDMGE            **
C**************************
C
C
       IMPLICIT NONE
C
C Purpose: Create .gen/.ge1/.ge2 files from an LDM file
C
C Author:  John W. Campbell, October 1993
C
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C
      LOGICAL CCPEXS
      INTEGER NXTTOK
      INTEGER INTFP
C
C====== LOCALS:
C
      LOGICAL QUIT
      LOGICAL EXISTS
      LOGICAL EOF
      LOGICAL NONLDM
      LOGICAL DISTOR
      LOGICAL WR_SPAT
      LOGICAL WR_CLOS
      INTEGER IERR
      INTEGER IFAIL
      INTEGER IUN_LDM
      INTEGER IUN_GEN
      INTEGER IUN_GE1
      INTEGER IUN_GE2
      INTEGER IEND
      INTEGER INDIRECT
      INTEGER KFLAG
      INTEGER MININT
      INTEGER PACK_ID
      INTEGER IFLAG
      INTEGER KEYPLATE
      INTEGER NUM_SPOTS
      INTEGER NUM_MULTIPLETS
      INTEGER NUM_SPATIALS
      INTEGER NUM_NODALS
      INTEGER NUM_BOTH
      INTEGER I
      INTEGER IPACK
      INTEGER NUMPACK
      CHARACTER*120 FILNAM
      CHARACTER*120 PAR_FILE
      CHARACTER*120 GEN_FILE
      CHARACTER*120 GE1_FILE
      CHARACTER*120 GE2_FILE
      CHARACTER*120 STR
      CHARACTER*80 BADTOK
      CHARACTER*80 ERRSTR
      CHARACTER*250 LIN

C     LOGICAL QUIT            !Quit program flag
C     LOGICAL EXISTS          !File exists flag
C     INTEGER I               !Temp/loop variable
C     INTEGER IPACK           !Pack number
C     INTEGER NUMPACK         !Number of packs
C     INTEGER J               !Temp/loop variable
C     INTEGER IERR            !Error flag
C     INTEGER IFAIL           !File open error flag
C     INTEGER IL              !Character counter for NXTTOK routine
C     INTEGER M1              !Start of field found by NXTTOK
C     INTEGER M2              !End of field found by NXTTOK
C     INTEGER IV              !Integer value returned from INTFP
C     CHARACTER*120 FILOUT    !Output file name
C     CHARACTER*120 GEN_FILE  !.gen file name
C     CHARACTER*120 STR       !Record read from .gen file
C     CHARACTER*2 ERRSTR      !Error string (dummy)
C

C
C-------------------------------------------------------------------------------
C
C====== Initialisations
C
      IUN_LDM = 1
      IUN_GEN = 2
      IUN_GE1 = 3
      IUN_GE2 = 4
      CALL TRWSET (5, 6, 0, 0, 0)
C
C====== Read in LDM file name
C
100   CALL TPROMP ('LDM input file name (default ext=.ldm): ')       
      CALL TREAD (FILNAM, QUIT)
      IF (FILNAM.EQ.' ') STOP
      CALL FILEXT (FILNAM, 'ldm', PAR_FILE, I)
      EXISTS = CCPEXS (PAR_FILE)
      IF (.NOT.EXISTS) THEN
         CALL TWRITE ('**File not found**')
         GO TO 100
      END IF
      IFAIL = 1
      CALL CCPDPN (-IUN_LDM, PAR_FILE, 'READONLY', 'F', 0, IFAIL)
      IF (IFAIL.LT.0) THEN
         CALL TWRITE ('**Error in opening file**')
         GO TO 100
      END IF
C
C====== Reset default LDM parameters
C
      CALL LDM_RESET (' ', 0, 0, IERR, ERRSTR)
C
C====== Read in parameters file
C
      NONLDM = .FALSE.
      INDIRECT = 0
150   CALL LDM_READLINE (IUN_LDM, 10, INDIRECT, LIN, EOF, IERR)
      IF (EOF) GO TO 200
      IF (IERR.GT.0) THEN
         STR = '**Error** LDM_READLINE error, ierr='
         CALL STR_ADDI (STR,IERR,0)
         CALL TWRITE (STR)
         GO TO 200
      END IF
      CALL LDM_PARSE (LIN, KFLAG, BADTOK, ERRSTR)
      IF (KFLAG.LT.0) THEN
         CALL TWRITE ('**Error** LDM_PARSE error')
         STR = 'Line:'
         CALL STR_ADDS (STR, LIN, 1)
         CALL TWRITE (STR)
         STR = 'Bad token:'
         CALL STR_ADDS (STR, BADTOK, 1)
         CALL TWRITE (STR)
         STR = 'Error message:'
         CALL STR_ADDS (STR, ERRSTR, 1)
         CALL TWRITE (STR)
         GO TO 150
      END IF
      IF (KFLAG.EQ.1) NONLDM = .TRUE.
      GO TO 150
C
C====== End of file (or fatal error)
C
200   IF (NONLDM) CALL TWRITE ('*Warning* Non-LDM data present')
C
C====== Get .gen/.ge1/.ge2 file names
C
250   CALL TPROMP ('Root name for .gen/.ge1/.ge2 files: ')
      CALL TREAD (FILNAM, QUIT)
      IF (FILNAM.EQ.' ') STOP
      CALL FILEXT (FILNAM, 'gen', GEN_FILE, I)
      EXISTS = CCPEXS (GEN_FILE)
      IF (EXISTS) THEN
         CALL TWRITE ('**.gen file already exists**')
         GO TO 250
      END IF
      CALL FILEXT (FILNAM, 'ge1', GE1_FILE, I)
      EXISTS = CCPEXS (GE1_FILE)
      IF (EXISTS) THEN
         CALL TWRITE ('**.ge1 file already exists**')
         GO TO 250
      END IF
      CALL FILEXT (FILNAM, 'ge2', GE2_FILE, I)
      EXISTS = CCPEXS (GE2_FILE)
      IF (EXISTS) THEN
         CALL TWRITE ('**.ge2 file already exists**')
         GO TO 250
      END IF
      IFAIL = 1
      CALL CCPDPN (-IUN_GEN, GEN_FILE,'NEW','F',0,IFAIL)
      IF (IFAIL.LT.0) THEN
         CALL TWRITE ('**Error in opening the .gen file**')
         GO TO 250
      END IF
      CALL WRITE_GENFIL (IUN_GEN)
      CLOSE (UNIT=IUN_GEN)
      IFAIL = 1
      CALL CCPDPN (-IUN_GE1, GE1_FILE,'NEW','DU',19,IFAIL)
      IF (IFAIL.LT.0) THEN
         CALL TWRITE ('**Error in opening the .ge1 file**')
         GO TO 250
      END IF
      IFAIL = 1
      CALL CCPDPN (-IUN_GE2, GE2_FILE,'NEW','DU',2,IFAIL)
      IF (IFAIL.LT.0) THEN
         CALL TWRITE ('**Error in opening the .ge2 file**')
         GO TO 250
      END IF
      CALL LDM_GET_NUMP (NUMPACK, IFLAG)
 260  IPACK = 1
      CALL TASKI ('Pack number for which .ge files are required',
     +            IPACK, QUIT)
      IF (QUIT) STOP
      IF (IPACK.LT.1.OR.IPACK.GT.NUMPACK) THEN
         CALL TWRITE('**Error** Pack number out of range')
         GO TO 260
      END IF
      WR_SPAT = .FALSE.
      CALL TASKQ ('Do you want to write spatials to the file',
     +            WR_SPAT, QUIT)
      WR_CLOS = .FALSE.
      MININT = 100
      CALL TASKI (
     +     'Enter Min I to by-pass measurements on subsequent plates',
     +      MININT, QUIT)
      DISTOR = .FALSE.
      CALL LDM_GET_PID (1, PACK_ID, IFLAG, IERR, ERRSTR)
      CALL LDM_GET_KEYP (1, KEYPLATE, IFLAG, IERR, ERRSTR)
      CALL WRITE_GEFILE (IPACK, IUN_GE1, IUN_GE2, WR_SPAT, WR_CLOS,
     +                         DISTOR, MININT, PACK_ID, KEYPLATE,
     +                         NUM_SPOTS, NUM_MULTIPLETS,
     +                         NUM_SPATIALS, NUM_BOTH, NUM_NODALS)
      WRITE (STR,6060) NUM_SPOTS
      CALL TWRITE (STR)
      WRITE (STR,6070) NUM_MULTIPLETS
      CALL TWRITE (STR)
      WRITE (STR,6080) NUM_SPATIALS
      CALL TWRITE (STR)
      WRITE (STR,6090) NUM_BOTH
      CALL TWRITE (STR)
      WRITE (STR,6100) NUM_NODALS
      CALL TWRITE (STR)
      STOP
C
C====== Format statements
C
 6060 FORMAT ('Summary:',I6,' spots output, of which')
 6070 FORMAT (8X,I6,' were multiplets,')
 6080 FORMAT (8X,I6,' were spatial overlaps,')
 6090 FORMAT (8X,I6,' were both, and')
 6100 FORMAT (8X,I6,' were nodals.')
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     WRITE_GENFIL     **
C**************************
C
C
      SUBROUTINE WRITE_GENFIL (IUN_GEN)
C
      IMPLICIT NONE
C
C Purpose: Write a .gen file 
C
C Author:  John W. Campbell, February 1992
C
C Note: File must be opened and closed outside this routine
C
C ****** Modified version from LAUEGEN for LDMGE program *****
C
C 
      INTEGER IUN_GEN
C
C IUN_GEN   (R)  Unit no. for .gen file
C
C====== Parameter statements
C 
C      INCLUDE 'limits.finc'
C
C====== GLOBALS:
C
C      INCLUDE 'view_handles.finc'
C      INCLUDE 'file_units.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER LENSTR
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NCP(7)
      INTEGER IFLAG
      INTEGER ND
      INTEGER IERR
      INTEGER ISYST
      INTEGER LATT
      INTEGER IBEAM
      INTEGER IROT
      INTEGER DISTOR_TYPE
      INTEGER NFID
      INTEGER N1OD
      INTEGER NXRAST
      INTEGER NYRAST
      INTEGER IFLIP
      INTEGER NSPDX
      INTEGER NSPDY
      REAL SPDXY
      REAL SPDX_MIN
      REAL SPDX_MAX
      REAL SPDY_MIN
      REAL SPDY_MAX
      REAL SPDX(100)
      REAL SPDY(100)
      REAL CELL(6)
      REAL RCELL(6)
      REAL CELLC(6)
      REAL PHIS(3)
      REAL SPINDLE
      REAL WMIN
      REAL WMAX
      REAL CTOF
      REAL DMIN
      REAL FILMR
      REAL X_C
      REAL Y_C
      REAL W_C
      REAL TWIST
      REAL TILT
      REAL BULGE
      REAL ROFF
      REAL TOFF
      REAL Y_SCALE
      REAL RMIN
      REAL FIDX1
      REAL FIDY1
      REAL FIDX2
      REAL FIDY2
      REAL FIDX3
      REAL FIDY3
      REAL FIDBOX
      REAL G1OD
      REAL BASEOD
      REAL CURV
      REAL XCENF
      REAL YCENF
      REAL RASTER
      REAL XLOW
      REAL XHIGH
      REAL YLOW
      REAL YHIGH
      CHARACTER*80 CELINF
      CHARACTER*12 FLIP(3)
      CHARACTER*5 SYSTYP(7)
      CHARACTER*4 AXTYP(7)
      CHARACTER*1 LATT_STR(7)
      CHARACTER*2 ERRSTR
      CHARACTER*80 STR

C     INTEGER I                     !Temp/loop variable
C     INTEGER NCP(7)                !No. of unique cell parameters
C     INTEGER IFLAG                 !Parameter status flag
C     INTEGER ND                    !No. of decimal places flag
C     INTEGER IERR                  !Error flag
C     INTEGER ISYST                 !Crystal system
C     INTEGER LATT                  !Lattice type
C     INTEGER IBEAM                 !Beam axis
C     INTEGER IROT                  !Rotation axis
C     INTEGER DISTOR_TYPE           !Distortion type
C     INTEGER NFID                  !No. of fiducials
C     INTEGER N1OD                  !Film n1od
C     INTEGER NXRAST                !No. x-rasters
C     INTEGER NYRAST                !No. y-rasters
C     INTEGER IFLIP                 !Film/image-plate flag
C     INTEGER NSPDX                 !Polynomial order spd correction
C     INTEGER NSPDY                 !Polynomial order spd correction
C     REAL SPDXY                    !xy spd correction
C     REAL SPDX_MIN                 !Chebyshev minimum for x
C     REAL SPDX_MAX                 !Chebyshev maximum for x
C     REAL SPDY_MIN                 !Chebyshev minimum for y
C     REAL SPDY_MAX                 !Chebyshev maximum for y
C     REAL SPDX(100)                !Polynomial spd correction
C     REAL SPDY(100)                !Polynomial spd correction
C     REAL CELL(6)                  !Real cell
C     REAL RCELL(6)                 !Reciprocal cell
C     REAL CELLC(6)                 !Compressed cell
C     REAL PHIS(3)                  !Missetting angles
C     REAL SPINDLE                  !Spindle angle
C     REAL WMIN                     !Lambda-min
C     REAL WMAX                     !Lambda-max
C     REAL CTOF                     !Crystal-film distance
C     REAL DMIN                     !'dmin'
C     REAL FILMR                    !Film radius
C     REAL X_C                      !Camera constant x_c
C     REAL Y_C                      !Camera constant y_c
C     REAL W_C                      !Camera constant w_c
C     REAL TWIST                    !Twist
C     REAL TILT                     !Tilt
C     REAL BULGE                    !Bulge
C     REAL ROFF                     !MAR Roff
C     REAL TOFF                     !MAR Toff
C     REAL Y_SCALE                  !'y_scale'
C     REAL RMIN                     !'rmin'
C     REAL FIDX1                    !Fiducials
C     REAL FIDY1                    ! "
C     REAL FIDX2                    ! "
C     REAL FIDY2                    ! "
C     REAL FIDX3                    ! "
C     REAL FIDY3                    ! "
C     REAL FIDBOX                   !Fiducials box search size
C     REAL G1OD                     !Film g1od
C     REAL BASEOD                   !Film baseod
C     REAL CURV                     !Film non-linearity
C     REAL XCENF                    !Image x-centre
C     REAL YCENF                    !Image y-centre
C     REAL RASTER                   !Raster size
C     REAL XLOW                     !Area limit xlow
C     REAL XHIGH                    !Area limit xhigh
C     REAL YLOW                     !Area limit ylow
C     REAL YHIGH                    !Area limit yhigh
C     CHARACTER*80 CELINF           !Cell information string
C     CHARACTER*12 FLIP(3)          !Film/image-plate flags
C     CHARACTER*5 SYSTYP(7)         !Crystal system strings
C     CHARACTER*4 AXTYP(7)          !Axis strings
C     CHARACTER*1 LATT_STR(7)       !Lattice type strings
C     CHARACTER*2 ERRSTR            !Error string (dummy)
C     CHARACTER*80 STR              !Temp string

      DATA NCP/6,4,3,2,2,2,1/
      DATA FLIP /'film','image-plate',' '/
      DATA SYSTYP/'TRICL','MONOC','ORTHO','TETRA','HEXAG','RHOMB',
     +            'CUBIC'/
      DATA AXTYP/'-cst','-bst','-ast',' ','+ast','+bst','+cst'/
      DATA LATT_STR /'P','A','B','C','I','F','R'/
C
C-------------------------------------------------------------------------------
C
      CALL LDM_GET_TITLE (STR, IFLAG)
      WRITE (IUN_GEN,'(A)') STR(1:LENSTR(STR))
      CALL LDM_GET_CRYST (STR, IFLAG)
      WRITE (IUN_GEN,'(A)') STR(1:LENSTR(STR))
      CALL LDM_GET_SYST (STR, ISYST, IFLAG)
      CALL LDM_GET_LATT (STR, LATT, IFLAG)
      CALL LDM_GET_ROT (STR, IROT, IFLAG)
      CALL LDM_GET_BEAM (STR, IBEAM, IFLAG)
      CALL LDM_GET_A (CELL(1), IFLAG, ND)
      CALL LDM_GET_B (CELL(2), IFLAG, ND)
      CALL LDM_GET_C (CELL(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELL(4), IFLAG, ND)
      CALL LDM_GET_BETA (CELL(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELL(6), IFLAG, ND)
      CELINF = ' '
      CALL RECCEL(RCELL,CELL,1.0)
      CALL CELLCX(ISYST,CELLC,RCELL,2,2)
      IF (ISYST.EQ.1) THEN
         WRITE (CELINF(1:60),'(3F10.7,3F10.5)')
     +         (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.2) THEN
         WRITE (CELINF(1:60),'(3F10.7,F10.5)')
     +         (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.3) THEN
         WRITE (CELINF(1:60),'(3F10.7)') (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.4) THEN
         WRITE (CELINF(1:60),'(2F10.7)') (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.5) THEN
         WRITE (CELINF(1:60),'(2F10.7)') (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.6) THEN
         WRITE (CELINF(1:60),'(F10.5,F10.5)')
     +         (CELLC(I),I=1,NCP(ISYST))
      ELSE IF (ISYST.EQ.7) THEN
         WRITE (CELINF(1:60),'(F10.7)') (CELLC(I),I=1,NCP(ISYST))
      END IF
      CELINF = CELINF(1:60)//' '//SYSTYP(ISYST)
     +         //' '//LATT_STR(LATT)
     +         //' '//AXTYP(IROT+4)//' '//AXTYP(IBEAM+4)
      WRITE (IUN_GEN,'(A)') CELINF
      CALL LDM_GET_PHIX (1, PHIS(1), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIY (1, PHIS(2), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIZ (1, PHIS(3), IFLAG, ND, IERR, ERRSTR)
      WRITE (IUN_GEN,*) PHIS
      CALL LDM_GET_LMIN (1, WMIN, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_LMAX (1, WMAX, IFLAG, ND, IERR, ERRSTR)
      WRITE (IUN_GEN,*) 1.0, WMIN , WMAX
      CALL LDM_GET_CTOF (1, 1, CTOF, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_DMIN (1, DMIN, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_RMAX (FILMR, IFLAG, ND)
      WRITE (IUN_GEN,*) CTOF , DMIN , FILMR
      CALL LDM_GET_SPIN (1, SPINDLE, IFLAG, ND, IERR, ERRSTR)
      WRITE (IUN_GEN,*) SPINDLE
      CALL LDM_GET_XC (1, 1, X_C, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_YC (1, 1, Y_C, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_WC (1, 1, W_C, IFLAG, ND, IERR, ERRSTR)
      WRITE (IUN_GEN,*) X_C, Y_C, W_C
      CALL LDM_GET_DSTOR(STR, DISTOR_TYPE, IFLAG)
      CALL LDM_GET_TWIST (1, 1, TWIST, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_TILT (1, 1, TILT, IFLAG, ND, IERR, ERRSTR)
      IF (DISTOR_TYPE.EQ.1) THEN
         CALL LDM_GET_BULGE (1, 1, BULGE, IFLAG, ND, IERR, ERRSTR)
         WRITE (IUN_GEN,*) TWIST, TILT, BULGE
      ELSE IF (DISTOR_TYPE.EQ.2) THEN
         CALL LDM_GET_ROFF (1, 1, ROFF, IFLAG, ND, IERR, ERRSTR)
         CALL LDM_GET_TOFF (1, 1, TOFF, IFLAG, ND, IERR, ERRSTR)
         WRITE (IUN_GEN,*) TWIST, TILT, ROFF, TOFF
      ELSE
         WRITE (IUN_GEN,*) TWIST, TILT, 0.0
      END IF
      CALL LDM_GET_YSCAL (1, 1, Y_SCALE, IFLAG, ND, IERR, ERRSTR)
      WRITE (IUN_GEN,*) Y_SCALE
      CALL LDM_GET_RMIN (RMIN, IFLAG, ND)
      WRITE (IUN_GEN,*) RMIN, 0.0
C
C====== Convert fid posns to original units, uncentred.
C
      CALL LDM_GET_FIDT (STR, NFID, IFLAG)
      IF (NFID.EQ.1) NFID = 3
      IF (NFID.EQ.3) THEN
         CALL LDM_GET_FIDX1 (FIDX1, IFLAG, ND)
         CALL LDM_GET_FIDY1 (FIDX1, IFLAG, ND)
         CALL LDM_GET_FIDX2 (FIDX2, IFLAG, ND)
         CALL LDM_GET_FIDY2 (FIDX2, IFLAG, ND)
         CALL LDM_GET_FIDX3 (FIDX3, IFLAG, ND)
         CALL LDM_GET_FIDY3 (FIDX3, IFLAG, ND)
         WRITE (IUN_GEN,'(I2,3(2F8.1))') NFID, FIDX1, FIDY1,
     +                     FIDX2, FIDY2, FIDX3, FIDY3
      ELSE
         WRITE (IUN_GEN,'(I2,3(2F8.1))') NFID, 0.0, 0.0, 0.0, 0.0, 
     +                                   0.0, 0.0
      END IF
      CALL LDM_GET_FBOX (FIDBOX, IFLAG, ND)
      CALL LDM_GET_N1OD (N1OD, IFLAG)
      CALL LDM_GET_G1OD (G1OD, IFLAG, ND)
      CALL LDM_GET_BASE (BASEOD, IFLAG, ND)
      CALL LDM_GET_NONL (CURV, IFLAG, ND)
      WRITE (IUN_GEN,'(F6.1,F4.1,I4,3E10.3)') 0.0, FIDBOX, N1OD, 
     +       BASEOD, G1OD, CURV
      CALL LDM_GET_NXRAS (NXRAST, IFLAG)
      CALL LDM_GET_NYRAS (NYRAST, IFLAG)
      CALL LDM_GET_XCENF (1, 1, XCENF, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_YCENF (1, 1, YCENF, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_RAST (RASTER, IFLAG, ND)
      CALL LDM_GET_XLOW (XLOW, IFLAG, ND)
      CALL LDM_GET_XHIGH (XHIGH, IFLAG, ND)
      CALL LDM_GET_YLOW (YLOW, IFLAG, ND)
      CALL LDM_GET_YHIGH (YHIGH, IFLAG, ND)
      CALL LDM_GET_IMTYP (STR, IFLIP, IFLAG)
      IF (IFLIP.GT.2) IFLIP = 3
      WRITE (IUN_GEN,'(2I5,F4.1,2F8.2,F7.2,4F7.1,X,A)') NXRAST, NYRAST, 
     +       2.0, XCENF, YCENF, RASTER,
     +       XLOW, XHIGH, YLOW, YHIGH, FLIP(IFLIP)
      IF (DISTOR_TYPE.EQ.3) THEN
         CALL LDM_SPDXY (1, 1, SPDXY,
     +                   NSPDX, SPDX_MIN, SPDX_MAX, SPDX,
     +                   NSPDY, SPDY_MIN, SPDY_MAX, SPDY,
     +                   100, IERR, ERRSTR)
         WRITE (IUN_GEN, '(2F10.2,I3,/(4G14.5)))') SPDX_MIN, SPDX_MAX,
     +          NSPDX, (SPDX(I),I=1,MIN(NSPDX,100))
         WRITE (IUN_GEN, '(2F10.2,I3,/(4G14.5)))') SPDY_MIN, SPDY_MAX,
     +          NSPDY, (SPDY(I),I=1,MIN(NSPDY,100))
         WRITE (IUN_GEN, '(G14.5)') SPDXY
      END IF
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     WRITE_GEFILE     **
C**************************
C
C
      SUBROUTINE WRITE_GEFILE (IPACK, IUN_GE1, IUN_GE2, WR_SPAT, 
     +                         WR_CLOS, DISTOR, MININT, PACK_ID, 
     +                         KEYPLATE, NUM_SPOTS, NUM_MULTIPLETS,
     +                         NUM_SPATIALS, NUM_BOTH, NUM_NODALS)
C
C      INCLUDE 'implicit_none.finc'
C
C Purpose: Write .ge1/.ge2 files
C
C Author:  John W. Campbell, February 1992
C
C Note:    Output .ge1/.ge2 files must be opened before calling this
C          routine. They will be closed within this routine.
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IUN_GE1
      INTEGER IUN_GE2
      LOGICAL WR_SPAT
      LOGICAL WR_CLOS
      LOGICAL DISTOR
      INTEGER MININT
      INTEGER PACK_ID
      INTEGER KEYPLATE
      INTEGER NUM_SPOTS
      INTEGER NUM_MULTIPLETS
      INTEGER NUM_SPATIALS
      INTEGER NUM_NODALS
      INTEGER NUM_BOTH
C
C IPACK          (R)    Pack no. for which .ge files to be written
C IUN_GE1        (R)    Unit no. for .ge1 file (already opened)
C IUN_GE2        (R)    Unit no. for .ge2 file (already opened)
C WR_SPAT        (R)    Write spatials flag
C WR_CLOS        (R)    Write too close spots flag
C DISTOR         (R)    Apply distortion corrections to coordinates flag
C MININT         (R)    Minimum intensity for subsequent films integration
C PACK_ID        (R)    Pack identifier
C KEYPLATE       (R)    Key plate no.
C NUM_SPOTS      (W)    No. of spots written
C NUM_MULTIPLETS (W)    No. of multiplets
C NUM_SPATIALS   (W)    No. of spatials
C NUM_NODALS     (W)    No. nodals
C NUM_BOTH       (W)    No. both spatial & multiple
C
C====== Parameter statements
C 
C
C====== LOCALS:
C
      LOGICAL    SPAT, CLOS, NODFLG
      INTEGER    NODAL_INDEX
      INTEGER    CURRENT_REFLEX, I, J
      INTEGER    II, JJ, KK
      INTEGER    IERR, KERR, IFLAG, ND
      INTEGER    ISPOT, NUMSPOTS
      INTEGER    IBEAM, IROT, DISTOR_TYPE
      INTEGER    LINK_REC
      INTEGER    BINTENS(6), BSIGS(6), PINTENS(6), PSIGS(6)
      INTEGER    M_FLAG, H, K, L, NODIDX, H2, K2, L2
      INTEGER    IREAL
      INTEGER    IH,IK,IL
      INTEGER    GE1_REC, GE2_REC
      INTEGER    MULTIPLICITY, IX, IY, HKLPAK
      INTEGER    JPACK, NFILM
      INTEGER    MULTI_BINS(21)
      INTEGER    BUFFER(19), BUFFER2(2), MM
      INTEGER    MHKL(4), HKL2(4)
      INTEGER    NODAL_BINS(12)
      INTEGER    BOX(5)
      INTEGER    LAMBDA_BINS(16)
      INTEGER    ISYLAT(2), ISL
      INTEGER    IFAIL, NSING, NSINGOV
      INTEGER    ISYST, LATT
      INTEGER    H_N, K_N, L_N, MINH, MAXH, INCH
      REAL       ALNOD, DTHRNOD
      REAL       CTOF_RATS(6), CTOF
      REAL       PHIS(3), SPINDLE, CTOF_KEYPLATE, DELTA, EPSILON
      REAL       X_C(6), Y_C(6), W_C(6), Y_SCALE(6)
      REAL       TWIST(6), TILT(6), BULGE(6)
      REAL       WMIN, WMAX, DMIN, FILMR
      REAL       RINT
      REAL       BX, BY, BDMIN2
      REAL       LAMBDA1, LAMBDA2
      REAL       X_CEN_FS(6)
      REAL       Y_CEN_FS(6)
      REAL       W_FS(6)
      REAL       RCELL(6), CELLC(6), CELL(6)
      REAL       MM_RAST_X, MM_RAST_Y
      CHARACTER*80 STR, COMPOUND
      CHARACTER*6 CRYST
      CHARACTER*2 VALSTR, ERRSTR

      EQUIVALENCE (IREAL, RINT)
      EQUIVALENCE (BUFFER(1),BX)
      EQUIVALENCE (BUFFER(2),BY)
      EQUIVALENCE (BUFFER(16),LAMBDA1)
      EQUIVALENCE (BUFFER(18),LAMBDA2)
      EQUIVALENCE (BUFFER(19),BDMIN2)

      EQUIVALENCE (MHKL(1),M_FLAG)
      EQUIVALENCE (MHKL(2),H)
      EQUIVALENCE (MHKL(3),K)
      EQUIVALENCE (MHKL(4),L)

      EQUIVALENCE (HKL2(1),NODIDX)
      EQUIVALENCE (HKL2(2),H2)
      EQUIVALENCE (HKL2(3),K2)
      EQUIVALENCE (HKL2(4),L2)

      EQUIVALENCE (BUFFER2(2),LAMBDA1)

      DATA  BINTENS, BSIGS, PINTENS, PSIGS /6*-9999, 6*-9999,
     +                6*-9999, 6*-9999/
C
C------------------------------------------------------------------------------
C
C
C====== Put initialised intensity/sigma values in reflection record buffer
C
      CALL CCPII2(BUFFER(3),BINTENS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(6),BSIGS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(9),PINTENS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(12),PSIGS,6,.TRUE.,.FALSE.)
C
C====== Zero histograms
C
      NUM_SPOTS = 0
      NUM_SPATIALS = 0
      NUM_MULTIPLETS = 0
      NUM_BOTH = 0
      NUM_NODALS = 0
      DO 20 I = 1 , 21
         MULTI_BINS(I) = 0
 20   CONTINUE
      NSING = 0
      NSINGOV = 0
      DO 30 I = 1 , 12
         NODAL_BINS(I) = 0
 30   CONTINUE
C
C====== Now write .ge1/.ge2 files 
C
      GE1_REC = 10
      GE2_REC = 0

      CALL LRL_GEN (IPACK, 1, DISTOR, 2, IERR, ERRSTR)
      CALL LRL_NUMSPOTS (NUMSPOTS)
      write(6,*) 'numspots = ', numspots

      DO 40 ISPOT = 1 , NUMSPOTS

         CALL LRL_GET (ISPOT, H, K, L, BX, BY, LAMBDA1, BDMIN2,
     +                 MULTIPLICITY, NODFLG, SPAT, KERR)
         CALL LRL_GET_XFYFD (ISPOT, BX, BY, KERR)
         CALL LRL_GET_OVLP (ISPOT, SPAT, CLOS, KERR)
         M_FLAG = 0
         H2 = 0
         K2 = 0
         L2 = 0
         CALL STBITS (M_FLAG,0,5,MULTIPLICITY)
         IF (SPAT.OR.CLOS) CALL STBITS (M_FLAG,6,1,1)
         IF (CLOS) CALL STBITS (M_FLAG,7,1,1)
         IF (NODFLG) CALL STBITS (M_FLAG,5,1,1)
         IF ((.NOT.SPAT .OR. WR_SPAT) .AND. 
     +       (.NOT.CLOS .OR. WR_CLOS)) THEN
            IF (MULTIPLICITY.EQ.1) THEN
C
C====== Now write out .ge1 record for the single
C
               GE1_REC = GE1_REC + 1
               IF (NODFLG) THEN
                  NUM_NODALS = NUM_NODALS + 1
                  IH = H
                  IK = K
                  IL = L
                  NODIDX = NODAL_INDEX(IH,IK,IL)
                  NODAL_BINS(NODIDX) = NODAL_BINS(NODIDX) + 1
               ELSE
                  NODIDX = 0
               END IF
               LAMBDA2 = 0.0
               CALL CCPIBY(BUFFER(15),MHKL,-4)
               CALL CCPIBY(BUFFER(17),HKL2,-4)
               WRITE (IUN_GE1,REC=GE1_REC) BUFFER

               MULTI_BINS(1) = MULTI_BINS(1) + 1
               NSING = NSING + 1
               IF (SPAT) NUM_SPATIALS = NUM_SPATIALS + 1
            ELSE
C
C====== Multiple
C
               CALL LRL_GET_HARM (ISPOT, H, K, L, H_N, K_N, L_N,
     +                            MULTIPLICITY, MINH, MAXH, INCH,
     +                            ALNOD, DTHRNOD, IERR)
               H2 = (MINH+INCH)*H_N
               K2 = (MINH+INCH)*K_N
               L2 = (MINH+INCH)*L_N

               LAMBDA2 = ALNOD/(MINH+INCH)
               BDMIN2 = DTHRNOD/(MINH)
C
C====== Update histogram
C
               I = MIN(MULTIPLICITY,21)
               MULTI_BINS(I) = MULTI_BINS(I) + 1
               NUM_MULTIPLETS = NUM_MULTIPLETS + 1
               IF (NODFLG) THEN
                  NUM_NODALS = NUM_NODALS + 1
                  IH = H
                  IK = K
                  IL = L
                  NODIDX = NODAL_INDEX(IH,IK,IL)
                  NODAL_BINS(NODIDX) = NODAL_BINS(NODIDX) + 1
               ELSE
                  NODIDX = 0
               END IF
               IF (SPAT) THEN
                  NUM_SPATIALS = NUM_SPATIALS + 1
                  IF (MULTIPLICITY.GT.1) NUM_BOTH = NUM_BOTH + 1
               END IF
C
C====== Now write out .ge1 record for the multiplet
C
               GE1_REC = GE1_REC + 1
               CALL CCPIBY(BUFFER(15),MHKL,-4)
               CALL CCPIBY(BUFFER(17),HKL2,-4)
               IF (MULTIPLICITY.LE.2) THEN
                  WRITE (IUN_GE1,REC=GE1_REC) BUFFER
               ELSE
                  LINK_REC = GE2_REC + 1
                  IREAL = LINK_REC
                  LAMBDA2 = RINT
                  WRITE (IUN_GE1,REC=GE1_REC) BUFFER
C
C====== Write out consecutive records to .ge2 file
C
                  DO 80 I = 2 , MULTIPLICITY
                     H = (MINH+(I-1)*INCH)*H_N
                     K = (MINH+(I-1)*INCH)*K_N
                     L = (MINH+(I-1)*INCH)*L_N
                     LAMBDA1 = ALNOD/(MINH+(I-1)*INCH)
                     GE2_REC = GE2_REC + 1
                     CALL CCPIBY(BUFFER2(1),MHKL,-4)
                     WRITE (IUN_GE2,REC=GE2_REC) BUFFER2
 80               CONTINUE
              END IF
            ENDIF
         END IF

 40   CONTINUE
C
C====== Finally write out header info, number of spots output
C
 90   JPACK = 0
      NFILM = 0
      NUM_SPOTS = GE1_REC - 10
      WRITE (IUN_GE1,REC=9) JPACK , NFILM , NUM_SPOTS
C
C====== Calc cumulative totals
C
      DO 110 I = 1 , 11
         NODAL_BINS(I+1) = NODAL_BINS(I+1) + NODAL_BINS(I)
 110  CONTINUE
C
C====== Now write out header
C
      CALL LDM_MMRAST (IPACK, 1, MM_RAST_X, MM_RAST_Y)

      DO 660 II=1,19
         BUFFER(I) = 0
660   CONTINUE

      DO 666 II=1,6
         CALL LDM_GET_XCENF (IPACK, II, X_CEN_FS(II), IFLAG, ND,
     +                       IERR, ERRSTR)
         CALL LDM_GET_YCENF (IPACK, II, Y_CEN_FS(II), IFLAG, ND,
     +                       IERR, ERRSTR)
         W_FS(II) = 0.0
666   CONTINUE
      CALL LDM_GET_TITLE(COMPOUND, IFLAG)
      WRITE (IUN_GE1,REC=1) COMPOUND(1:76)
      CALL LDM_GET_CRYST (CRYST, IFLAG)
      WRITE (IUN_GE1,REC=2) CRYST , 
     +                      (X_CEN_FS(II)/MM_RAST_X,II=1,6),
     +                      (Y_CEN_FS(JJ)/MM_RAST_Y,JJ=1,6),
     +                      (W_FS(I),I=1,5)

      CALL LDM_GET_SYST (VALSTR, ISYST, IFLAG)
      CALL LDM_GET_LATT (VALSTR, LATT, IFLAG)
      ISYLAT(1) = ISYST
      ISYLAT(2) = LATT
      CALL CCPII2(ISL,ISYLAT,2,.TRUE.,.FALSE.)
      CALL LDM_GET_A (CELL(1), IFLAG, ND)
      CALL LDM_GET_B (CELL(2), IFLAG, ND)
      CALL LDM_GET_C (CELL(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELL(4),IFLAG, ND)
      CALL LDM_GET_BETA (CELL(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELL(6), IFLAG, ND)
      CALL RECCEL(RCELL,CELL,1.0)
      CALL CELLCX(ISYST,CELLC,RCELL,2,2)
      CALL LDM_GET_BEAM (VALSTR, IBEAM, IFLAG)
      CALL LDM_GET_ROT (VALSTR, IROT, IFLAG)
      WRITE (IUN_GE1,REC=3) W_FS(6) , (BUFFER(I),I=1,9) , CELLC , 
     +                      ISL, IBEAM, IROT
C
C====== Initialise integration box parameters
C
      CALL LDM_GET_SPOTD (IPACK, KEYPLATE, DELTA, IFLAG, ND, IERR, 
     +                    ERRSTR)
      BOX(1) = DELTA*MM_RAST_X + 1
C
C====== Force box size to be odd.
C
      IF (MOD(BOX(1),2).EQ.0) BOX(1) = BOX(1) + 1
      BOX(2) = BOX(1)
C
C====== Cut off a corner 1/4 way along.
C
      BOX(3) = BOX(1)/4
      BOX(4) = 1
      BOX(5) = BOX(4)
      CALL LDM_GET_PHIX (IPACK, PHIS(1), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIY (IPACK, PHIS(2), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIZ (IPACK, PHIS(3), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_LMIN (IPACK, WMIN, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_LMAX (IPACK, WMAX, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_DMIN (IPACK, DMIN, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_SPIN (IPACK, SPINDLE, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_RMAX (FILMR, IFLAG, ND)
      WRITE (IUN_GE1,REC=4) 1.0/MM_RAST_X, BOX , 
     +                      (BUFFER(I),I=1,4) , PHIS , 
     +                      1.0 , WMIN , WMAX , DMIN , 
     +                      FILMR , SPINDLE
      DO 777 II = 1, 6
         CALL LDM_GET_XC (IPACK, II, X_C(II), IFLAG, ND,
     +                       IERR, ERRSTR)
         CALL LDM_GET_YC (IPACK, II, Y_C(II), IFLAG, ND,
     +                       IERR, ERRSTR)
         CALL LDM_GET_WC (IPACK, II, W_C(II), IFLAG, ND,
     +                       IERR, ERRSTR)
777   CONTINUE
      WRITE (IUN_GE1,REC=5) (X_C(II),II=1,6), 
     +                      (Y_C(JJ),JJ=1,6),
     +                      (W_C(KK),KK=1,6),
     +                      DELTA
C
C====== Calculate cryst-film ratios
C
      CTOF_RATS(KEYPLATE) = 1.0
      CALL LDM_GET_CTOF (IPACK, KEYPLATE, CTOF_KEYPLATE, IFLAG, ND, 
     +                   IERR, ERRSTR)
      DO 330 I = 1 , 6
         CALL LDM_GET_CTOF (1, I, CTOF, IFLAG, ND, IERR, ERRSTR)
         CTOF_RATS(I) = CTOF/CTOF_KEYPLATE
330   CONTINUE
      CALL LDM_GET_DSTOR (VALSTR, DISTOR_TYPE, IFLAG)
      DO 333 I = 1, 6
         CALL LDM_GET_TWIST (IPACK, I, TWIST(I), IFLAG, ND, IERR, 
     +                       ERRSTR)
         CALL LDM_GET_TILT (IPACK, I, TILT(I), IFLAG, ND, IERR, ERRSTR)
         IF (DISTOR_TYPE.EQ.1) THEN
            CALL LDM_GET_BULGE (IPACK, I, BULGE(I), IFLAG, ND, 
     +                          IERR, ERRSTR)
         ELSE
            BULGE(I) = 0.0
         END IF
333   CONTINUE
      WRITE (IUN_GE1,REC=6) (TWIST(II),II=1,6),
     +                      (TILT(JJ),JJ=1,6),
     +                      (BULGE(KK),KK=1,6),
     +                      CTOF_KEYPLATE
      DO 334 I = 1, 6
         CALL LDM_GET_YSCAL (IPACK, I, Y_SCALE(I), IFLAG, ND, IERR, 
     +                       ERRSTR)
334   CONTINUE
      CALL LDM_GET_EPS (EPSILON, IFLAG, ND)
      WRITE (IUN_GE1,REC=7) (Y_SCALE(II),II=1,6), 
     +                      CTOF_RATS , MININT , EPSILON
      WRITE (IUN_GE1,REC=8) (MULTI_BINS(I),I=1,19)
      CALL LAM_HIST(IPACK, LAMBDA_BINS)
      WRITE (IUN_GE1,REC=9) (MULTI_BINS(I),I=20,21) , 
     +                       (LAMBDA_BINS(I),I=1,16) , 2
      WRITE (IUN_GE1,REC=10) PACK_ID , NUM_SPOTS , NUM_MULTIPLETS , 
     +                       NUM_SPATIALS, NUM_BOTH , NUM_NODALS,
     +                       NODAL_BINS
C
C====== Close .ge1/.ge2 files
C
      CLOSE (UNIT=IUN_GE1)
      CLOSE (UNIT=IUN_GE2)
      RETURN
      END
C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       NODAL_INDEX        **
C******************************
C
C
      INTEGER   FUNCTION NODAL_INDEX (H, K, L)
C
C      INCLUDE 'implicit_none.finc'
C
C Purpose: To calculate the "nodal spot selection index" at which a
C          nodal would be selected.
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER   H
      INTEGER   K
      INTEGER   L
C
C H    (R)   h index
C K    (R)   k index
C L    (R)   l index
C
C====== FUNCTION TYPES:
C
      INTEGER HCF
C
C====== LOCALS:
C
      INTEGER   HCF1, HCF2
      INTEGER   H1, K1, L1
      INTEGER   HI, LO
C
C------------------------------------------------------------------------------
C
      H1 = H
      K1 = K
      L1 = L
      HCF1 = HCF(H1,K1)
      HCF2 = HCF(HCF1,L1)
      H1 = H1/HCF2
      K1 = K1/HCF2
      L1 = L1/HCF2
C
      HI = MAX0(H1,K1,L1)
      LO = MIN0(H1,K1,L1)
      NODAL_INDEX = MAX0(HI,-LO)
C
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LAM_HIST         **
C**************************
C
C
      SUBROUTINE LAM_HIST (IPACK, LAMBDA_BINS)
C
C     INCLUDE 'implicit_none.finc'
C
C Purpose: Get wavelength histogram (16 bins)
C
C Author:  John W. Campbell, February 1992 (based on LAMBDA_HIST from GENLAUE)
C
C Arguments:
C
      INTEGER IPACK
      INTEGER LAMBDA_BINS(16)
C
C IPACK           (R) Pack number
C LAMBDA_BINS(16) (W) Returned histogram
C
C====== Parameter statements
C 
C
C====== GLOBALS:
C
C
C====== FUNCTION TYPES:
C
      
C
C====== LOCALS:
C
      LOGICAL SPAT
      LOGICAL NODFLG
      INTEGER I
      INTEGER J
      INTEGER IFLAG
      INTEGER ND
      INTEGER IERR
      INTEGER NUMSPOTS
      INTEGER JH
      INTEGER JK
      INTEGER JL
      INTEGER MULT
      REAL DTHR2
      REAL XFD
      REAL YFD
      REAL BIN_CUTOFF(16)
      REAL LAMBDA
      REAL WMIN
      REAL WMAX
      CHARACTER*2 ERRSTR

C     LOGICAL SPAT          !Spatial overlap flag
C     LOGICAL NODFLG        !OK nodal flag
C     INTEGER I             !Loop variable
C     INTEGER J             !Loop variable
C     INTEGER IFLAG         !Parameter status flag
C     INTEGER ND            !No. of decimal places flag
C     INTEGER IERR          !Error flag
C     INTEGER NUMSPOTS      !No. of spots in list
C     INTEGER JH            !'h' index
C     INTEGER JK            !'k' index
C     INTEGER JL            !'l' index
C     INTEGER MULT          !Multiplicity
C     REAL DTHR2            !dmin-thresh**2
C     REAL XFD              !Spot 'x' position
C     REAL YFD              !Spot 'y' position
C     REAL BIN_CUTOFF(16)   !Upper Lambda limits for each bin
C     REAL LAMBDA           !Wavelength of current reflection
C     REAL WMIN             !Lambda-min
C     REAL WMAX             !Lambda-max
C     CHARACTER*2 ERRSTR    !Error string (dummy)
C
C-------------------------------------------------------------------------------
C
C
C====== Set up bins
C
      CALL LDM_GET_LMIN (IPACK, WMIN, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_LMAX (IPACK, WMAX, IFLAG, ND, IERR, ERRSTR)
      DO 20 I = 1 , 16
         BIN_CUTOFF(I) = WMIN + I*(WMAX-WMIN)/15.0
         LAMBDA_BINS(I) = 0
 20   CONTINUE
      BIN_CUTOFF(16) = 100.0
C
C====== Work out histogram
C
      CALL LRL_NUMSPOTS (NUMSPOTS)
      DO 50 J = 1 , NUMSPOTS
         CALL LRL_GET (J, JH, JK, JL, XFD, YFD, LAMBDA, DTHR2, MULT,
     +                       NODFLG, SPAT, IERR)
         DO 30 I = 1 , 16
            IF (LAMBDA.LT.BIN_CUTOFF(I)) GO TO 40
 30      CONTINUE
         I = 16
 40      LAMBDA_BINS(I) = LAMBDA_BINS(I) + 1
 50   CONTINUE
      RETURN
      END
C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       HCF           **
C*************************
C
C
      INTEGER    FUNCTION  HCF (II,JJ)
C
C      INCLUDE 'implicit_none.finc'
C
C Purpose: To calculate h.c.f. of 2 integers
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER    II
      INTEGER    JJ
C 
C II    (R)   First integer
C JJ    (R)   Second integer
C
C====== LOCALS:
C
      INTEGER    I, J, ISAV, K
C
C------------------------------------------------------------------------------
C
      I = II
      J = JJ
      IF (I.EQ.J) THEN
         HCF = I
         RETURN
      END IF

      IF (I.LT.J) THEN
         ISAV = I
         I = J
         J = ISAV
      END IF

      IF (J.EQ.0) THEN
         HCF = I
         RETURN
      END IF

 20   K = MOD(I,J)

      IF (K.EQ.0) THEN
         HCF = J
         RETURN
      END IF

      I = J
      J = K
      GO TO 20

      END
