C
C==== CAD4(Version 5) Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 5th February 1997
C
      subroutine refrd(itype)
C**** =======================================================================
C     Routine to interpret reflection data from a CAD4.
C     All data is output by calling REFOUT
C**** =======================================================================
      integer itype,i
      real pmt
      logical oka
C
      include 'common.for'
C---- Reflections with Goniometer status C have no intensity measurement
      if (( code6(2:2) .ne. 'C') .and. ( code6(2:2) .ne. 'X')) then
C
C---- In the DIFRAC module for the CAD4 version 5.0 diffractometer, the
C---- UB, UBINV and UBDET variables contain the values reset for a negative
C---- determinant. So for the psi calculation the transformed hkls are used.
C---- This is possible because the CAD4 keeps the setting angles with each
C---- reflection.
C
C---- Transform h k l for negative UB
C
        if (trans) then
          ih = -ih
          ik = -ik
          il = -il
        endif
C---- Calculate crystal based psi value
        call psic
C---- Set standard flag
        ptstd = -999
        if ( code6(1:1) .eq. 'N') then
          ptstd = 0
          stdseq = 0
        elseif (code6(1:1) .eq. 'I') then
          stdseq = stdseq + 1
          if ((ih.eq.href(stdseq)).and.(ik.eq.kref(stdseq)).and.
     1        (il.eq.lref(stdseq))) then
            ptstd = stdseq
          else
            do 100 i = 1, numref
              if ((ih.eq.href(i)).and.(ik.eq.kref(i)).and.
     1          (il.eq.lref(i))) ptstd = i
 100        continue
          endif
        endif
        if(ptstd.eq.-999) write(*,
     1   '('' Unlisted standard '',3i5)') ih,ik,il
C---- Totalize the standards
        if(ptstd.gt.0) then
          avint(ptstd) = avint(ptstd) + peak
          sdint(ptstd) = sdint(ptstd) + peak*peak
          count(ptstd) = count(ptstd) + 1.0
        endif
C---- Check up that time increases
        if ( (time + bastim) .le. prvtim)
     1     bastim = prvtim - time + 1.0
        time   = time + bastim
        prvtim = time
C---- Set attenuator flag and speed
        if ( spe . ge. 0.0) then
          ptfilt = 0
        else
          ptfilt = 1
        endif
        spe = 16.48/abs(spe)
C
C---- Horizontal detector(aperture) slit opening
C
C---- Get sin theta
        sint  = sin(dar*abs(theta))
        sintl = sint/wavel
        dethz = apta + aptb * sint / sqrt(1.0 - sint*sint)
        if (dethz .lt. apmin) dethz = apmin
        if (dethz .gt. apmax) dethz = apmax
        dethz = rad * dethz / radius
C
C---- Scan width
C
        swidth = 0.66666667*abs(swidth)
C---- Calculation of transformation coefficients
        if (statbk) then
C---- Stationary background measuring technique
          cnp = 0.24*timp/swidth
          cnb = 0.24*timp
          cnl = 0.24*timl
          cnr = 0.24*timr
          backsn = 0
        else
C---- Scanned background measuring technique
C       Pmt - Peak-with-background measuring time
          pmt = swidth/spe
C       Bgt - left and right background counting time.
          bgt = 0.25 * pmt
          cnp = pmt/swidth
          cnb = pmt
          cnl = bgt
          cnr = bgt
          backsn = 1
        endif
C---- Count the reflections
        nrefl = nrefl + 1
C---- Keep min/max values of h, k, l
        hmin = min(hmin, ih)
        kmin = min(kmin, ik)
        lmin = min(lmin, il)
        hmax = max(hmax, ih)
        kmax = max(kmax, ik)
        lmax = max(lmax, il)
C---- Keep min/max values of sintl
        stlmin = min(sintl, stlmin)
        stlmax = max(sintl, stlmax)
C---- Keep on eye on the filter codes that have been used
        oka = .false.
        do 2010 i = 1, nfilt
          if (ptfilt .eq. nfltcd(i)) oka = .true.
 2010   continue
        if (.not. oka) then
          nfilt = nfilt + 1
          nfltcd(nfilt) = ptfilt
        endif
C
        call refout
      endif
      return
      end
C
      subroutine dtflrd(jj)
C
C===== Subroutine to read a CAD4 diffractometer output file
C
      integer i,j,jj, mrefl, trim
      integer mdatln, mdlim
      real sinsqa
C
      include 'common.for'
C
C---- Initialisations
C
      reoren = .false.
      mrefl = 0
      parclp = 0
      mdatln = 0
      mdlim  = 3
C
C---- Get the data file name and open it
C
      if (jj .eq. 0) then
        ptscl = 1
        open( unit = datafl, file = compid(1:trim(compid))//datxtn,
     1  status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     2  ,readonly)
CVMS 2  ,readonly)
CUNX 2  )
        write(*,'(/'' Reading file: '',A)')
     1   compid(1:trim(compid))//datxtn
      else
        ptscl = jj
        i = jj/10
        j = mod(jj, 10)
        open( unit = datafl, 
     1  file = compid(1:trim(compid))//datxto(1:trim(datxto))
     2    //digit(i)//digit(j), status = 'OLD', access = 'SEQUENTIAL', 
     3    form = 'FORMATTED'
     4  ,readonly)
CVMS 4  ,readonly)
CUNX 4  )
        write(*,'(/'' Reading file: '',A)') compid(1:trim(compid))
     1    //datxto(1:trim(datxto))//digit(i)//digit(j)
      endif
C****
C     Interpretation of the .DAT file
C****
 1000 continue
      line = ' '
      read(datafl, '(a64)', end = 2000) line
C---- Shift the buffer if necessary
      if (slidc4.gt.0) then
        do 1010 i = 80-slidc4,1,-1
          line(i+slidc4:i+slidc4)=line(i:i)
 1010   continue
      else if (slidc4.lt.0) then
        do 1015 i = 1-slidc4,80
          line(i+slidc4:i+slidc4) = line(i:i)
 1015   continue
      endif
C---- Find the line type
      if ( line(2:6) .eq. '     ') then
C---- Do nothing - we need to ignore entirely blank lines
C
C---- Type 0 line. Get the check reflection indices
C
      elseif (line(2:4) .eq. '   ') then
        if (((line(16:16).eq.'H').or.(line(16:16).eq.'N')).and.
     1     ((line(20:20).eq.'I').or.(line(20:20).eq.'F'))) then
          numref = numref + 1
          read(line,'(6X,3I3)')
     1      href(numref), kref(numref), lref(numref)
        endif
C
C---- Type 21 line.
C
      elseif (line(2:4) .eq. ' 21') then
        numref = 0
        stdseq = 0
        parclp = 1
        read(line,'(46x,f4.0,3f4.1)')
     1    radius, slit, apmin, apmax
C
C---- Type 22 line.
C
      elseif (line(2:4) .eq. ' 22') then
        parclp = parclp + 2
        read(line,'(26x,2f5.2,f2.0,6x,f6.3)')
     1    apta, aptb, scnrat, csigma
C
C---- Type 23 line.
C
      elseif (line(2:4) .eq. ' 23') then
        parclp = parclp + 4
C
C---- Type 24 line.
C
      elseif (line(2:4) .eq. ' 24') then
        parclp = parclp + 8
C
C---- Type 25 line.
C
      elseif (line(2:4) .eq. ' 25') then
        parclp = parclp + 16
C
C---- Type 26 line.
C
      elseif (line(2:4) .eq. ' 26') then
        read(line,'(4X,f7.0,24x,2f10.7)')
     1    sinter, sinsqa, cosalp
        sinalp = sqrt(sinsqa)
        parclp = parclp + 32
C
C---- Type 27 line.
C
      elseif (line(2:4) .eq. ' 27') then
        parclp = parclp + 64
C
C---- Type 31 line
C
      elseif (line(2:4) .eq. ' 31') then
        parclp = parclp + 128
        read(line,'(6x,3f9.6,2x,3f9.6)')
     1   ub(1,1),ub(1,2),ub(1,3),ub(2,1),ub(2,2),ub(2,3)
C
C---- Type 32 line
C
      elseif (line(2:4) .eq. ' 32') then
        parclp = parclp + 256
        read(line,'(6x,3f9.6,2x,3f9.6)')
     1   ub(3,1),ub(3,2),ub(3,3),flam1,flam2,attn
C
C---- Only allow the reflection lines if control data present
C
      elseif (parclp .eq. 1023) then
C
C---- Type 1 line
C
        if (line(2:4) .eq. '  1') then
          read(line,'(4x,i6,3i5,1x,a6,f7.2,f4.0,f6.0,f7.0,f6.0)')
     1      mrefl, ih, ik, il, code6, psi, spe, backl, peak, backr
          mdatln = 1
C     Allow for possible line type 19 output
          if (code6(1:1) .eq. 'S')  then
            mdlim = 7
            statbk = .true.
          else
            mdlim = mdlmz
            statbk = .false.
          endif
C
C---- Type 2 line
C
        elseif (line(2:4) .eq. '  2') then
          read (line,'(4x,i6,f8.3,3f9.3,f7.3,f7.0,i2,i1)')
     1      i, theta, phi, omega, kappa, swidth, time, nfridl, apc
          if (i.eq.mrefl) mdatln = mdatln + 2
C
C---- Type 19 line
C
        elseif (line(2:4) .eq. ' 19') then
          read (line,'(4x,i6,33x,f6.0,f7.0,f6.0)')
     1      i, timl, timp, timr
          if ( i .eq. mrefl) mdatln = mdatln + 4
        endif
      endif
C
      if (parclp .eq. 511) then
        parclp = 1023
        call initlz
      endif
      if (parclp .eq. 1407) then
        reoren = .true.
        parclp = 1023
        call initlz
      endif
C
C---- Process reflection if all lines have been found
C
      if (mdatln .eq. mdlim) then
        mdatln = 0
        call refrd(0)
      endif
C
C---- Bottom of line reading loop
C
      go to 1000
C
C---- Finish off this data file
C
 2000 continue
      close (unit = datafl, status = 'KEEP')
      return
      end
C
      function junpak(val)
      integer*2 i,junpak, val, tmp, zmp, xmp
      tmp=val
      zmp=0
      do 147 i=1,4
      xmp=tmp-8*(tmp/8)
      tmp=tmp/8
  147 zmp=(zmp/10)+1000*xmp
      junpak=zmp
      return
      end
C
      subroutine choos
C
C----- Set up fixed characteristics for CAD4 diffractometer
C
      include 'common.for'
C
C==== Set Logical Unit Numbers for DIFRAC
C
C     The program uses extensively the standard input and output files
C     as accessed by the FORTRAN instructions read(*,.. and write(*,.. .
C
C     datafl is the unit used to read the diffractometer data files
C     scfsfl is the unit used to write the standardized output file
C     parmfl is the unit used to read the diffractometer parameter files
C     addrfl is the unit used to read the address list file 
C     qaofl  is the unit used to write/read an existing saved answer file
C     qanfl  is the unit used to read/write a new saved answer file
C
      datafl =  9
      scfsfl = 10
      parmfl = 11
      addrfl = 12
      qaofl  = 13
      qanfl  = 14
C
C---- Tell them its an Enraf-Nonius diffractometer
      difnam = 'E'
      device = 'Enraf-Nonius CAD4 4circle diffractometer'
      softwr = 'Enraf-Nonius CAD4 version 5.0           '
      cellpg = '?                                       '
      dettyp = 'scintillation counter                   '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = '[flack.cristal]address.lst                         '
CNORMAL addnam = 'address.lst                                        '
        addnam = '[flack.cristal]address.lst                         '
C
C---- File extension for data file name
      datxtn = '.dat'
      datxto = '.d'
C
C---- Set dihedral angle
      fmnhed = 90.0
C
C---- Default wavelength pointer for cell dimension determination, alpha 1
      cdefwv = 1
C
C---- Set the buffer shift for .dat files of CAD4. Enraf-Nonius has
C---- chosen to have a very funny ASCII file structure for their *.dat
C---- files. A consequence of this is that when the files are shifted from
C---- machine to machine the character get shifted one place to the right or
C---- left. The variable slidc4 sets the shift to get back to the standard
C---- position usually obtained with a DIFRAC running on the machine which
C---- produced the file.
C
      slidc4 = 0
      return
      end
C
      subroutine initlz
C
C**** Subroutine to interpret a complete 21-26 information in memory
C     for the CAD4.
C****
      integer i,j,k
      real ubdet
C
      include 'common.for'
C
C---- Get UB**-1
C
      ubinv(1,1) = +ub(2,2)*ub(3,3) - ub(3,2)*ub(2,3)
      ubinv(2,1) = -ub(2,1)*ub(3,3) + ub(3,1)*ub(2,3)
      ubinv(3,1) = +ub(2,1)*ub(3,2) - ub(3,1)*ub(2,2)
      ubdet=ub(1,1)*ubinv(1,1)+ub(1,2)*ubinv(2,1)+ub(1,3)*ubinv(3,1)
      ubinv(1,1) = ubinv(1,1)/ubdet
      ubinv(2,1) = ubinv(2,1)/ubdet
      ubinv(3,1) = ubinv(3,1)/ubdet
      ubinv(1,2) = (-ub(1,2)*ub(3,3) + ub(3,2)*ub(1,3))/ubdet
      ubinv(2,2) = (+ub(1,1)*ub(3,3) - ub(3,1)*ub(1,3))/ubdet
      ubinv(3,2) = (-ub(1,1)*ub(3,2) + ub(3,1)*ub(1,2))/ubdet
      ubinv(1,3) = (+ub(1,2)*ub(2,3) - ub(2,2)*ub(1,3))/ubdet
      ubinv(2,3) = (-ub(1,1)*ub(2,3) + ub(2,1)*ub(1,3))/ubdet
      ubinv(3,3) = (+ub(1,1)*ub(2,2) - ub(2,1)*ub(1,2))/ubdet
C
C**** Calculate (approximate) direct space cell dimensions from UB matrix
C
C---- Get real metric in GM from ubinv
      do 400 i = 1,3
      do 400 k = 1,3
        gm(i,k) = 0.0
      do 400 j = 1,3
        gm(i,k) = gm(i,k) + ubinv(i,j)*ubinv(k,j)
  400 continue
C---- Get cell dimensions from ub
      do 401 i = 1, 3
       cellub(i) = sqrt(gm(i,i))
 401  continue
      do 402 i = 1,3
       j = 1 + i/3
       k = i -j + 2
       cellub(7-i) = rad*acos(gm(j,k)/(cellub(j)*cellub(k)))
 402  continue
       volub=abs(1.0/ubdet)
C---- Make angles 90 and 120 and cell lengths equal
      do 404 i = 4, 6
        if(abs(cellub(i)-90.) .lt. 0.4) cellub(i) = 90.
 404  continue
      if(abs(cellub(6)-120.) .lt. 0.4) cellub(6) = 120.
      if(abs(cellub(1)-cellub(2)).lt.0.01) then
        if(abs(cellub(1)-cellub(3)).lt.0.01) then
          cellub(1) = (cellub(1)+cellub(2)+cellub(3))/3.0
          cellub(2) = cellub(1)
          cellub(3) = cellub(1)
        else
          cellub(1) = 0.5*(cellub(1)+cellub(2))
          cellub(2) = cellub(1)
        endif
      endif
C
C**** Check for negative Ub mark to transform UB and hkl if necessary
C
      if (ubdet .gt. 0.0) then
        trans = .false.
      else
        trans = .true.
        write(*,'('' ***** Negative ub matrix.'',
     1  '' Transforming UB and h k ls*****'')')
        ubdet = - ubdet
        do 100 i = 1, 3
        do 100 j = 1, 3
          ub(i,j)    = - ub(i,j)
          ubinv(i,j) = - ubinv(i,j)
 100    continue
        do 105 i = 1, numref
          href(i) = - href(i)
          kref(i) = - kref(i)
          lref(i) = - lref(i)
 105    continue
      endif
C
C---- Do not do the rest of the calculations for a reorientation only
C
      if (reoren) then
C---- Reset reoren so that INITLZ will not be called again
        reoren = .false.
      else
C
C---- Convert the scan type parameter
C
        scnrat = scnrat/3.0
        mstnn = 3
        if(abs(scnrat).lt. 0.1)       mstnn = 1
        if(abs(scnrat - 2.0).lt. 0.1) mstnn = 2
C---- Wavelength
        weigh1 = 0.6666667
        weigh2 = 0.3333333
        weigh3 = 0.0000000
        flam3  = 0.0
        wavel = weigh1*flam1 + weigh2*flam2
C
C---- Set up the wavelength for output
C
        do 201 i = 1, navtab
          if (abs(wavtab(1,i)-wavel).le.0.005) then
            flam1  = wavtab(2,i)
            flam2  = wavtab(3,i)
            weigh1 = wavtab(4,i)
            weigh2 = wavtab(5,i)
            tubtar = wavnam(1,i)
            wl0nam = wavnam(2,i)
            wl1nam = wavnam(3,i)
            wl2nam = wavnam(4,i)
          endif
 201    continue
        cwave = flam1
        clrnam = wl1nam
C---- Vertical detector opening
        detvt = rad*slit/radius
C---- Other filter variables
        if (attn .gt. 1.0) then
          nfilt= 2
          nfltcd(2)=1
          attfac(2)= attn
          ettfac(2)= 0.01*attn
        endif
C
C---- Set mdlmz: value 3 means expect lines type 1 & 2 only,
C                value 7 means expect lines type 1, 2 & 19.
        mdlmz = 3
C---- Learnt profile output
        if (csigma .lt. 0.0) mdlmz = 7
C
      endif
C
      return
      end
C
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Enraf-Nonius CAD4 machine.
C
C     The equations used for the calculation of the crystal-based azimuth
C     are given in "On the definition and Practical Use of Crystal-based
C     Azimuthal Angles" by D.Schwarzenbach and H.D.Flack (1989)
C     J. Appl. Cryst. 22, 601-605; D.Schwarzenbach and H.D.Flack (1992)
C     J. Appl. Cryst. 25, 69. hereafter DSHDF.
C     The definition of the reference vector u is given in equation (3) of
C     DSHDF and equation (9) gives the psi-crystal.
C
C     The CAD-4 corresponds to case (a) of section 3) "Dealing with
C     Manufacturers" design of DSHDF since 2-theta, omega, phi and kappa
C     are found on the diffractometer output file as well as UB each time
C     it changes. The CAD4 is thus a DSHDF type I diffractometer and hence
C     W = 1 according to DSHDF equation (11). (W and 1 being 3x3 matrices).
C
C     The CAD4 uses the kappa geometry described in detail in the CAD4
C     User's Manual, Chapter II, "CAD4 Geometry", hereafter CAD4MII.
C     At page 6 of CAD4MII one finds the matrix expression for the setting
C     of the goniometer to angles omega, kappa and phi as:
C         Z(omega) * Y(-alpha) * Z(Kappa) * Y(+alpha) * Z(phi).
C     Thus for the kappa geometry, equation (9) of DSHDF becomes:
C
C    T        T     -1
C Psi  = s * u  * UB   * Z(-phi) * Y(-alpha) * Z(-kappa) * Y(+alpha) * Z(-omega)
C
C               -2     T     -1     -1T
C        where s   =  u  * UB   * UB    * u    and
C
C                (  1          0           0          )
C     X(alph) =  (  0          cos(alph)   sin(alph)  )
C                (  0         -sin(alph)   cos(alph)  )
C
C                (  cos(beta)  0          -sin(beta)
C     Y(beta) =  (  0          1           0          )
C                (  sin(beta)  0           cos(beta)  )
C
C                (  cos(gamm)  sin(gamm)   0          )
C     Z(gamm) =  ( -sin(gamm)  cos(gamm)   0          )
C                (  0          0           1          )
C
C
      integer  i,j
      real r(3), s(3), sinq, cosq
      include 'common.for'
C
C---- Form the reference vector u: equation (3) of DSHDF
C
      if ((ih.eq.ik).and.(ik.eq.il)) then
        r(1) = real(ih)
        r(2) = real(-ih)
        r(3) = 0.0
      else
        r(1) = real(ik-il)
        r(2) = real(il-ih)
        r(3) = real(ih-ik)
      endif
C---- u(transpose)*ubinv
      do 110 i = 1 ,3
        s(i) = 0.0
        do 105 j = 1, 3
          s(i) = s(i) + r(j)*ubinv(j,i)
 105    continue
 110  continue
C---- *z(-phi)
      sinq = sin(dar*phi)
      cosq = cos(dar*phi)
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *y(-alph)
      s(1) =  r(1)*cosalp - r(3)*sinalp
      s(2) =  r(2)
      s(3) =  r(1)*sinalp + r(3)*cosalp
C---- *z(-kappa)
      sinq = sin(dar*kappa)
      cosq = cos(dar*kappa)
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *y(+alph)
      s(1) =  r(1)*cosalp + r(3)*sinalp
      s(2) =  r(2)
      s(3) = -r(1)*sinalp + r(3)*cosalp
C---- *z(-omega)
      sinq = sin(dar*(omega-theta))
      cosq = cos(dar*(omega-theta))
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *W (DSHDF type I, W is a unit matrix, nothing to be done here
C---- Measurements at negative theta correspond to 
C         (1    0    0)
C       W (0   -1    0)
C         (0    0   -1)
C     so we must take care of that
      if (theta .lt. 0.0) then
        r(2) = -r(2)
        r(3) = -r(3)
      endif
C---- Crystal based azimuth
      psicry = rad*atan2(r(1),r(3))
C
C==== r(2) = r(2)/ SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
C==== WRITE(12,'(/'' psi(2) '', e15.6)') r(2)
C==== Look at them to see if they agree
C==== WRITE(12,'(1x,3i5,6f10.2)')
C====1 ih,ik,il, omega, 2*theta, kappa, phi, psi, psicry
C
      return
      end
C
      subroutine prflrd
C
C===== Subroutine to read a CAD4 .PAR and .CRY files
C
      include 'common.for'
      logical     zresnt, print
      integer trim
C
C---- Read the GONCA0 file of a CAD4 diffractometer (CAD4 version 5.0)
C
C==== Variable declarations for GONCAn.PAR
C
      real        thpos, thneg, zudius, zlit, con1, con2, con3
      real        zapmax, zapmin, pdist
      real        fshv, fsll, fswd, ftime, baud, fspa, fsmt, fsmp, fsmo
      real        fsmk, fsmd, fsem
      integer*2   junpak
      integer*2   maxvar, minvar, upperh, lowerh, negsl, possl
      integer*2   vslit,  hslit,  hole
      integer*2   shv,  sll,  swd, sdth, sdtl, sbrl, sbrh, scsl, scsh
      integer*2   spa,  smt,  smp,  zmo,  smk,  smd, sflags, sem
      character*64  cry,dat,hkl
C
C==== Variable declarations  <compid>.CRY
C
      integer     nrec, i, l, m
      integer*2   iscant, npipre, itmax, idumpf, nbalf, mfridl, neqfl
      integer*2   nrpsi, indzz(3,3), nfms, nabs, labs(3,6), ichkt
      integer*2   norchk, nrhkl, npid, nics, nors, npsis, npifl
      integer*2   nnrefl
      integer*2   irsind,irsth,irsang,irscor,irscin,irsscn
      integer*2   listnr,irnpi
      real        r(3,3), det, d(3,3), flam(2), zattn, params(7)
      real        intca, dvect(3), dhkl(3), dints(6)
      real        fval(7,9)
      real        rhhh, rkkk, rlll, rtheta, rphik, romk, rkappa
      real        rscang, rscint, dthmin, dthmax, doma, domb, sigpre
      real        sigma, psifl, psist, delpsi, seq, sigmi, fading
      real        dang, xrtime, dtheta, domk, datpsi, apt
      character*36 dtext
      character*6 cption(8), chkl(6), cpsifl, cseq, cnfms
      real        rption(8)
      integer     rhkl(6)
      data cption(1) /'BISECT'/
      data cption(2) /'AZIMUT'/
      data cption(3) /'FLAT  '/
      data cption(4) /'NEEDLE'/
      data cption(5) /'VECHOR'/
      data cption(6) /'ZIGZAG'/
      data cption(7) /'SEPHKL'/
      data cption(8) /'FILE  '/
      data chkl(1)   /'HKL   '/
      data chkl(2)   /'HLK   '/
      data chkl(3)   /'KHL   '/
      data chkl(4)   /'KLH   '/
      data chkl(5)   /'LHK   '/
      data chkl(6)   /'LKH   '/
CVAX
CVAX---- The R format is only allowed in DATA statements, is very specific
CVAX     to VMS/VAX FORTRAN, packs three characters into two bytes and is only
CVAX---- because the Enraf-Nonius software uses it.  Why? Why? Why? Why?
CVAX
CVAX      data rption(1) /6rbisect/
CVAX      data rption(2) /6razimut/
CVAX      data rption(3) /6rflat  /
CVAX      data rption(4) /6rneedle/
CVAX      data rption(5) /6rvechor/
CVAX      data rption(6) /6rzigzag/
CVAX      data rption(7) /6rsephkl/
CVAX      data rption(8) /6rfile  /
CVAX      data rhkl(1)   /3rhkl   /
CVAX      data rhkl(2)   /3rhlk   /
CVAX      data rhkl(3)   /3rkhl   /
CVAX      data rhkl(4)   /3rklh   /
CVAX      data rhkl(5)   /3rlhk   /
CVAX      data rhkl(6)   /3rlkh   /
C
C     But the following from Michael Estermann works on VMS/ALPHA and lets hope
C     it also works for VMS/VAX
C
      DATA rption(1) /  3.8691810E-31 /
      DATA rption(2) /  2.6160208E-33 /
      DATA rption(3) /  7.8231963E-16 /
      DATA rption(4) /  2.2551539E+14 /
      DATA rption(5) / -2.4459881E-33 /
      DATA rption(6) / -6.4953328E-18 /
      DATA rption(7) /  1.5435138E+33 /
      DATA rption(8) /  4.2348429E-16 /
      DATA   rhkl(1) /       13252    /
      DATA   rhkl(2) /       13291    /
      DATA   rhkl(3) /       17932    /
      DATA   rhkl(4) /       18088    /
      DATA   rhkl(5) /       19531    /
      DATA   rhkl(6) /       19648    /
C
C---- Inquire GONCA0.PAR and open it if present
C
      inquire(file='gonca0.par',exist=zresnt)
      if(.not.zresnt) then
        write(*,'(''   File: GONCA0.PAR not available'')')
      else
        write(*,'(/'' Reading file: GONCA0.PAR'')')
        open( unit=parmfl,file='gonca0.par',status='OLD',
     1   recl=16, access='DIRECT', err=60, readonly)
C
        read(unit=parmfl,rec=1) cry
        read(unit=parmfl,rec=2) dat
        read(unit=parmfl,rec=3) hkl
        read(unit=parmfl,rec=4) thpos,thneg,zudius,zlit,con1,con2,con3,
     1   zapmax,zapmin,pdist,maxvar,minvar,upperh,lowerh,negsl,possl,
     2   vslit,hslit,hole
C---- Convert the packed values!?
        maxvar=junpak(maxvar)
        minvar=junpak(minvar)
        upperh=junpak(upperh)
        lowerh=junpak(lowerh)
        negsl =junpak(negsl)
        possl =junpak(possl)
        vslit =junpak(vslit)
        hslit =junpak(hslit)
        hole  =junpak(hole)
C
        read(unit=parmfl,rec=5)
        read(unit=parmfl,rec=6)
        read(unit=parmfl,rec=7)
        read(unit=parmfl,rec=8) shv,sll,swd,sdth,sdtl,sbrl,sbrh,scsl,
     1   scsh,spa,smt,smp,zmo,smk,smd,sflags,sem
C==== Convert values
        fshv=real(255-iand(shv,255))*3.0+255.0
        fsll=real(255-iand(sll,255))*5.0
        fswd=real(255-iand(swd,255))*5.0
        ftime=real((dble(float(sdth))*(2.0**15)+dble(float(sdtl)))
     1   *(1000000.0/400.0)/(2.0**32))
        baud=115200.0/(real(iand(sbrh,255))*2**8+float(iand(sbrl,255)))
        fspa=real(spa)
        fsmt=real(smt)
        fsmp=real(smp)
        fsmo=real(zmo)
        fsmk=real(smk)
        fsmd=real(smd)
        fsem=(iand(sem,15))*10.0/3.0+10.0
C
C---- Use the values if possible
C
C---- Pick up the dead time
        if (ftime .gt. 0.0) then
          fdeadt = ftime
          efdead = 0.0
          ask(3) = .false.
        endif
C
C---- Write out the values if needed
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
        write(*,'('' Do you wish to display the contents of '',
     1   ''this file (Y,N)? '')')
        call fread(000040,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'Y') then
          write(*,'(/'' File names:'')')
          write(*,'(''    CRY Name '',A64)') cry
          write(*,'(''    DAT Name '',A64)') dat
          write(*,'(''    HKL Name '',A64)') hkl
          write(*,'(/'' Goncon parameters:'')')
          write(*,'(4X,F11.7,'' Goniometer constant CON1'')') con1
          write(*,'(4X,F11.7,'' Goniometer constant CON2'')') con2
          write(*,'(4X,F11.7,'' Goniometer constant CON3'')') con3
          write(*,'(4X,F11.2,'' degrees Max theta'')') thpos
          write(*,'(4X,F11.2,'' degrees Min theta'')') thneg
          write(*,'(4X,F11.2,'' mm. Centre gonio to screen'')') pdist
          write(*,'(4X,F11.2,'' mm. Centre gonio to aperture'')') zudius
          write(*,'(4x,f11.2,
     1      '' mm. Vertical height insertable slit'')') zlit
          write(*,'(4x,f11.2,
     1      '' mm. Max horizontal size of variable slit'')') zapmax
          write(*,'(4x,f11.2,
     1      '' mm. Min horizontal size of variable slit'')') zapmin
          write(*,'(/'' Encoder values for:'')')
          write(*,'(4X,I11,'' Max variable aperture '')') maxvar
          write(*,'(4X,I11,'' Min variable aperture'')') minvar
          write(*,'(4X,I11,'' Upper half-aperture'')') upperh
          write(*,'(4X,I11,'' Lower half-aperture'')') lowerh
          write(*,'(4X,I11,'' Negative slanted slit'')') negsl
          write(*,'(4X,I11,'' Positively slanted slit'')') possl
          write(*,'(4X,I11,'' Vertical slit'')') vslit
          write(*,'(4X,I11,'' Horizonal slit'')') hslit
          write(*,'(4X,I11,'' 9mm hole'')') hole
          write(*,'(/'' Caspar values:'')')
          write(*,'(4X,F11.1,'' volts Detector setting'')') fshv
          write(*,'(4X,F11.1,'' Detector lower level'')') fsll
          write(*,'(4X,F11.1,'' Detector window setting'')') fswd
          write(*,'(4X,F11.4,'' microsecs. Detector dead time'')') ftime
          write(*,'(4X,F11.1,'' baud for interface terminal'')') baud
          write(*,'(4X,F11.1,'' positionning accuracy'')') fspa
          write(*,'(4X,F11.1,'' Emission level setting'')') fsem
          write(*,'(4X,I11,''  Flags'')') sflags
          write(*,'(4X,I11,''  Clock speed low byte'')') scsl
          write(*,'(4X,I11,''  Clock speed high byte'')') scsh
          write(*,'(/'' Motor gain values:'')')
          write(*,'(4X,F11.1,'' Theta motor'')') fsmt
          write(*,'(4X,F11.1,'' Phi motor'')') fsmp
          write(*,'(4X,F11.1,'' Omega motor'')') fsmo
          write(*,'(4X,F11.1,'' Krappa motor'')') fsmk
          write(*,'(4X,F11.1,'' Dial motor'')') fsmd
        endif
C---- Finish off
 60   close(unit=parmfl,status='KEEP')
      endif
C
C==== Read the <Compid>.CRY file of a CAD4 diffractometer (CAD4 version 5.0)
C
C---- Inquire <Compid>.CRY and open it if present
C
      inquire(file=compid(1:trim(compid))//'.cry',
     1 exist=zresnt)
      if(.not.zresnt) then
        write(*,'(''   File: '',A,''.CRY not available'')')
     1   compid(1:trim(compid))
      else
        write(*,'(/'' Reading file: '',A,''.CRY'')')
     1   compid(1:trim(compid))
        open( unit=parmfl,file=compid(1:trim(compid))//
     1   '.CRY', status='OLD', recl=128, access='DIRECT', readonly,
     2   err = 70)
        write(*,'('' Do you wish to display the contents of '',
     1   ''this file (Y,N)? '')')
        call fread(000060,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        print = .false.
        if (inchar .eq. 'Y') print = .true.
C
C---- Treat 25 reflection records
C
        do 2000 nrec = 1, 25
          read(unit=parmfl,rec=nrec)
     1     irsind,irsth,irsang,irscor,irscin,irsscn,listnr,irnpi,
     2     rhhh,rkkk,rlll,rtheta,rphik,romk,rkappa,rscang,rscint,
     3     fval
C
C---- Use the values if possible
C
C---- Pick up the two theta angles for "Good" angle reflections
          if((irsth.eq.84).or.(irsang.eq.83)) then
            numang = numang + 1
            hang(numang) = nint(rhhh)
            kang(numang) = nint(rkkk)
            lang(numang) = nint(rlll)
            twoang(numang) = 2.0*rtheta
            wang(numang) = cdefwv
            ctthn = amin1(twoang(numang),ctthn)
            ctthx = amax1(twoang(numang),ctthx)
          endif
          nmcell = numang
          if (print) then
            write(*,'(/'' List reflection no:'',I3)') nrec
            write(*,'(I11,'' List sequence no.'')') listnr
            write(*,'(10X,A1,'' Indices status'')') char(irsind)
            write(*,'(I11,'' h'')') nint(rhhh)
            write(*,'(I11,'' k'')') nint(rkkk)
            write(*,'(I11,'' l'')') nint(rlll)
            write(*,'(10X,A1,'' Theta angle status'')') char(irsth)
            write(*,'(10X,A1,'' Angle status'')') char(irsang)
            write(*,'(F11.3,'' degrees. Theta'')') rtheta
            write(*,'(F11.2,'' degrees. Phi'')') rphik
            write(*,'(F11.2,'' degrees. Omega'')') romk
            write(*,'(F11.2,'' degrees. Kappa'')') rkappa
            write(*,'(10x,A1,'' Orientation control indicator'')')
     1       char(irscor)
            write(*,'(10X,A1,'' Intensity control indicator'')')
     1       char(irscin)
            write(*,'(10X,A1,'' Scan information status'')')
     1       char(irsscn)
            write(*,'(I11,'' Attenuator/scan speed param'')') irnpi
            write(*,'(F11.2,'' degrees. Scan angle'')') rscang
            write(*,'(F11.1,'' Net Intensity'')') rscint
            write(*,'(/''  The assignment of the values in this '',
     1       '' dump of FVAL is dubious!'')')
            write(*,'(''     theta      phi    omega      chi'',
     1       ''      Npi    Scang    Scint'')')
            write(*,'(9(1X,7F9.2,I2/))') ((fval(l,m),l=1,7),m,m=1,9)
          endif
 2000   continue
C---- Treat Records 26 27 and 28
        read(unit=parmfl,rec=26)
     1   iscant, npipre, itmax, idumpf, nbalf, mfridl, neqfl, nrpsi,
     2   indzz, nfms, nabs, labs, ichkt, norchk, nrhkl, npid, nics,
     3   nors, npsis, npifl, nnrefl,
     4   dtext, dthmin, dthmax, doma, domb, apta, aptb, sigpre, sigma,
     5   psifl, psist, delpsi, dvect, seq, sigmi, fading, dang,
     6   intca, xrtime, dtheta, domk, dhkl, datpsi, dints
        read(unit=parmfl,rec=27) r, det, d, flam, zattn, params
        read(unit=parmfl,rec=28)
C
C---- Convert from Radix-50 to Character values
C
        cnfms = 'Unknow'
        do 201 i=1,6
          if (nfms.eq.rhkl(i)) cnfms=chkl(i)
 201    continue
        cseq  = 'Unknow'
        do 202 i=1,8
          if (seq.eq.rption(i)) cseq=cption(i)
 202    continue
        cpsifl = 'Unknow'
        do 203 i=1,8
          if (psifl.eq.rption(i)) cpsifl=cption(i)
 203    continue
C
C---- Use the values if possible
C
          nmcell = numang
          cwave  = flam(1)
C
C---- Write out the values for records 26, 27 and 28
C
        if (print) then
          write(*,'(/'' Data collection parameters:'')')
          write(*,'(/''   Text:'')')
          write(*,'(1X,A36)') dtext
          write(*,'(/''   Thlim - theta limits:'')')
          write(*,'(F11.2,'' degrees. Theta minimum'')') dthmin
          write(*,'(F11.2,'' degrees. Theta maximum'')') dthmax
          write(*,'(/''   Scan  - Scan parameters:'')')
          write(*,'(F11.2,'' degrees. Omega scan width param DOMA'')')
     1     doma
          write(*,'(F11.2,'' degrees. Omega scan width param DOMB'')')
     1     domb
          write(*,'(F11.2,'' degrees. Aperture width param APTA'')') apta
          write(*,'(F11.2,'' degrees. Aperture width param APTB'')') aptb
          write(*,'(I11,'' Scan type'')') iscant
          write(*,'(/''   Int   - intensity data collection '',
     1     ''parameters:'')')
          write(*,'(F11.3,'' Prescan acceptance parameter'')') sigpre
          write(*,'(F11.3,'' Final scan sigma(I)/I'')') sigma
          write(*,'(I11,'' Prescan speed parameter NPIPRE'')') npipre
          write(*,'(I11,'' Maxmimum time'')') itmax
          write(*,'(I11,'' Output profile flag'')') idumpf
          write(*,'(/''   Flag:'')')
          write(*,'(I11,'' Balanced filter unit use'')') nbalf
          write(*,'(I11,'' Measure Friedel pairs'')') mfridl
          write(*,'(I11,'' Non-equal test'')') neqfl
          write(*,'(/''   Psi   - Azimuthal setting parameters:'')')
          write(*,'(5X,A6,'' Psi data collection'')') cpsifl
          write(*,'(I11,'' Number of final scans'')') nrpsi
          write(*,'(F11.2,'' degrees. Starting psi value'')') psist
          write(*,'(F11.2,'' degrees. Psi increment'')') delpsi
          write(*,'(F11.6,'' Shape vector V1'')') dvect(1)
          write(*,'(F11.6,'' Shape vector V2'')') dvect(2)
          write(*,'(F11.6,'' Shape vector V3'')') dvect(3)
          write(*,'(/''   Mode  - Hkl sequence flags:'')')
          write(*,'(5X,A6,'' Sequence descriptor'')') cseq
          write(*,'(I11,'' Minimum h'')') indzz(1,1)
          write(*,'(I11,'' Maximum h'')') indzz(3,1)
          write(*,'(I11,'' Minimum k'')') indzz(1,2)
          write(*,'(I11,'' Maximum k'')') indzz(3,2)
          write(*,'(I11,'' Minimum l'')') indzz(1,3)
          write(*,'(I11,'' Maximum l'')') indzz(3,3)
          write(*,'(I11,'' Starting h'')') indzz(2,1)
          write(*,'(I11,'' Starting k'')') indzz(2,2)
          write(*,'(I11,'' Starting l'')') indzz(2,3)
          write(*,'(5X,A6,'' Sequence descriptor'')') cnfms
          write(*,'(/''   Absent- Systematic absences flags:'')')
          write(*,'(I11,'' Number of absent conditions'')') nabs
          write(*,'(18I4,'' LABS'')') labs
          write(*,'(/''   Intcr - Intensity control parameters:'')')
          write(*,'(I11,'' secs. Intensity control parameter'')') ichkt
          write(*,'(F11.3,'' Sigma(I)/I control value'')') sigmi
          write(*,'(F11.3,'' Intensity cut-off point'')') fading
          write(*,'(/''   Orient- Reorientation parameters:'')')
          write(*,'(I11,'' Orientation control parameter'')') norchk
          write(*,'(F11.2,'' degrees Max allowable deviation'')') dang
          write(*,'(/''   Backg - Background measurement times:'')')
          write(*,'(I11,'' Static background selection'')')
     1     nint(dints(1))
          write(*,'(F11.2,'' secs. Static background time'')') dints(2)
          write(*,'(F11.2,'' Background postion parameter'')') dints(3)
          write(*,'(/''   Current values of data collection:'')')
          write(*,'(I11,'' Pointer to next refln in SEPHKL file'')')
     1     nrhkl
          write(*,'(I11,'' Scan speed parameter for this reflection'')')
     1     npid
          write(*,'(I11,'' Intensity control counter'')') nics
          write(*,'(I11,'' Orientation control counter'')') nors
          write(*,'(I11,'' Azimuth parameter'')') npsis
          write(*,'(I11,'' Intensity sequence control'')') npifl
          write(*,'(I11,'' Reflection sequence no.'')') nnrefl
          write(*,'(F11.4,'' Current intensity'')') intca
          write(*,'(F11.1,'' X-ray exposure time'')') xrtime
          write(*,'(F11.2,'' degrees. Centre of theta scan'')') dtheta
          write(*,'(F11.2,'' degrees. Centre of omega scan'')') domk
          write(*,'(F11.0,'' Current h index'')') dhkl(1)
          write(*,'(F11.0,'' Current k index'')') dhkl(2)
          write(*,'(F11.0,'' Current l index'')') dhkl(3)
          write(*,'(F11.2,'' degrees. Current azimuth angle'')') datpsi
          write(*,'(F11.4,'' DINTS4'')') dints(4)
          write(*,'(F11.4,'' DINTS5'')') dints(5)
          write(*,'(F11.4,'' DINTS6'')') dints(6)
C
          write(*,'(/'' General orientation and scan parameters:'')')
          write(*,'(/'' Orientation:''
     1     /3f11.5,'' reciprocal r11 r12 r13'',
     2     /3f11.5,'' reciprocal r21 r22 r23'',
     3     /3F11.5,'' Reciprocal R31 R32 R33'')')
     4     r(1,1),r(1,2),r(1,3),r(2,1),r(2,2),r(2,3),
     5     r(3,1),r(3,2),r(3,3)
          write(*,'(F11.5,'' Determinant'')') det
          write(*,'(/3f11.5,'' Direct D11 D12 D13'',
     1     /3f11.5,'' Direct D21 D22 D23'',
     2     /3F11.5,'' Direct D31 D32 D33'')')
     4     d(1,1),d(1,2),d(1,3),d(2,1),d(2,2),d(2,3),
     5     d(3,1),d(3,2),d(3,3)
          write(*,'(/'' Wavelengths:'')')
          write(*,'(F11.5,'' Angstroms. Wavelength 1 '')') flam(1)
          write(*,'(F11.5,'' Angstroms. Wavelength 2 '')') flam(2)
          write(*,'(/'' Attenuator:'')')
          write(*,'(F11.5,'' Attenuator factor'')') zattn
          write(*,'(/'' Scan parameters:'')')
          write(*,'(F11.2,'' degrees. Omega scan param SWOMA'')')
     1     params(1)
          write(*,'(F11.2,'' degrees. Omega scan param SWOMB'')')
     1     params(2)
          write(*,'(F11.2,'' degrees. Aperture param DIAFRA'')')
     1     params(3)
          write(*,'(F11.2,'' degrees. Aperture param DIAFRB'')')
     1     params(4)
          write(*,'(F11.2,'' Scan ratio SCANR'')') params(5)
          write(*,'(F11.2,'' Quality factor QFACT'')') params(6)
          write(*,'(F11.2,'' Poisson Fraction PRFRACT'')') params(7)
          write(*,'(/'' Learnt profile information:'')')
          write(*,'(''   Record 28 of the file has information '',
     1      ''on the learnt profile.'')')
          write(*,'(''   The problem is that i do not have any idea '',
     1      ''of the format!''/
     2      '' See Chapter XI page 3 of the Cad4 Version 5 Manual.'')')
          write(*,'(//)')
        endif
 70     close(unit=parmfl,status='KEEP')
      endif
      return
      end
COTH
COTH   subroutine prflrd
C
C===== Subroutine to read a CAD4 .PAR and .CRY files
C
COTH  include 'common.for'
COTH  integer trim
C
C==== In fact this is dummy subroutine that does nothing. It is provided
C     for those who have to run DIFRAC on a machine other than a VMS
C     used to run the diffractometer. This is because the .PAR and .CRY
C     files are very VMS specific in their use of FORTRAN.
C 
COTH  write(*,'(''   File: gonca0.par not available'')')
COTH  write(*,'(''   File: '',A,''.cry not available'')')
COTH 1   compid(1:trim(compid))
COTH  return
COTH  end
