
                          READ ME FIRST
                          -------------

  The DIRAC software is supplied on the Xtal3.2 tape as a set of concatenated
  files. These are separated by a line starting with '-----cut here------'.
  To implement DIFRAC, first separate the files and follow the procedure
  given in the first file "Anciliary Material" following immediately below.


------cut here------------------------------------------------------------------




      DIFRAC Single-Crystal Diffractometer Output-Conversion Software
      ===============================================================

                By H.D. Flack, E. Blanc and D. Schwarzenbach


Date:    7th. October 1992

Contact: H.D. Flack  E-mail: flack@sc2a.unige.ch  or FLACK@CGEUGE52.bitnet

Reference:  J. Appl. Cryst. (1992) 25, 455-459

Distribution: By E-mail from H.D. Flack


                              Program Modules
                              ===============

Files            Content                      Update
=====            =======                      ======
ancill.mat    Ancillary material          7th. Oct. 1992
gn.for        General Module              7th. Oct. 1992
cf.for        CIF Output Module           7th. Oct. 1992
sf.for        SCFS Output Module          7th. Oct. 1992
c4v4.for      CAD-4 (v4.x) Module         7th. Oct. 1992
c4v5.for      CAD-4 (v5.0) Module         7th. Oct. 1992
niv411.for    Syntex/Nicolet/Siemens      7th. Oct. 1992
pw.for        PW1100 Module               7th. Oct. 1992
common.tlb    Common include Module       7th. Oct. 1992
address.lst   Sample address list         7th. Oct. 1992



                    Installation Instructions for DIFRAC
                    ====================================

1) Read the paper "DIFRAC, single-crystal diffractometer output-conversion
   software" by H.D. Flack, E. Blanc and D. Schwarzenbach (1992), 25, 455-459.
   (A copy of this paper in computer-transmissible form is included below)

2) Read the whole contents of this file

3) It is on the computer controlling the diffractometer that DIFRAC should be
   installed if at all possible.

4) For each diffractometer, one has to choose the appropriate source codes from
   which to make an executable module of DIFRAC. One needs:
   a)  the general module, gn.for,
   b)  an output module, i.e. either cf.for for CIF output or sf.for for
       SCFS 90 output,
   c)  a diffractometer module,
   d)  and the COMMON include file common.tlb .

5) Edit the source codes of the chosen general, output and diffractometer
   modules for installation dependent factors:
   a) Set the directory paths and logical unit numbers appropriate for your
      installation. Information is given in the section "Files and Logical Units
      Used by DIFRAC".
   b) Look at the open statements by searching for open( . It may be necessary
      to change these since they are very often installation dependent.
   c) If you are using CAD-4 version 5.0 with c4v5.for and you will NOT be
      running it on VAX/VMS, you will need to follow the following steps.
      Search for the subroutine prflrd by using the search word prflrd . Delete
      the existing subroutine between the line subroutine prflrd and its
      terminating line end. Just after the real subroutine prflrd (which will
      only work on VAX/VMS you will find a short dummy subroutine. Activate
      this dummy subroutine by changing COTH into 4 blank spaces (    ) globaly.

6) If you are using a Enraf-Nonius CAD-4 diffractometer, you will need to
   set the integer variable slidc4 in the choos subroutine in the diffractometer
   module. A value of zero for slidc4 works when DIFRAC is installed on the
   computer controlling the diffractometer. As Enraf-Nonius has chosen to have
   a rather unusual structure for their files, it often happens that the
   characters in the lines in the file become shifted by one character either
   to the left or to the right when the files are shifted to another machine.
   It is possible to compensate for this effect by setting slidc4 to +1 or -1.
   If an incorrect value is chosen, DIFRAC will not manage to read the
   information in the data file.

7) Compile and link into an executable module, the chosen edited source codes.

8) Make your own local address list file address.lst following the pattern of
   the distributed sample address list file.

-------------------------------------------------------------------------------

                    Files and Logical Units Used by DIFRAC
                    ======================================

Files for which a specific directory path is not given should be placed in the
directory from which the installed DIFRAC program is executed.

Definition of symbols:
----------------------

<compid>                     The 6 character compound identification code of
                             the studied compound. This code is set at the very
                             beginning of a DIFRAC run by question-and-answer.

<datxtn>                     The extension given to the primary (manufacturers')
                             diffractometer output data files. <datxtn> is set
                             by the variable datxtn in the subroutine choos of
                             the appropriate diffractometer module.


<addnam>                     Directory path and name of the address list file
                             for the diffractometer. <addnam> is set by the
                             variable addnam in the subroutine choos of the
                             appropriate diffractometer module.

<extn>                       The extension given to output standardized
                             diffractometer file. <extn> is set by the variable
                             extn in the subroutine set of the appropriate
                             output module.


Files:
------

<addnam>                     File containing a list of postal and e-mail
                             addresses, telephone and fax numbers of the local
                             users of DIFRAC. A sample list of the Geneva
                             installation is provided to show the format of this
                             file.

answer.sty                   File containing prepared (compound-independent, 
                             user-specific) answers for use by DIFRAC.

<compid><datxtn>             The primary (manufacturers') diffractometer output
                             data file when the measurement has produced a
                             single file.

<compid>nn<datxtn>           The primary (manufacturers') diffractometer output
nn = 01, 02, 03, ...         data files when the measurements have produced
                             several files.

<compid><extn>               The output standardized diffractometer data file.

<compid>.tth                 For the Philips PW1100 only: a diffractometer
                             parameter file containing h k l and omega values
                             of centered reflections.

gonca0.par                   For the CAD-4 version 5.x only:  the goniometer
                             constant file.

<compid>.cry                 For the CAD-4 version 5.x only: the crystal
                             parameter file.


<compid>_REF.DCF             For the Siemens P21 only: the crystal data file.

<compid>_PAR.DCF             For the Siemens P21 only: the crystal parameter
                             file.

<compid>.qao                 File containing answers typed in reply to questions
                             asked in the penultimate run of DIFRAC for the 
                             compound <compid>.

<compid>.qan                 File containing answers typed in reply to questions
                             asked in the latest run of DIFRAC for the compound
                             <compid>.



Logical Units:
--------------

The program uses extensively the standard input and output files as accessed
by the FORTRAN instructions read(*,.. and write(*,.. . The following logical
unit integer variables are set in the subroutine choos of the appropriate
diffractometer module.

     datafl is the unit used to read the diffractometer data files
     scfsfl is the unit used to write the standardized output file
     parmfl is the unit used to read the diffractometer parameter files
     addrfl is the unit used to read the address list file
     qaofl  is the unit used to write/read an existing saved answer file
     qanfl  is the unit used to read/write a new saved answer file

-------------------------------------------------------------------------------

                            Saved Answer Procedure
                            ======================

  A system of saving the answers given to the questions asked by DIFRAC has been
installed. This makes it possible to vastly speed up the use of DIFRAC if one
has to rerun the program. Expert users might like to use this system to reduce
the time necessary to run the question-and-answer session even in the case of
a single run of DIFRAC.

  In every run of DIFRAC during the question-and-answer session, answers to
questions are written to a file <compid>.qan , where <compid> is the compound
identification of the compound studied. On each line of this file is an integer
code of 6 digits identifying the question, followed by up to 80 characters as
typed in reply to the question. The table below gives the correspondence
between the integer code and the question asked.

  At the beginning of a run, DIFRAC looks to see if a file called <compid>.qan
is present:
a) If <compid>.qan is present, the contents of this file are copied to a file
   <compid>.qao
b) If <compid>.qan is not present, DIFRAC looks to see if a file answer.sty is
   present. If answer.sty is present, its contents are copied to <compid>.qao .
  During the question-and-answer session, saved answers are retrieved from the
newly created <compid>.qao file and displayed to the user with a message
** Saved answer - Return to accept. It is thus only necessary to type a
"Return" to reply to the question with the saved answer.

  For expert users wishing to create a-priori or edit existing <compid>.qan
files, it is essential that the answers in this file appear in ascending order
of their integer code. The integer code must occupy columns 1 to 6 of the file
and the answer commences directly in column 7.

 000020  Give a title - max 66 characters:
 000040 Do you wish to display the contents of this file (Y,N)? gonca0.par
 000060 Do you wish to display the contents of this file (Y,N)? <compid>.cry
 000080 Do you wish to display the contents of this file (Y,N)? <compid>_PAR.DCF
 010020  Do you have refined cell dimension data (Y,N)?
 010040  Enter a, b, c, alpha, beta, gamma in Angstroms and degrees:
 010060  Enter e.s.d.s in a, b, c, alpha, beta, gamma:
 010080  Enter cell volume and e.s.d. in cubic Angstroms:
 010100  How were these cell dimensions measured?  enter wavelength,
         temperature(k), 2theta-min, 2theta-max, number of reflections:
 010120  Cell dimension data OK (Y,N)?
 050020  Have you measured the density (Y/N)?
 050040  Enter density, e.s.d. and temperature(K) of measurement:
 050060  Density data OK (Y,N)?
 100020  Make a choice for the shape of the crystal:
 100040  Give sphere radius and e.s.d in mm:
 100055  How many indexed faces does the crystal have?
 100060  Input one face per line: h,k,l,d,esd(d).
         d & esd(d) in mm.  Finish with CTRL Z:
 100080  Crystal shape OK (Y,N)?
 150020  Enter Z, the number of formula units per cell:
 150040  Enter chemical formula
 150060  <chemical formula second line>
 150080  Do you have a chemical analysis (Y/N)?
 150100  Wt % & esd of element:
 150120  Do you have an observed value of the formula weight (Y/N)?
 150140  Enter observed formula weight and e.s.d.(Daltons):
 150160  Chemical formula data OK (Y,N)?
 200020  Enter the melting-point temperature of the compound in K:
 200040  Enter a description of the source of the chemical compound
         under study. If not known leave blank:
 200060  Enter the trivial name of the compound. If not known leave blank:
 200080  Enter the mineral name of the compound
 200100  Chemical compound data OK (Y,N)?
 250020  Enter the colour of the crystal:
 250040  Enter a description of the crystal quality and habit
 250060  Enter special details concerning the preparation or mounting of
         the crystal otherwise blank:
 250080  Crystal data OK (Y,N)?
 300020  Filter index xxx:
 300040  Diffractometer attenuator data OK (Y,N)?
 350020  Did you use a beta-filter(b) or a monochromator(m) (B,M)?
 350040  Have you a measured value of the monochromator polarization
         ratio (Y,N)?
 350060  Give polarization ratio and e.s.d.:
 350080  Give monochromator 2-theta angle in degrees:
 350100  Have you a measured value of the monochromator inhomogeneity
         half-width (Y,N)?
 350120  Give inhomogeneity half-width(mm):
 350140  Give beta-filter absorption edge wavelength(a):
 350160  Diffractometer monochromator data OK (Y,N)?
 400020  Have you a measured value of the counter dead time (Y,N)?
 400040  Give measured dead-time and e.s.d. in microsecs:
 400060  Counter dead time data OK (Y,N)?
 450020  Enter temperature (k) and pressure during intensity measurement:
 450040  Environmental conditions data OK (Y,N)?
 500020  Enter the number of your name in the list:
 500040  Enter your name as: family name(s) first name(s) or initials:
 500060  Enter your address (5 lines):
 500080  Enter your e-mail address:
 500100  Enter your telefax number:
 500120  Enter your telephone number:
 500140  Name address data OK (Y,N)?
-------------------------------------------------------------------------------

      DIFRAC, Single-Crystal Diffractometer Output-Conversion Software
      ================================================================

      J. Appl. Cryst. (1992) 25, 455-459

By H.D. Flack
Laboratoire de Cristallographie, University of Geneva, 24 quai Ernest-Ansermet,
CH-1211 Geneva 4, Switzerland
Tel: [+ [41] 22] 702 62 49  Fax: [+ [41] 22] 781 21 92
E-mail: flack@sc2a.unige.ch

E. Blanc and D.Schwarzenbach
Institut de Cristallographie, University of Lausanne, BSP Dorigny,
CH-1015 Lausanne, Switzerland

Abstract
========
  Software is described which converts single-crystal diffractometer output
files as produced by manufacturers' software into a standardized instrument-
independent form consisting of a clear, complete, well-documented record of the
sample and the diffraction measurements performed upon it. Information not
already available in the manufacturers' diffractometer files is obtained in an
interactive question-and-answer session. The software is written in a modular
way. Available modules can treat Enraf-Nonius CAD-4, Philips PW1100 and
Siemens P21 single-crystal diffractometers, and produce output in CIF or SCFS
format.


Introduction
============
  Users of several models of 4-circle single-crystal diffractometers may well
have been struck by the diversity of form and content of the data files
generated by diffractometer-manufacturers' software. The form of the file can
create difficulties in its transfer to other computing equipment and further
renders the corresponding computer programs specific to a certain type or
types of diffractometer. The difficulties in the content of the files manifest
themselves principally by the paucity of the available information
necessitating additional input to the data treatment software (e.g. type of
radiation, wavelength of radiation, scan width, sense of rotation of angles on
the single-crystal diffractometer, etc.). The problem has become
aggravated in recent times by the rapid development in electronic data
exchange. It is now frequent to exchange files by computer networking and
electronic mail rather than by the standard means of the recent past, the
1/2in magnetic tape. The advent of machine-readable submission for publication,
data-base and supplementary-material deposition further highlights the problem
of missing data.
  Our own motivation for producing this software arose from the work of the
IUCr Subcommittee on Statistical Descriptors (Schwarzenbach, Abrahams 
et al., 1989). We wished to be able to perform structure refinements on the 
directly-observed quantities such as peak and background intensities 
(Schwarzenbach & Flack, 1991), crystal dimensions (Schwarzenbach, 1991; 
Blanc, Schwarzenbach & Flack, 1991) etc. which should not be corrected for any 
systematic effects. In this project we also had to face the problem of 
producing refinement and data-treatment software which would function smoothly 
on data produced by three different types of single-crystal diffractometer. 
Further it became clear that with the advent of machine-readable submission of 
structural papers to journals and data bases it would be highly advantageous to 
obtain a clear, complete, instrument-independent record of the diffractometer 
measurements and conditions in the form of a data file at the earliest possible 
opportunity i.e. while the crystal was still mounted on the instrument.

Storing Raw Intensities
=======================
  Raw observed intensity data are output to the converted file. In the current
version of the program we have dealt only with low-angle background -
integrated peak - high-angle background counts, thus leaving the problem of the
transfer of profile data to a later date. Conversion of raw counts to net
intensities Pobs (i.e. with background subtracted and allowing for differences
in effective scan speed) or to |Fobs|**2 values (i.e. applying a Lorentz-
polarization correction to the net intensity) is not undertaken by this
software. As the raw data are conserved in the converted file, a more detailed
modelling and refinement of background than that normally undertaken at present
will be feasible, whilst retaining the possibility of a classic data reduction.
As manufacturers' diffractometer-control software specifies and records its
scan conditions (i.e. scan time, scan width, scan speed and number of scans)
in many different forms, we have chosen to standardize this information into
four intensity-measurement coefficients cp, cb, cl, and ch most convenient for
least-squares calculations.
  Let Qobs be the observed integrated-intensity count measured by scanning
over the reflection profile, and Lobs and Hobs the observed low- and high-angle
background counts respectively. The model relating the raw observations of a
diffraction experiment to derived parameters may be written as
Qcalc = cp Pcalc (vm) + cb Bcalc ~ Qobs, Lcalc = cl Bcalc ~ Lobs and
Hcalc = ch Bcalc ~ Hobs. Pcalc is the calculated net integrated intensity per
unit inverse scan speed and is a function of vm adjustable parameters
describing the crystal structure, the scale factor, secondary extinction and
possibly other physical quantities. Bcalc is the calculated count rate per unit
time of the mean background and is itself an adjustable parameter. One can see
that in general, cp is equal to the time to count the peak divided by the scan
width i.e. the effective scan speed allowing for multiple scans over the same
peak, cb to the time to count the peak, cl to the time to count the low-angle
background and ch to the time to count the high-angle background. It follows
that the scan width is equal to cb / cp.
  For a classical treatment, the net intensity Pobs may be calculated from the
raw counts by using Pobs = [Qobs - 0.5 cb (Lobs / cl + Hobs / ch)] / cp.

System Design
=============
  Table 1 presents the list of items that we considered essential to include
in the converted single-crystal diffractometer output file. This information is
intended to be complete and independent of the particular instrument used. It
serves for the calculations of structure determination and refinement and 
contains the items considered necessary for machine submission of documents for 
publication or data bases. Some of the items we have chosen to include or 
exclude from our output list require special comment:
  (a) The chemical formula has been included in the converted output file
although in some cases this might have to be revised as a result of the crystal
structure analysis. The formula is useful in identifying the compound, and the
element names and composition are necessary for running structure-solution and
refinement software. Moreover, the subroutine (see below) which undertakes the
interpretation of the chemical formula was already available to us.
  (b) On the other hand, the space group of the crystal has not been included.
For a crystal of unknown structure, the exact space group will not be
determined reliably until after structure refinement. Further it was thought
that flexible software to test the integrity of the input of space group
information in any form (e.g. a symbol, generators, or individual symmetry
operations) would have greatly lengthened the time necessary to produce this
version of DIFRAC.
  (c) The crystal orientation matrix (or matrices in the case of resetting) for
the diffractometer (often called the UB matrix) and the diffractometer setting
angles (omega, two-theta, chi, phi or kappa) for each reflection have not been
included in the converted output file. All orientation information has been
condensed into a crystal-based azimuthal angle psi for each reflection
measurement. This psi is independent of the instrument geometry and axial
definitions of the instrument. Full details are to be found in Schwarzenbach &
Flack, (1989, 1992).
  (d) An index or code is included to identify the reference reflections. In
our view it is not possible to correctly identify a reference reflection from
its Miller indices alone. The inclusion of an index eases the work and
programming of data-treatment and refinement software analysing the variations
in the reference reflections.
  (e) The file includes a fairly detailed description of the incident-beam
characteristics. This has been done to allow the calculation of polarization
and beam-inhomogeneity corrections (Harkema, Dam, Van Hummel & Reuvers, 1980;
Markov, Fetisov & Zhukov, 1990).
  A choice had to be made concerning the format of the converted output file.
The necessity of being able to transmit the file over networks or by electronic
mail imposes a clear-text (ASCII) form with no more than 80 characters per line.
Additionally we required that the content and format of the file be clearly 
defined and relatively compact. Both the SCFS-90 - Standard Crystallographic 
File Structure-90 (Brown, 1988) and the more recent CIF - Crystallographic
Information File (Hall, Allen & Brown, 1991) are available for just this
purpose and we have adopted both of them.
  The information to be included in the converted output file is obtained from
the following sources: (a) fixed geometrical details of the single-crystal
diffractometer, (b) information in the diffractometer parameter file,
(c) information in the diffractometer data file and (d) interactive questions
and answers. An underlying design principle in this arrangement is that
information from the diffractometer files should be used, if available,
otherwise a question will be asked just once only. This is most important to
avoid confusion and frustration on the part of the user. It further explains
our preference for installing the executable programs on the computer
controlling the diffractometer so that each installation may be customized for
its instrument and have access to those files which are not specific to the
user.

DIFRAC
======
  This is a set of free-standing programs written in FORTRAN-77 consisting of
modules having precisely defined functions described below. Entry points to
each module and cross-reference calling and linking information are given in
Table 2. Each program is diffractometer and output specific and is generated
from the general module (GN.FOR), from one of the diffractometer modules
(C4V4.FOR, C4V5.FOR, NIV411.FOR, PW.FOR) , from one of the output modules
(CF.FOR, SF.FOR) and the included file of COMMON definitions (COMMON.TLB). Due
to the modular construction of DIFRAC, it is easy to produce additional
programs for other types of single-crystal diffractometer and output styles
by writing diffractometer-specific or output-specific modules having the same
entry points and functionality as the existing modules. The salient features
of the modules and individual implementations are given below.

General Module
--------------
  This contains the main controlling section of the program and the
interactive free-format question-and-answer subroutine. DIFRAC uses a file
ADDRESS.LST containing the names, addresses, telephone numbers, telefax numbers
and e-mail addresses of the users on the local installation. The subroutine
CHFORM for the conversion and verification of chemical formulae is adapted from
the program CASTOR (Brett, Miebach, Boese, & Eyl, 1990) with the authors'
permission.

Diffractometer Modules
----------------------
  Currently, modules are provided for Philips PW1100, Enraf-Nonius CAD-4
(versions 4.x and 5.0) and Syntex/Nicolet/Siemens P21 (versions 4.11 or lower)
single-crystal diffractometers.

  PW1100 - Philips: This is a modified system, similar to others described in
the literature, (Grigg & Barnea, 1990; Van Hummel & Graafsma, 1989) which no
longer uses the original paper tape, magnetic tape or Teletype output. The
Teletype has been replaced by an IBM-compatible PC running KERMIT which
captures all of the output into a data file of name <compid>.DAT
(or <compid>01.DAT, <compid>02.DAT, ....., <compid>nn.DAT for a sequence of nn
data files on the same compound). The .DAT file should start with the output of
the PW1100 instruction PAR followed by the intensity collection itself. A
parameter file called <compid>.TTH may also be prepared containing the results
of the accurate centring of reflections intended for cell-dimension
determination. For compatibility with other software, this file should start
with one line which may contain anything, followed by one line containing an
integer code defining the default wavelength of the centred reflections,
followed by lines containing h, k, l, +omega, -omega in degrees, (the latter
being the omega values of the centring of the reflection at positive and
negative two-theta) and the wavelength code of the measurement. The wavelength
codes are defined as follows: 0 or blank indicates a default value, 1 means 
alpha1 line, 2 means alpha2 line and any other value indicates the mean alpha 
value. Presently the DIFRAC software has been installed on a DEC VAX 8700/VMS. 
We have not yet attempted to install it on the IBM-compatible PC serving as 
interface to the diffractometer.

  CAD-4 - Delft Instruments (Enraf - Nonius). Modules for versions 4.x and 5.0
of the manufacturer's software are available. The module for version 4.x reads
only a crystal data file <compid>.DAT  whereas that for version 5.0 can read
the goniometer constant file GONCA0.PAR, the crystal parameter file
<compid>.CRY and the crystal data file <compid>.DAT (or <compid>01.DAT,
<compid>02.DAT, ....., <compid>nn.DAT for a sequence of nn data files on the
same compound). The DIFRAC software has been installed on the DEC MicroVAX
3100/VMS controlling the diffractometer running version 5.0 of the control
software.

  P21 Siemens (Previously Syntex - Nicolet). This module was written for a
machine running version 4.11 of the control software but takes account of the
specifications of version 4.12 and higher, up to but not including version 5.1.
DIFRAC can read the crystal parameter file <compid>_PAR.DCF, the crystal data
file <compid>_REF.DCF and has been installed on the MicroVax II controlling
the diffractometer. There are a few known restrictions in the code. The P-1
version of the SYNTEX machine had the monochromator oriented at 90 degrees to 
the current arrangement and this condition is not tested for by the software. 
(It would need to be set manually by the variable FMNHED in the CHOOS 
subroutine.) Further, the current software does not treat data measured with 
Wyckoff scans.

Output Modules
--------------
  In the current version of DIFRAC, modules are provided for CIF (Hall, Allen
& Brown, 1991) and SCFS90 output (Brown, 1988). Both the FORTRAN source code of
the CIF output module and the CIF output files it produces have been checked
for conformity with the CIF definitions contained in the CIF dictionary (Core
version 1991) as contained in the electronic file cifdic.c91 by using the
checking program CYCLOPS of Hall (1991). The only non-standard names are
_diffrn_refln_fbcoeff_peak_z, _diffrn_refln_fbcoeff_bg_peak_z,
_diffrn_refln_bg_1_z and _diffrn_refln_fbcoeff_bg_2_z for the intensity-
measurement coefficients cp, cb, cl and ch defined above and
_cell_measurement_refln_wl_id_z giving the wavelength code of each reflection
centred for cell dimension measurement.
  In the XTAL 3.0 system (Hall & Stewart, 1990) the SCSF 90 file is read by the
program REFM90 (Flack,1990) to form an XTAL bdf which may then be processed
by STARTX in update mode. Following STARTX the treatment of the intensity
measurements is undertaken by REFCAL (Flack, 1992). This arrangement has been
functioning in the Geneva laboratory to our thorough satisfaction for about 1
year with the PW1100 and CAD-4 diffractometers. For the input of a CIF file into
XTAL, the same type of arrangement will use the program CIFIO (Hall, 1990).
At the present time the reading of the CIF file is much less reliable than the
mature SCFS-file-reading software.

Concluding Remarks
==================
  It must be emphasised that the program modules described in this paper are
no more than a stop-gap solution on the way to standardised, complete,
documented and easily transmissible single-crystal diffractometer output files.
There can be no substitute for the diffractometer-control software itself being
upgraded to produce such data files. Nevertheless the conversion of
diffractometer output by way of DIFRAC is a useful step in the right direction.
Program packages for structure solution, refinement and publication can thus
be updated to become instrument independent with all information on the
experiment coming from the data file itself. The user interaction with his
favourite software can thus be reduced entirely to the task of specifying the
calculations to be undertaken rather than being a mixture of data and
operational commands.
  The current paper is the prime source of documentation concerning DIFRAC. A
computer file of ancillary material contains updates and corrections to this
paper. The source code and ancillary material of DIFRAC are available from
flack@sc2a.unige.ch or flack@cgeuge52.bitnet by e-mail.

References
==========
BLANC, E., SCHWARZENBACH, D. & FLACK, H. D. (1991). J. Appl. Cryst. 24,
  1035-1041.
BRETT, W.A., MIEBACH, T., BOESE, R., & EYL, J. (1990). CASTOR. University of
  Bonn.
BROWN, I.D. (1988). Acta Cryst. A44, 232.
FLACK, H.D. (1990). REFM90: XTAL 3.0 Crystallographic Program System (Eds.
  S.R. Hall & J.M. Stewart) Universities of Western Australia and Maryland.
FLACK, H.D. (1992). REFCAL: XTAL 3.2 Crystallographic Program System (Eds.
  S.R. Hall & J.M. Stewart) Universities of Western Australia and Maryland. In
  preparation.
GRIGG, M.W. & BARNEA, Z. (1990). J.Appl. Cryst. 23, 432-433.
HALL, S.R. (1990). CIFIO: XTAL 3.0 Crystallographic Program System (Eds.
  S.R. Hall & J.M. Stewart) Universities of Western Australia and Maryland.
HALL, S.R. (1991) CYCLOPS: A Program for validating data names in a text file.
  University of Western Australia.
HALL, S.R., ALLEN, F.R. & BROWN, I.D. (1991). Acta Cryst. A47, 655-685.
HALL, S.R. & STEWART, J.M. (1990). Eds. XTAL 3.0 Reference Manual. Universities
  of Western Australia and Maryland.
HARKEMA, S., DAM, J. VAN HUMMEL, G.J. & REUVERS, A.J. (1980). Acta Cryst. A36,
  433-435.
HUMMEL, G.J. VAN  & GRAAFSMA, H. (1989), J. Appl. Cryst. 22, 75.
MARKOV, V.T., FETISOV, G.V. & ZHUKOV, S.G. (1990). J.Appl. Cryst. 23, 94-98.
SCHWARZENBACH, D. (1991). In: Crystallographic Computing 5: From Chemistry to
  Biology, edited by  D. MORAS, A.D. PODJARNY & J.C. THIERRY, Oxford University
  Press
SCHWARZENBACH, D., ABRAHAMS, S.C., FLACK, H.D., GONSCHOREK, W., HAHN, TH.,
  HUML, K., MARSH, R.E., PRINCE, E., ROBERTSON, B.E., ROLLETT, J.S. &
  WILSON, A.J.C. (1989). Acta Cryst. A45, 63-75.
SCHWARZENBACH, D. & FLACK, H.D. (1989). J. Appl. Cryst. 22, 601-605.
SCHWARZENBACH, D. & FLACK, H.D. (1991). Acta Cryst. A47, 134-137.
SCHWARZENBACH, D. & FLACK, H.D. (1992). J. Appl. Cryst. 25, 69.


Table 1: Content of the Converted Output File
=============================================

Compound, measurement and experimenter identification
  General descriptive title
  Compound/measurement identification code
  Date and method of creation of the converted output file
  General experimental technique and special details
  Name and postal, e-mail, telephone and telefax addresses of the experimenter

Instrument specification
  Type of single-crystal diffractometer
  Specification of software used to control the instrument
  Detector dead time and e.s.d.

Radiation description
  Type of radiation (i.e. neutrons, electrons, X-rays)
  Mean wavelength
  Wavelengths and intensity weights of radiation components
  Polarization ratio and e.s.d.
  Angle between the normals to the diffraction planes of the monchromator and
    sample
  Incident beam spatial half-width at the crystal

Filter and attenuator information
  Wavelength of beta filter absorption edge
  For each attenuator filter: its reciprocal transmission factor with e.s.d and
    index

Environmental conditions
  Temperature with e.s.d.
  Pressure with e.s.d.

Compound and crystal specification
  Measured density with e.s.d. and temperature
  Melting point
  Colour and description of the crystal
  Source
  Common chemical and mineral names
  Measured formula mass with e.s.d.
  Chemical formula
  Chemical analysis in mass%
  Number of formula units per cell

Cell dimension and its determination
  Cell dimensions with e.s.d.s
  Cell volume with e.s.d.
  Temperature
  Specification of radiation used
  Minimum and maximum two-theta of reflections
  Number of reflections
  h, k, l, two-theta and wavelength code of each reflection

Crystal-shape description: either (a), (b) or (c)
  (a) h, k and l of each face and its distance with e.s.d. from centre
  (b) Radius with e.s.d. of a spherical crystal
  (c) Radius and length with e.s.d.s of a cylindrical crystal

Intensity measurement procedure
  Minimum and maximum sin(theta)/lambda
  Total number of reflections measured
  Minimum and maximum values of h, k and l

Reference-reflection description
  Time interval between measurement of reference reflections
  Number of reference reflections measured
  Foe each reference reflection: its h, k and l Miller indices and reference
    code

Reflection-specific information
  Miller indices h, k and l
  Counts: peak, low- and high-angle background counts (Qobs, Lobs and Hobs)
  Intensity-measurement coefficients: cp, cb, cl and ch (defined in text)
  Crystal-based azimuthal angle
  Elapsed time of measurement
  Attenuator-filter index
  Scale-factor index
  Reference-reflection index
  Background scanning mode
  Scan type
  Total horizontal and vertical detector aperture


Table 2: DIFRAC Construction
============================

General Module (GN.FOR)
  Subroutine  Called by    Use
  DIFRAC                   Main Program
  BLOCKD                   Block data constants and initialisations
  CHFORM      QUESAN       Conversion of chemical formula into components
  COMPNM      DIFRAC       Ask for compound name
  FREAD       Many         Free format interactive read
  QUESAN      DIFRAC       Interactive question and answer
  TRIM        Many         Trim non-significant blanks from strings

Output Modules
  File Name   Output type
  CF.FOR      CIF - Crystallographic Information File
  SF.FOR      SCFS 90 - Standard Crystallographic File Structure - 90

  Subroutine  Called by    Use
  BASOUT      DIFRAC       Output all other information
  CONOUT      INITLZ       Output for changed conditions
  DTFLND      DIFRAC       Finish off converted output file
  HEDOUT      DIFRAC       Output starting section
  REFOUT      REFRD        Output reflection information
  SET         DIFRAC       Set up output style
  VALOUT      BASOUT       Convert value and e.s.d. into string

Diffractometer Modules
  File Name   Diffractometer type
  C4V4.FOR    Enraf-Nonius CAD-4 (Version 4.x)
  C4V5.FOR    Enraf-Nonius CAD-4 (Version 5.0)
  NIV411.FOR  Nicolet P21 (Version 4.11 and earlier)
  PW.FOR      Philips PW1100

  Subroutine  Called by    Use
  CHOOS       DIFRAC       Set up diffractometer type options
  DTFLRD      DIFRAC       Read diffractometer data file
  INITLZ      DTFLRD       Interpret non-intensity data
  JUNPAK      PRFLRD       Packing function for CAD-4
  PRFLRD      DIFRAC,FREAD Read diffractometer parameter file
  PSIC        REFRD,INITLZ Crystal-based azimuth calculation
  REFRD       DTFLRD       Interpret diffractometer intensity datum

Include Module (COMMON.TLB)
  Include file of COMMONs used by all modules

Linking Information
  Program     Modules Needed
  CAD4V4CF    GN, C4V4, CF
  CAD4V4SF    GN, C4V4, SF
  CAD4V5CF    GN, C4V5, CF
  CAD4V5Se    GN, C4V5, SF
  NICOCe      GN, NIV411, CF
  NICOSF      GN, NIV411, SF
  PW11CF      GN, PW, CF
  PW11SF      GN, PW, SF



------cut here------------------------------------------------------------------


C
C==== GENERAL Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
C
      program difrac
C---------------------------------------------------------------------------C
C**** Convert PW1100, CAD4, STADI-4 or P21  output into                     C
C     an SCFS 90 file.                                                      C
C---------------------------------------------------------------------------C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
      integer jj
C
      include 'common.tlb'
C
C---- Set type of Output
C
      call set
C
C---- Choose fixed diffractometer parameters
C
      call choos
C
C---- Enquiry of compid and check data file existence
C
      call compnm
C
C---- Read the parameter file
C
      call prflrd
C
C---- Output Header records to output file
C
      call hedout
C
C---- Loop over the individual data collection files
C
      do 9005 jj = sdtfil, ndtfil
        call dtflrd(jj)
        call dtflnd
 9005 continue
C
C---- Interactive question and answer section for data values
C
      call quesan
C
C---- Output the basic non-intensity information
C
      call basout
C
C---- Close the Output file
C
      close(unit=scfsfl, status='KEEP')
C
      stop ' '
      end
C
      subroutine compnm
C
C---- Get compound id and check for data file existence
C
      integer trim
      integer  i, j, n
      logical  zresnt, valid, lansw
C
      include 'common.tlb'
C
C---- Get the COMPOUND ID
C
 100  valid = .false.
      compid = ' '
      write(*,'('' Give the Compound ID - max 6 characters:'')')
      call fread(0,'CNNNNNNNNN')
      compid = charin(1)(1:6)
C---- Set up DSK
      dsk = compid(1:4)
C
C---- Make sure there are data files of this name
C
      write(*,'(/'' Data files for compound: '',A6)') compid
C---- Try a <compid> data  file first
      inquire(file= compid(1:trim(compid))//datxtn, exist = zresnt)
      if (zresnt) then
C---- One data file of name <compid>//datxtn
        sdtfil = 0
        ndtfil = 0
        valid = .true.
        write(*,'(''   File: '',A,A8)') compid(1:trim(compid)),datxtn
      else
        sdtfil = 1
        ndtfil = 0
 200    i = (ndtfil+1)/10
        j = mod((ndtfil+1), 10)
        inquire(file= compid(1:trim(compid))//digit(i)//digit(j)
     1   //datxtn, exist = zresnt)
        if (zresnt) then
          ndtfil = ndtfil + 1
          valid = .true.
          write(*,'(''   File: '',A,2I1,A8)')
     1     compid(1:trim(compid)),i,j,datxtn
        endif
        if (zresnt) go to 200
      endif
      if (.not. valid) then
        write(*,'(''   There are no data files for this name.'')')
        write(*,'(''   Y to reinput compound ID, N to stop (Y,N)? '')')
        call fread(0,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'N') stop ' '
      endif
      if (.not. valid) go to 100
C
C---- Open the Output File
C
      open( unit=scfsfl, file = compid(1:trim(compid))//extn,
     1 status = 'UNKNOWN', access = 'SEQUENTIAL', form = 'FORMATTED'
     2 )
CVAX 2 ,carriagecontrol='LIST')
CUNX 2 )
C
C---- Set up option to go high speed through the saved answers
C
      hspeed = .false.
      write(*, '(//'' Is ancilliary data to be entered without'',
     1 '' user interaction (Y,N)? '')')
      call fread(0,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') hspeed = .true.
C
C---- File control for saved answers.
C
C---- Is either a compound specific or a general saved answer file there?
      inquire(file= compid(1:trim(compid))//'.qan', exist = lquan)
      lansw=.false.
      if (.not.lquan) then
        inquire(file= 'answer.sty', exist = lansw)
        if (lansw) lquan=.true.
      endif
      if (lquan) then
C---- Copy existing saved answer file to <compid>.qao
        open( unit=qaofl, file = compid(1:trim(compid))//'.qao',
     1  status = 'UNKNOWN', access = 'SEQUENTIAL', form = 'FORMATTED'
     2  )
CVAX 2  ,carriagecontrol='LIST')
CUNX 2 )
        if (lansw) then
          open( unit=qanfl, file = 'answer.sty',
     1    status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     2    )
CVAX 2    ,readonly)
CUNX 2    )
        else
          open( unit=qanfl, file = compid(1:trim(compid))//'.qan',
     1    status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     2    )
CVAX 2    ,readonly)
CUNX 2    )
        endif
        do 210 i=1,10000
          read (qanfl,'(i6,a80)', end=220) n, qline
          write(qaofl,'(i6,a)') n, qline(1:trim(qline))
 210    continue 
 220    rewind  qaofl
        nqaocd = -1
        close(qanfl, status='KEEP')
      endif
C---- Open the new saved answer file
      open( unit=qanfl, file = compid(1:trim(compid))//'.qan',
     1 status = 'UNKNOWN', access = 'SEQUENTIAL', form = 'FORMATTED'
     2 )
CVAX 2 ,carriagecontrol='LIST')
CUNX 2 )
C
C---- Get the TITLE
C
      if (ask(4)) then
        write(*,'(//'' Give a title - max 66 characters:'')')
        call fread(000020,'LNNNNNNNNN')
        title = line(1:66)
      endif
      return
      end
C
      subroutine quesan
C---------------------------------------------------------------------------C
C     Question and answer subroutine for interactive input of diffractometerC
C     information.  The sections of the questions which are asked in the    C
C     interactive session are controlled by the Boolean array Ask.  The     C
C     values in Ask are set up in SET and CHOOS, so the questions are       C
C     tailored to the information available (or not) from the diffractometerC
C     and for the output SCFS or another.                                   C
C                                                                           C
C     Ask(2)    Cell dimensions of crystal.  Not needed by Nicolet          C
C     Ask(3)    Deadtime. May not be needed by CAD4 5.0                     C
C     Ask(4)    TITLE.  Not needed by Nicolet                               C
C     Ask(6)    Attenuator filter values. Only need to enquire if filters   C
C               have really been used during the data collection and the    C
C               attenuator coefficients not found in the files.             C
C                                                                           C
C---------------------------------------------------------------------------C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
      integer i, k
      logical oka
      integer trim
C
      include 'common.tlb'
C
C---- Input of cell dimensions
C
      if (ask(2)) then
 2005   write(*,'(//
     1  '' ***** Input of the cell dimensions of the crystal'')')
        write(*,'(
     1  '' Do you have refined cell dimension data (Y,N)? '')')
        call fread(010020,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'N') then
          do 2011 i = 1, 6
            cell(i) = cellub(i)
            ecell(i) = 0.01*cell(i)
 2011     continue
          volum = volub
          evolum = 0.03*volum
          write(*,'('' No refined cell dimension data available.'')')
          write(*,'(//'' Cell dimensions from ub matrix:''/12x,''a'',
     1    7x,''b'',7x,''c'',3x,''alpha'',4x,''beta'',3x,''gamma''/
     2    '' cell'',6f8.2//7x,''volume'',f8.0)') cellub, volub
        endif
        if (inchar .eq. 'N') go to 2033
        if (inchar .ne. 'Y') go to 2005
C---- Read the cell dimensions
        write(*,'('' Enter a, b, c, alpha, beta, gamma '',
     1  ''in Angstroms and degrees:'')')
        call fread(010040,'RRRRRRNNNN')
        do 2020 i = 1, 6
 2020   cell(i) = realin(i)
        write(*,
     1  '('' Enter e.s.d.s in a, b, c, alpha, beta, gamma:'')')
        call fread(010060,'RRRRRRNNNN')
        do 2025 i = 1, 6
 2025   ecell(i) = realin(i)
        write(*,
     1  '('' Enter cell volume and e.s.d. in cubic Angstroms:'')')
        call fread(010080,'RRNNNNNNNN')
        volum = realin(1)
        evolum = realin(2)
        write(*,
     1  '('' How were these cell dimensions measured?  Enter''/
     2  '' wavelength, temperature(k), 2theta-min, 2theta-max,'',
     3  '' number of reflections:'')')
        call fread(010100,'RRRRINNNNN')
        cwave = realin(1)
        ctemp = realin(2)
        ctthn = realin(3)
        ctthx = realin(4)
        nmcell = intin(5)
C---- End of cell dimension input
        write(*,'(//'' Cell dimension information:''/12x,''a'',
     1  7x,''b'',7x,''c'',3x,''alpha'',4x,''beta'',3x,''gamma''/
     2  '' cell'',6f8.4/''  esd'',6f8.4//7x,''volume'',5x,''esd''/
     3  5x,2f8.2//8x,''wavel'',4x,''temp'',3x,''2tmin'',3x,
     4  ''2Tmax'',4x,''Nref''/5x,f8.5,3f8.1,i8)')
     5  cell, ecell, volum, evolum, cwave, ctemp, ctthn, ctthx, nmcell
        oka = .true.
        do 2032 i = 1 , 3
          if(abs(cell(i)-cellub(i)).gt.0.2) oka = .false.
          if(abs(cell(i+3)-cellub(i+3)).gt.2.0) oka = .false.
 2032   continue
        if(.not. oka) then
          write(*,'(/'' ***** Error: Input cell dimensions are too '',
     1    ''different from those calculated ''/13x,''from ub matrix:''
     2    /12x,''a'',7x,''b'',7x,''c'',3x,''alpha'',4x,''beta'',3x,
     3    ''gamma''/'' cell'',6f8.2)') cellub
          go to 2005
        endif
 2033   write(*,'(//'' Cell dimension data OK (Y,N)? '')')
        call fread(010120,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .ne. 'Y') go to 2005
      endif
C
C---- Input of Density Measurement
C
 2035 write(*,'(//
     1  '' ***** Density measurement information'')')
      write(*,'('' Have you measured the density (Y/N)?'')')
      mdensi = .false.
      call fread(050020,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') then
        write(*,'(
     1'' Enter density, e.s.d. and temperature(K) of measurement:'')')
        call fread(050040,'RRRNNNNNNN')
        densit = realin(1)
        edensi = realin(2)
        tdensi = realin(3)
        mdensi = .true.
      elseif (inchar .ne. 'N') then
        go to  2035
      endif
C---- End of density input
 2040 if (mdensi) then
        write(*,'(/9x,''dobs'',5x,''esd'',4x,''temp''/
     1  5x,2f8.3,f8.1)') densit, edensi, tdensi
      else
        write(*,'('' No density measurements'')')
      endif
      write(*,'(//'' Density data OK (Y,N)? '')')
      call fread(050060,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2035
      if (inchar .ne. 'Y') go to 2040
C
C**** Input the physical dimensions of the crystal (faces or sphere)
C
 2105 write(*,'(//
     1  '' ***** Input of the physical dimensions of the crystal'')')
      write(*,'(''    N means no measurements of crystal shape '')')
      write(*,'(''    S means crystal is spherical in shape '')')
      write(*,'(''    F means polyhedral shape '')')
      write(*,'('' Make a choice for the shape of the crystal: '')')
      call fread(100020,'XNNNNNNNNN')
      crytyp = charin(1)(1:1)
      if (crytyp .eq. 'N') then
C---- No crystal shape measurements
        crytyp = 'N'
        write(*,'(//'' Physical dimensions of crystal not measured.'')')
      elseif (crytyp .eq. 'S') then
C---- Spherical shape
        crytyp = 'S'
        write(*,'('' Give sphere radius and e.s.d in mm:'')')
        call fread(100040,'RRNNNNNNNN')
        dface(1)  = realin(1)
        deface(1) = realin(2)
        write(*,'(//'' Sphere radius '',f6.3,''mm.   e.s.d. '',
     1  f6.3,''mm'')') dface(1), deface(1)
C---- Polyhedral shape
      elseif (crytyp .eq. 'F') then
        crytyp = 'F'
 2125   write(*,
     1    '('' How many indexed faces does the crystal have? '')')
        call fread(100055,'INNNNNNNNN')
        if (intin(1) .le. 0) go to 2125
        nface = intin(1)
        write(*,
     1  '('' Input one face per line: h, k, l, d, esd(d) in mm. '')')
        do 2130 i = 1, nface
          call fread(100060,'IIIRRNNNNN')
          hface(i) = intin(1)
          kface(i) = intin(2)
          lface(i) = intin(3)
          dface(i) = realin(4)
          deface(i) = realin(5)
 2130   continue
        write (*,'(//'' Crystal faces:''/
     1  ''  Num   h   k   l     d     esd(d)'')')
        write (*,'(1x,4i4,2f8.3)')
     1  (i, hface(i),kface(i),lface(i),dface(i),deface(i),i=1,nface)
      else
C---- Wrong input
        go to 2105
      endif
 2190 write(*,'(//'' Crystal shape OK (Y,N)? '')')
      call fread(100080,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2105
      if (inchar .ne. 'Y') go to 2190
C
C---- Chemical formula
C
 2205 write(*,'(//
     1 '' ***** Input of the chemical formula information'')')
      natmtp = 0
      do 2215 i = 1, 24
        peratm(i) = 0.0
        eperat(i) = 0.0
 2215 continue
      do 2210 i = 1, 5
 2210 lform(i) = .false.
C---- Z
      write(*,'('' Enter Z, the number of formula units per cell:'')')
      call fread(150020,'INNNNNNNNNN')
      zcell = intin(1)
C---- Assumed formula
      write(*,'('' Enter chemical formula  '',/,
     1'' all numbers are interpreted as subscripts.'',
     2'' separate element symbols (with ''/'' their coefficients '',
     3''and parentheses) from each other by a blank. '',/,
     4'' e.g.: (C6 H5) C H3  or  (C5 H5) Fe (C O)2  :'',/,
     5'' maximum 2 lines of 70 chars each, '',
     6''end each line with <CR> : '')')
      call fread(150040,'LNNNNNNNNN')
      chemfm(1:70)=line(1:70)
      call fread(150060,'BNNNNNNNNN')
      chemfm(71:140)=line(1:70)
C---- Convert formula to tabular sum form
      call chform
      if (natmtp .le. 0) then
        if (natmtp.eq.0)  write(*,'('' No known elements in formula'')')
        if (natmtp.eq.-1) write(*,'('' Position error in formula'')')
        if (natmtp.eq.-2) write(*,'('' Wrong special character  '')')
        if (natmtp.eq.-3) write(*,'('' Unknown element symbol   '')')
        if (natmtp.eq.-4) write(*,'('' Unbalanced parentheses   '')')
        if (natmtp.eq.-5) write(*,'('' Conversion error         '')')
        if (natmtp.eq.-6) write(*,'('' Overflow                 '')')
        go to 2205
      endif
C---- Chemical analysis
 2225 write(*,'('' Do you have a chemical analysis (Y/N)? '')')
      call fread(150080,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') then
        lform(1) = .true.
        write(*,'('' Enter observed relative weight % of each'',
     1  '' element after the element symbol as:''/
     2  '' observed % composition, e.s.d.:'')')
        do 2240 i = 1, natmtp
          write(*,'('' Wt % & esd of '',a2,'' : '')') namtyp(i)
          call fread(150100,'RRNNNNNNNN')
          peratm(i) = realin(1)
          eperat(i) = realin(2)
 2240   continue
      elseif (inchar .ne. 'N') then
        go to 2225
      endif
C---- Formula weight
 2250 write(*,'('' Do you have an observed value of the formula '',
     1 '' weight (Y/N)? '')')
      call fread(150120,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') then
        lform(2) = .true.
        write(*,'(
     1  '' Enter observed formula weight and e.s.d.(Daltons):'')')
        call fread(150140,'RRNNNNNNNN')
        forwg  = realin(1)
        eforwg = realin(2)
      elseif (inchar .ne. 'N') then
        go to 2250
      endif
C---- End of chemical formula
 2290 write(*,'(//'' Formula information:'',20x,
     1 ''Z number of formula units per cell '',i4//
     2 '' Chemical formula'',/,3x,a70,/,3x,a70,//
     3 5x,''el'',5x,''num'',5x,''wt%'',5x,''esd''/
     4 (5x,a2,f8.3,2f8.2))')
     5 zcell, chemfm(1:70),chemfm(71:140),
     6 (namtyp(i),fnamtp(i),peratm(i),eperat(i), i = 1 ,natmtp)
      if (lform(2)) write(*,'(//'' Formula weight & esd:'',2f10.3)')
     1 forwg, eforwg
      write(*,'(//'' Chemical formula data OK (Y,N)? '')')
      call fread(150160,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2205
      if (inchar .ne. 'Y') go to 2290
C
C---- Input of chemical compound information
C
 2750 write(*,'(//
     1  '' ***** Chemical compound information'')')
      write(*,'('' Enter the melting-point temperature of '',
     1  ''the compound in K:'')')
      call fread(200020,'RNNNNNNNNN')
      mptemp=realin(1)
      write(*,'('' Enter a description of the source of the '',
     1  ''chemical compound under study.''/
     2  '' If not known leave blank:'')')
      call fread(200040,'BNNNNNNNNN')
      cpdsou=line(1:76)
      write(*,'('' Enter the trivial name of the compound.'',
     1  '' If not known leave blank:'')')
      call fread(200060,'BNNNNNNNNN')
      namecm=line(1:76)
      write(*,'('' Enter the mineral name of the compound'',
     1  '' (if any - otherwise leave blank):'')')
      call fread(200080,'BNNNNNNNNN')
      nameml=line(1:76)
C---- End of name input
 2755 write(*,'(/9x,''Melting-point temperature:'',5x,f8.2,'' K'')')
     1 mptemp
      write(*,'( 9x,''Source of compound:''/1x,a76)') cpdsou
      write(*,'( 9x,''Trivial name of compound:''/1x,a76)') namecm
      write(*,'( 9x,''Mineral name of compound:''/1x,a76)') nameml
      write(*,'(//'' Chemical compound data OK (Y,N)? '')')
      call fread(200100,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2750
      if (inchar .ne. 'Y') go to 2755
C
C---- Input of crystal and crystal-related information
C
 2760 write(*,'(//
     1  '' ***** Crystal information'')')
      write(*,'('' Enter the colour of the crystal:'')')
      call fread(250020,'BNNNNNNNNN')
      cryscl=line(1:40)
      write(*,'('' Enter a description of the crystal '',
     1  ''quality and habit if known, otherwise blank:'')')
      call fread(250040,'BNNNNNNNNN')
      crysds=line(1:76)
      write(*,'('' Enter special details concerning the preparation'',
     1  '' or mounting of the crystal,''/'' otherwise blank:'')')
      call fread(250060,'BNNNNNNNNN')
      specdt=line(1:76)
C---- End of crystal input
 2765 write(*,'(/9x,''Crystal colour:'',5x,a40)') cryscl
      write(*,'( 9x,''Quality and habit of crystal:''/1x,a76)') crysds
      write(*,'( 9x,''Special preparation details:''/1x,a76)') specdt
      write(*,'(//'' Crystal data OK (Y,N)? '')')
      call fread(250080,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2760
      if (inchar .ne. 'Y') go to 2765
C
C---- Diffractometer attenuator information
C
C---- Have filters been used for which transmission factors were not given?
      ask(6) = .false.
      do 2410 i = 1, nfilt
        if (attfac(i) .lt. 1.0) ask(6) = .true.
 2410 continue
C---- Reciprocal attenuation factors and e.s.d.'s
      if (ask(6)) then
 2305   write(*,'(//
     1   '' ***** Attenuator filter information''/
     2   ''   Enter attenuation factors for the filter index numbers '',
     3   ''indicated as:'',/
     4   ''   Reciprocal attenuation factor (>1.0) and its e.s.d.'')')
        do 2321 i = 2, nfilt
 2320     write(*,'('' Filter index '',i3,'':'')') nfltcd(i)
          call fread(300020,'RRNNNNNNNN')
	  if(realin(1) .ge. 1.0) then
            attfac(i) = realin(1)
            ettfac(i) = realin(2)
          else
            Write(*,'('' Reciprocal attenuation factor should be '',
     1        '' greater than 1.0.  Reinput last entry.'')')
            go to 2320
          endif
 2321   continue
C---- End of diffractometer attenuator information
 2325   write(*,'(//'' Diffractometer attenuator information:'',//,
     2   '' attenuators: code'',4x,''factor'',4x,''e.s.d.''/
     3   (13x,i5,2f10.6))')
     5   (nfltcd(i), attfac(i), ettfac(i), i = 1, nfilt)
        write(*,'(//'' Diffractometer attenuator data OK (Y,N)? '')')
        call fread(300040,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'N') go to 2305
        if (inchar .ne. 'Y') go to 2325
      endif
C
C---- Diffractometer monochromator information
C
C---- Beta-filter or monochromator information
 2330 write(*,'(//
     1  '' ***** Monochromator information'')')
      write(*,'('' Did you use a beta-filter(B) or a monochromator(M)'',
     1 '' (B,M)? '')')
      call fread(350020,'XNNNNNNNNN')
      monoch = charin(1)(1:1)
      if (monoch .eq. 'M') then
C---- Monochromator
C---- Polarization ratio
 2340   write(*,'('' Have you a measured value '',
     1  '' of the monochromator polarization ratio (Y,N)? '')')
        call fread(350040,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'Y') then
C---- Measured values of polarization ratio
          write(*,'('' Give polarization ratio and e.s.d.:'')')
          call fread(350060,'RRNNNNNNNN')
          ratpol = realin(1)
          ertpol = realin(2)
          fmn2th = rad*acos(sqrt(abs(ratpol)))
          cmix = 0.0
        elseif (inchar .eq. 'N') then
C---- Values chosen by the programmme
C---- 2-theta value
          write(*,'('' Give monochromator 2-theta angle in degrees:'')')
          call fread (350080,'RNNNNNNNNN')
          fmn2th = realin(1)
          cmn2th = abs(cos(dar*fmn2th))
          cmix = 0.5
C---- Give an estimate based on c(mixing of K's) = 0.5 with e.s.d.(c) 0.2
          ratpol = cmix*cmn2th*cmn2th + cmix*cmn2th
          ertpol = 0.2*(cmn2th - cmn2th*cmn2th)
        else
C---- Wrong input
          go to 2340
        endif
 2350   write(*,'('' Have you a measured value of the monochromator'',
     1  '' inhomogeneity half-width (Y,N)? '')')
        call fread(350100,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'Y') then
C---- Measured values of inhomogeneity half-width
          write(*,'('' Give inhomogeneity half-width(mm):'')')
          call fread(350120,'RNNNNNNNNN')
          fmhomw = realin(1)
        elseif (inchar .eq. 'N') then
          fmhomw = 0.0
        else
          go to 2350
        endif
      elseif (monoch .eq. 'B') then
C----- Beta-filter
        write(*,'('' Give beta-filter absorption edge wavelength(A):''
     1  )')
        call fread (350140,'RNNNNNNNNN')
        bedwv = realin(1)
        ratpol = 1.0
        ertpol = 0.0
        fmn2th = 0.0
        cmix = 0.0
      else
C---- Wrong input
        go to 2330
      endif
C---- End of diffractometer monochromator information
 2390 write(*,'(//'' Diffractometer monochromator information:'',//,
     1 '' monochromator:'',4x,''polarz'',7x,''esd'',5x,
     2 ''dihed'',5x,''width''/15x,2f10.4,f10.1,f10.3//7x,''abs.e'',
     3 /1x,a1,f10.5)')
     4 ratpol, ertpol, fmnhed, fmhomw,
     5 monoch, bedwv
      write(*,'(//'' Diffractometer monochromator data OK (Y,N)? '')')
      call fread(350160,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2330
      if (inchar .ne. 'Y') go to 2390
C
C---- Counter dead time information
C
      if (ask(3)) then
 2405   write(*,'(//
     1  '' ***** Input of counter dead time information'')')
 2460   write(*,'('' Have you a measured value of the counter'',
     1  '' dead time (Y,N)? '')')
        call fread(400020,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'Y') then
C---- Measured values of counter dead-time and esd
          write(*,'('' Give measured dead-time and '',
     1    ''e.s.d. in microsecs:'')')
          call fread(400040,'RRNNNNNNNN')
          fdeadt = realin(1)
          efdead = realin(2)
        elseif (inchar .eq. 'N') then
          fdeadt = 0.0
          efdead = 0.0
        else
          go to 2460
        endif
C---- End of counter dead time information
 2490   write(*,'(//'' Diffractometer dead time information:'',//,
     1  5x,''Dead T'',7x,''Esd''/1x,2f10.2)') fdeadt, efdead
        write(*,'(//'' Counter dead time data OK (Y,N)? '')')
        call fread(400060,'XNNNNNNNNN')
        inchar = charin(1)(1:1)
        if (inchar .eq. 'N') go to 2405
        if (inchar .ne. 'Y') go to 2490
      endif
C
C---- Environmental Conditions information
C
 2470 write(*,'(//
     1 '' ***** Environmental conditions input'')')
      write(*,'('' Enter temperature (k) and pressure during '',
     1  ''intensity measurement:'')')
      call fread(450020,'RRNNNNNNNN')
      temint = realin(1)
      presnt = realin(2)
C---- End of environmental conditions information
 2475 write(*,'(//'' Environmental conditions information:'',//,
     1 '' Temperature(K)  Pressure(Units)''/1x,f14.3,f16.3)')
     4 temint, presnt
      write(*,'(//'' Environmental conditions data OK (Y,N)? '')')
      call fread(450040,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2470
      if (inchar .ne. 'Y') go to 2475
C
C---- Input of diffractometer operator's name, address etc
C
      open( unit=addrfl, file = addnam(1:trim(addnam)),
     1 status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     2 )
CVAX 2 ,readonly)
CUNX 2 )
 2710 WRITE(*,'(//
     1  '' ***** Personal information'')')
C---- Loop over names only
      rewind (addrfl)
      do 2718 i = 1, 200
        read(addrfl,'(A30)', end = 2719) line(1:30)
        write(*,'('' Name '',i2,'': '',a30)') i, line(1:30)
        do 2717 k = 1, 9
          read(addrfl,'(a1)', end = 2719)line(1:1)
 2717   continue
 2718 continue
 2719 continue
      Write(*,'('' name  0: not in list - enter manually''/
     1'' Enter the number of your name in the list:'')')
      call fread(500020,'INNNNNNNNN')
      if (intin(1).gt.0) then
C---- Fetch requested name and address
        rewind (addrfl)
        do 2727 k = 1, 10*(intin(1)-1)
          read(addrfl,'(a1)', end = 2729)line(1:1)
 2727   continue
        read(addrfl,'(A30)', end = 2729) family
        read(addrfl,'(a)', end = 2729) addres(1)
        read(addrfl,'(a)', end = 2729) addres(2)
        read(addrfl,'(a)', end = 2729) addres(3)
        read(addrfl,'(a)', end = 2729) addres(4)
        read(addrfl,'(a)', end = 2729) addres(5)
        read(addrfl,'(a)', end = 2729) email
        read(addrfl,'(a)', end = 2729) fax
        read(addrfl,'(a)', end = 2729) phone
C       Blank separating line
        read(addrfl,'(a1)', end = 2729)line(1:1)
 2729   continue
      else
        write(*,'('' Enter your name as: family name(s) , '',
     1    ''first name(s) or initials:'')')
        call fread(500040,'LNNNNNNNNN')
        family=line(1:30)
        write(*,'('' Enter your address (5 lines):'')')
        do 2720 i = 1,5
          call fread(500060,'BNNNNNNNNN')
          addres(i) = line(1:40)
 2720   continue
        write(*,'('' Enter your e-mail address:'')')
        call fread(500080,'BNNNNNNNNN')
        email=line(1:40)
        write(*,'('' Enter your telefax number:'')')
        call fread(500100,'BNNNNNNNNN')
        fax=line(1:40)
        write(*,'('' Enter your telephone number:'')')
        call fread(500120,'BNNNNNNNNN')
        phone=line(1:40)
      endif
C---- End of name/address input
 2715 write(*,'(/9x,''Name:'',t20,a30)') family
      write(*,'(9x,''Address:'',t20,a40,4(/t20,a40))') (addres(i),i=1,5)
      write(*,'(9x,''Email:'',t20,a40)') email
      write(*,'(9x,''Telefax:'',t20,a40)') fax
      write(*,'(9x,''Telephone:'',t20,a40)') phone
      write(*,'(//'' Name address data OK (Y,N)? '')')
      call fread(500140,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'N') go to 2710
      if (inchar .ne. 'Y') go to 2715
C---- Close file
      close(unit=addrfl, status = 'KEEP')
C
      return
      end
C
      subroutine fread(ncode,formt)
C
C**** Free format read for interactive input
C
      character formt*10, symb(0:13)*1
      integer   i,j,nval,mode,nchar,jchar,idecm
      integer ncode, trim
      logical   decm, minus, first
      include   'common.tlb'
C---- Formt contains the format of the line. I means an integer, R a real
C     and N nothing in the order of input.  All values are stored in
C     order in Intin and Realin. C and X mean character data separated by
C     blanks and stored in charin. For X the characters are converted to 
C     upper case and C they are left alone. L means a whole line of character
C     data to be found in line which must not be blank. B means a whole
C     line of character to be found in line which may be blank
      data  symb/
     1 '0','1','2','3','4','5','6','7','8','9','+','-','.',' '/
C       0   1   2   3   4   5   6   7   8   9   10 11  12  13
C
C---- Get the line
 95   line = ' '
      do 100 i = 1, 10
        intin(i) = 99999
        realin(i) = 99999.0
        charin(i) = '999999999999999'
 100  continue
C---- Defaults from saved answer file
      if (lquan) then
 110    if (nqaocd.ge.ncode) go to 130 
        read(qaofl,'(i6,a80)',end=120) nqaocd, qline
        go to 110
 120    nqaocd=999999
 130    if (nqaocd.eq.ncode) then
C---- A saved answer for this question
          write(*,'('' ** Saved answer - Return to accept: '',a)') 
     1     qline(1:trim(qline))
          if (hspeed) then
            line=qline        
            read(qaofl,'(i6,a80)',end=140) nqaocd, qline
            go to 145
 140        nqaocd=999999
 145        continue
          else
            read(*,'(a80)',end=9999) line
            if (line(1:10).eq.'          ') then
              line=qline        
              read(qaofl,'(i6,a80)',end=150) nqaocd, qline
            go to 155
 150        nqaocd=999999
 155        continue
            endif
          endif
        else
          read(*,'(a80)',end=9999) line
        endif
      else
C---- No saved answer file
        read(*,'(a80)',end=9999) line
      endif
      nval = 0
C---- mode 0 = looking for start, mode = 1 integer read, mode = 2 real read
C     mode 3 = character read, mode = 4 non-blank line read, mode = 5 line
C     (may be blank read), mode = 6 is character read converted to upper
C     case
      mode = 0
      if (formt(1:1) .eq. 'B') go to 9999
      do 8675 nchar = 1 ,80
C---- Check the input character
        jchar = -1
        first = .false.
        do 8605 i = 0, 13
          if (line(nchar:nchar) .eq. symb(i)) jchar = i
 8605   continue
C
C**** Treat mode = 0
C
        if (mode .eq. 0) then
C----   Looking for start
          if( jchar .ne. 13 ) then
C----     Got a start
            nval = nval + 1
            if (formt(nval:nval) .eq. 'I') mode = 1
            if (formt(nval:nval) .eq. 'R') mode = 2
            if (formt(nval:nval) .eq. 'C') mode = 3
            if (formt(nval:nval) .eq. 'L') mode = 4
            if (formt(nval:nval) .eq. 'X') mode = 6
            intin(nval) = 0
            charin(nval) = ' '
            ipchar = 0
            idecm = 0
            decm  = .false.
            minus = .false.
            first = .true.
          endif
        endif
C
C**** Treat mode = 1  - Integer
C
        if ( mode .eq. 1) then
          if (jchar. lt. 0) then
            go to 9950
          elseif (jchar .le. 9) then
C----     A digit
            intin(nval) = 10*intin(nval) + jchar
          elseif (jchar .eq. 10) then
C----     A +
            if ( .not. first) go to 9950
          elseif (jchar .eq. 11) then
C----     A -
            if ( .not. first) go to 9950
            minus = .true.
          elseif (jchar .eq. 12) then
C----     A .
            go to 9950
          elseif (jchar .eq. 13) then
C----     A space
            if (minus) intin(nval) = -intin(nval)
            mode = 0
          endif
        endif
C
C****  Treat mode = 2  - Real
C
        if (mode .eq. 2) then
          if (jchar. lt. 0) then
            go to 9950
          elseif (jchar .le. 9) then
C----     A digit
            intin(nval) = 10*intin(nval) + jchar
            if ( decm ) idecm = idecm + 1
          elseif (jchar .eq. 10) then
C----     A +
            if ( .not. first) go to 9950
          elseif (jchar .eq. 11) then
C----     A -
            if ( .not. first) go to 9950
            minus = .true.
          elseif (jchar .eq. 12) then
C----     A .
            if (decm) go to 9950
            decm = .true.
          elseif (jchar .eq. 13) then
C----     A space
            if (minus) intin(nval) = -intin(nval)
            realin(nval) = real(intin(nval))/10.0**idecm
            intin(nval) = 99999
            mode = 0
          endif
        endif
C
C**** Treat mode =  3 or 6 (without or with conversion to upper case)
C
        if ((mode .eq. 3).or.(mode .eq. 6)) then
          if (jchar .ne. 13) then
            ipchar = ipchar + 1
            if (ipchar .gt. 10) go to 9950
            inchar = line(nchar:nchar)
C----     Change to capitals
            if (mode .eq. 6) then
              do 8670 i = 1, 26
                if ( inchar .eq. alphl(i)) inchar = alphc(i)
 8670         continue
            endif
            charin(nval)(ipchar:ipchar) = inchar
          else
            mode = 0
          endif
        endif
C
C**** Treat mode = 4
C
        if (mode .eq. 4) then
C---- No need to do anything for mode 4. Whole line is the input
        endif
C
C**** Treat mode = 5
C
        if (mode .eq. 5) then
C---- No need to do anything for mode 5. Whole line is the input
        endif
 8675 continue
C---- Did we get everything we wanted
 8680 do 8685 i = 1 ,10
        if (formt(i:i) .eq. 'N') go to 8695
        j = i
 8685 continue
 8695 if (j .eq. nval) go to 9999
C---- Incorrect Input
 9950 write(*,'('' **** Wrong input - Try again ****'')')
      go to 95
C---- Copy input line to output
 9999 continue
      if(ncode.gt.0) 
     1  write(qanfl,'(i6,a)') ncode, line(1:trim(line))
C---- Subroutine exit
      return
      end

C
      integer function trim(line)
      character*(*) line
      do 100 trim = len(line),1,-1
         if (line(trim:trim) .ne. ' ') go to 200
 100  continue
 200  if (trim.le.0) trim=1
      return
      end
C
      subroutine chform
C
C===== Conversion of a string record of the chemical formula into the
C      component sums of elements. Adapted from a subroutine by B.Kilger
C===== 07/83 modified by R.Sievers 9.85 and 7.88
C
C    BEDEUTUNG DER ARGUMENTE:
C    CHEMFM:  INPUTSTRING (CHARACTER*(140)) MIT DER STRUKTURIERTEN
C             FORMEL
C    140:     LAENGE DES INPUTSTRINGS
C    namtyp:  CHARACTER*2-ARRAY (LAENGE: 24) MIT natmtp GEFUNDENEN
C             ELEMENT-SYMBOLEN
C    NELIST:  RUECKGABE DER ELEMENTNUMMERN, ALPHABETISCH
C    fnamtp:  RUECKGABE DER ENTSPRECHENDEN INDIZES
C    24:      MAXIMALLAENGE DER ARRAYS namtyp, NELIST UND fnamtp
C    natmtp:>0 - ZAHL DER GEFUNDENEN ELEMENTE
C           <0 - FEHLER. AM HAEUFIGSTEN: -1 : SYNTAX-FEHLER (Z.B. ZAHL
C                                             VOR ELEMENT STATT HINTER)
C                                        -2 : FALSCHES SONDERZEICHEN
C                                        -3 : FALSCHES ELEMENT
C                                        -4 : KLAMMERN NICHT AUSGEWOGEN
C
C
      include 'common.tlb'
C
      character*1 idmxx
      character*1 nb,flank,nc2,nz
      character*2 nel4
      character*2 nele(106)
C
      real f, ftot, fsav
      real fak(20)
C
      integer iendfl, ilpar, irpar, idig, iele, idezpt, istp, ipnt
      integer izifer, nv, nh, nr, i, jbeg, jl, ja, nesav, jend, klein
      integer ji, kleinp, j
      integer nelist(24)
C
C----VORGABE EINER ELEMENT-HIERARCHIE. DIE SPEZIFIZIERTEN
C     6 ELEMENTE WERDEN NACH VORN GEZOGEN, DIE UEBRIGEN FOLGEN
C     ALPHABETISCH.        Z. ZT.: C H O N S P
C
      integer ielord(6)
      data ielord /15,38,65,57,82,67/
C
      data nele /'AC','AG','AL','AM','AR','AS','AT','AU', 'B ','BA',
C                 1    2    3    4    5    6    7    8     9    10
     1 'BE','BI','BK','BR', 'C ','CA','CD','CE','CF','CL','CM','CO',
C       11   12   13   15    15   16   17   18   19   20   21   22
     2 'CR','CS','CU', 'D ','DY', 'ER','ES','EU', 'F ','FE','FM','FR',
C       23   24   25    26   27    28   29   30    31   32   33   34
     3 'GA','GD','GE', 'H ','HE','HF','HG','HO', 'I ','IN','IR',
C       35   36   37    38   39   40   41   42    43   44   45
     4 'K ','KR','KU', 'LA','LI','LU','LW', 'MD','MG','MN','MO',
C       46   47   48    49   50   51   52    53   54   55   56
     5 'N ','NA','NB','ND','NE','NI','NO','NP', 'O ','OS',
C       57   58   59   60   61   62   63   64    65   66
     6 'P ','PA','PB','PD','PM','PO','PR','PT','PU', 'RA','RB','RE',
C       67   68   69   70   71   72   73   74   75    76   77   78
     7 'RH','RN','RU', 'S ','SB','SC','SE','SI','SM','SN','SR','T',
C       79   80   81    82   83   84   85   86   87   88   89   90
     8 'TA','TB','TC','TE','TH','TI','TL','TM', 'U ', 'V ', 'W ',
C       91   92   93   94   95   96   97   98    99   100    101
     9 'XE', 'Y ','YB', 'ZN','ZR'/
C       102   103  104   105  106
C
      data nb/' '/
      data flank/' '/
      data nel4/' '/
      data nc2/' '/
      data idmxx/' '/
C
      natmtp = 0
      f = 1.0
      iendfl = 0
      ilpar = 0
      irpar = 0
      idig = 0
      iele = 0
      idezpt = 0
      istp = 0
      ipnt = 141
      nc2 = '0'
C
C            ***  LIES INPUTSTRING VON RECHTS NACH LINKS  ***
  100 ipnt = ipnt - 1
      if (ipnt.lt.1) goto 7000
      nb(1:1) = chemfm(ipnt:ipnt)
      nz(1:1) = chemfm(ipnt:ipnt)
C
C             A...Z     a...z     0...9    Blank  Sonderz.
C    -----------------------------------------------------
C    EBCDIC 193..233  129...169 240..249    64    65..128
C    ASCII   65...90   97...122  48...57    32    33...47
C
C----KLEIN-, GROSSBUCHSTABE ODER ZIFFER?
      if(lge(nz,'0')) go to 200
C
C----SONDERZEICHEN ?
      if(nz .eq.' ') go to 5000
      if(nz .eq.'(') go to 3000
      if(nz .eq.')') go to 4000
      if(nz .eq. '.') go to 4500
      go to 9902
C
C----KLEINBUCHSTABE, GROSSBUCHSTABE ODER ZIFFER
C     (GROSSBUCHSTABE ...)
 200  if(lle(nz,'Z') .and. lge(nz,'A')) go to 2000
      if(lle(nz,'9') .and. lge(nz,'0')) go to 1000
C     (KLEINBUCHSTABE ...)
      if(lgt(nz,'z') .or.  llt(nz,'a')) go to 9902
C    KLEIN->GROSS-KONVERSION
      if (lge(nb,'a') .and. lle(nb,'z')) then
         nb = char(ichar(nb) - ichar('a') + ichar('A'))
      end if
      go to 2000
C
C             ***  ZIFFER  ***
 1000 if (idig.gt.0) goto 1100
      if (iele.gt.0.or.ilpar.gt.0) goto 9901
      if (idezpt.le.0) f = 0.0
C
 1100 izifer = ichar(nz) - ichar('0')
      f = f + float(izifer)*10.**(idig-idezpt)
      idig = idig + 1
      irpar = 0
      goto 100
C
C             ***  BUCHSTABE  ***
 2000 if (ilpar.gt.0) goto 9901
      irpar = 0
      idig = 0
      idezpt = 0
      iele = iele + 1
      goto (2020,2030,9903),iele
 2020 nel4(1:1) = nb
      nel4(2:2) = flank
      goto 100
 2030 idmxx      = nel4(1:1)
      nel4(2:2) = idmxx
      nel4(1:1) = nb
      goto 100
C
C             ***  LINKE KLAMMER  ***
 3000 if (idig.gt.0.or.irpar.gt.0) goto 9901
      ilpar = 1
      if (iele.gt.0) goto 6000
 3010 istp = istp - 1
      if (istp.lt.0) goto 9904
      goto 100
C
C             ***  RECHTE KLAMMER  ***
 4000 if (iele.gt.0.or.ilpar.gt.0) goto 9901
      idig = 0
      idezpt = 0
      irpar = 1
      istp = istp + 1
      if (istp.gt.24) goto 9904
      fak(istp) = f
      f = 1.0
      goto 100
C
C             ***  DEZIMALPUNKT  ***
 4500 if (idezpt.gt.0) goto 9905
      f = f*10.**(-idig)
      idig = idig + 1
      idezpt = idig
      goto 100
C
C             ***  BLANK  ***
 5000 if (idig.gt.0) goto 9901
      ilpar = 0
      if (iele.lt.1) goto 100
C
C **********************************************************************
 6000 continue
C
      nv=1
      nh = 106
      nr=65
      go to 25
  10  nv=nr+1
  20  if (nv .gt. nh) go to 30
      nr=(nv+nh)/2
  25  if (nel4 .eq. nele(nr)) go to 40
      if (lgt(nel4,nele(nr))) go to 10
      nh=nr-1
      go to 20
C
  30  nr = 0
  40  continue
C **********************************************************************
C
      if (nr.le.0) goto 9903
      ftot = f
      if (istp.le.0) goto 6050
      do 6020 i=1,istp
 6020 ftot = ftot * fak(i)
 6050 if (natmtp.le.0) goto 6100
      do 6070 i=1,natmtp
      if (nr.ne.nelist(i)) goto 6070
      fnamtp(i) = fnamtp(i) + ftot
      goto 6150
 6070 continue
 6100 natmtp = natmtp + 1
      if(natmtp .gt. 24) go to 9906
      fnamtp(natmtp) = ftot
      nelist(natmtp) = nr
 6150 f = 1.0
      iele = 0
      if (ilpar.gt.0) goto 3010
      if (iendfl.le.0) goto 100
C
C            ***  ENDE DES SCANS VON RECHTS NACH LINKS  ***
 7000 iendfl = 1
      if (istp.ne.0) goto 9904
      if (iele.gt.0) goto 6000
C
C             ***  VERTAUSCHEN DER REIHENFOLGE  ***
      if (natmtp .le. 0) go to 9999
      jbeg = 1
C
C----DIE SPEZIFIZIERTEN ELEMENTE WERDEN NACH VORN GEZOGEN.
      do 7030 jl = 1,6
C----ELEMENT IELORD(JL) IN DER FORMEL SUCHEN:
      do 7010 ja = jbeg,natmtp
      if(nelist(ja) .eq. ielord(jl)) go to 7020
7010  continue
C----NICHT IN DER FORMEL.
      go to 7030
C
7020  nesav = nelist(ja)
      nelist(ja) = nelist(jbeg)
      nelist(jbeg) = nesav
      fsav = fnamtp(ja)
      fnamtp(ja) = fnamtp(jbeg)
      fnamtp(jbeg) = fsav
C
      jbeg = jbeg + 1
      if(jbeg .eq. natmtp) go to 9000
7030  continue
C
      jend = natmtp - 1
      do 7100 ja = jbeg,jend
      klein = 999
      do 7060 ji = ja,natmtp
      if (nelist(ji) .ge. klein) go to 7060
      klein = nelist(ji)
      kleinp = ji
7060  continue
      fsav = fnamtp(kleinp)
      fnamtp(kleinp) = fnamtp(ja)
      fnamtp(ja) = fsav
      nesav = nelist(kleinp)
      nelist(kleinp) = nelist(ja)
      nelist(ja) = nesav
7100  continue
C
9000  if (natmtp .gt. 0) then
         do 9020 j = 1,natmtp
            namtyp(j) = nele(nelist(j))
C
C              2. Buchstaben nach klein konvertieren
            if (namtyp(j)(2:2) .ne. ' ') then
               namtyp(j)(2:2) = char(ichar(namtyp(j)(2:2)) -
     1         ichar('A') + ichar('a'))
            end if
9020     continue
      end if
C
      go to 9999
C
C             ***  ALLG. STELLUNGSFEHLER  ***
 9901 natmtp = -1
      go to 9999
C             ***  FALSCHES SONDERZEICHEN  ***
 9902 natmtp = -2
      go to 9999
C             ***  FALSCHES ELEMENTSYMBOL  ***
 9903 natmtp = -3
      go to 9999
C             ***  KLAMMERN NICHT AUSGEWOGEN  ***
 9904 natmtp = -4
      go to 9999
C             ***  KONVERTIERFEHLER  ***
 9905 natmtp = -5
      go to 9999
C             ***  UEBERLAUF  ***
 9906 natmtp = -6
C
 9999 continue
      return
      end
C
      block data blockd
C
C**** COMMON and DATA statements for Diffractometer Conversion Progs
C
      include 'common.tlb'
C
C---- TYPE STATEMENTS
      integer i
C
C
C---- DATA STATEMENTS
C
      data digit/'0','1','2','3','4','5','6','7','8','9','-',' '/
C
      data alphc/'A','B','C','D','E','F','G','H','I','J','K','L','M',
     1           'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
      data alphl/'a','b','c','d','e','f','g','h','i','j','k','l','m',
     1           'n','o','p','q','r','s','t','u','v','w','x','y','z'/
C
      data includ /'INCLUDE '/
      data exclud /'EXCLUDE '/
C
      data dar, rad/ .017453293, 57.29577951/
C
      data nseq, nrefl/ 0, 0/
C
      data time, bastim, prvtim/ 0.0, 0.0, -99.0/
C
      data condon/ .false./
C
      data ask /16*.true./
C
      data meas(1)/'Omega_scans           '/
      data meas(2)/'Omega-twotheta_scans  '/
      data meas(3)/'Omega-X.X_theta_scans '/
      data meas(4)/'Variable_scan_type    '/
C
      data lquan  /.false./
      data hklopn /.false./
C
      data hmin, hmax, kmin, kmax, lmin, lmax
     1 /+9999,-9999,+9999,-9999,+9999,-9999/
C
      data stlmin, stlmax/50.0, 0.0/
C
      data (avint(i), sdint(i), count(i), i = 1,10)/30*0.0/
C
      data swdb /0.0/
      data sct  /0.0/
      data mxn  /0/
      data smo  /0/
      data bmo  /0/
      data bgt  /0.0/
      data parclp  /0/
      data numref  /0/
      data backsn  /0/
      data stdseq  /0/
      data fdeadt  /0.0/
      data efdead  /0.0/
      data numang  /0/
C
      data nfilt     /1/
      data nfltcd    /20*0/
      data attfac    /1.0,19*0.0/
      data ettfac(1) /0.0/
      data attn      /0.0/
C
      data nmcell   /0/
      data ctthn    /180./
      data ctthx    /0.0/
      data ctemp    /293.0/
C
      data navtab  /3/
      data wavtab  /
     1  1.54183, 1.540562,  1.544390, 0.6680, 0.3320,
     2  0.71073, 0.709300,  0.713590, 0.6671, 0.3329,
     3  0.56087, 0.5594075, 0.563798, 0.6671, 0.3329
     4 /
C
      data mstnn  /0/
C
      data sinter /0.0/
C
      end


------cut here------------------------------------------------------------------



C
C==== CIF Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992.
C
      subroutine conout
C
C**** Subroutine to output changes in condition during data collection
C
      include 'common.tlb'
C
C---- Close the Current reflection output section if open
C
      if (hklopn) then
        hklopn = .false.
      endif
C
      return
      end
C
      subroutine hedout
C
C**** Subroutine to output starting lines of the CIF file
C
      include 'common.tlb'
      integer im,id,iy
C
C==== Output the CIF data_ line
C
      write(scfsfl,'(/''data_ '',t6,a6)') compid
C
C==== Output the CIF audit items
C
      write(scfsfl,'(/)')
      call idate(im,id,iy)
      write(scfsfl,
     1  '(''_audit_creation_date '',t36,i2.2,''-'',i2.2,''-'',i2.2)')
     2  iy,im,id
      write(scfsfl,
     1  '(''_audit_creation_method '',t36,
     2  ''from_Flack_and_Blanc_DIFRAC '')')
C
      return
      end
C
      subroutine refout
C
C---- Output a reflection with the CIF diffrn items
C
      include 'common.tlb'
C
C---- Open up the CIF difrn loop
C
      if(.not. hklopn) then
        hklopn = .true.
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_diffrn_refln_index_h '')')
        write(scfsfl,'(''_diffrn_refln_index_k '')')
        write(scfsfl,'(''_diffrn_refln_index_l '')')
        write(scfsfl,'(''_diffrn_refln_counts_peak '')')
        write(scfsfl,'(''_diffrn_refln_counts_bg_1 '')')
        write(scfsfl,'(''_diffrn_refln_counts_bg_2 '')')
        write(scfsfl,'(''_diffrn_refln_angle_psi '')')
        write(scfsfl,'(''_diffrn_refln_elapsed_time '')')
        write(scfsfl,'(''_diffrn_refln_attenuator_code '')')
        write(scfsfl,'(''_diffrn_refln_scale_group_code '')')
        write(scfsfl,'(''_diffrn_refln_standard_code '')')
        write(scfsfl,'(''_diffrn_refln_scan_mode_backgd '')')
        write(scfsfl,'(''_diffrn_refln_fbcoeff_peak_z '')')
        write(scfsfl,'(''_diffrn_refln_fbcoeff_bg_peak_z '')')
        write(scfsfl,'(''_diffrn_refln_fbcoeff_bg_1_z '')')
        write(scfsfl,'(''_diffrn_refln_fbcoeff_bg_2_z '')')
        write(scfsfl,'(''_diffrn_refln_scan_width '')')
        write(scfsfl,'(''_diffrn_refln_scan_mode '')')
        write(scfsfl,
     1  '(''_diffrn_refln_detect_slit_horiz '')')
        write(scfsfl,
     1  '(''_diffrn_refln_detect_slit_vert '')')
      endif
C
      write(scfsfl,'(3i5,3f10.0,f9.3,f10.0,/,4i2,6f8.3,2f7.3)')
     2ih, ik, il, peak, backl, backr, psicry, time,
     3ptfilt, ptscl, ptstd, backsn, cnp, cnb, cnl, cnr,
     4swidth, scnrat, dethz, detvt
C
      return
      end
C
      subroutine basout
C
C**** Subroutine to output to the CIF all non-intensity information
C
      integer trim
      integer i,j,k
      integer ipnt
      character*20 valout
C
      include 'common.tlb'
C
C==== Output the CIF cell items
C
      write(scfsfl,'(/)')
      write(scfsfl,'(''_cell_length_a '',t36,a20)')
     1  valout(cell(1),ecell(1),3)
      write(scfsfl,'(''_cell_length_b '',t36,a20)')
     1  valout(cell(2),ecell(2),3)
      write(scfsfl,'(''_cell_length_c '',t36,a20)')
     1  valout(cell(3),ecell(3),3)
      write(scfsfl,'(''_cell_angle_alpha '',t36,a20)')
     1  valout(cell(4),ecell(4),3)
      write(scfsfl,'(''_cell_angle_beta '',t36,a20)')
     1  valout(cell(5),ecell(5),3)
      write(scfsfl,'(''_cell_angle_gamma '',t36,a20)')
     1  valout(cell(6),ecell(6),3)
      write(scfsfl,'(''_cell_formula_units_Z '',t36,a20)')
     1  valout(real(zcell),-1.0,1)
      write(scfsfl,'(''_cell_measurement_radiation '',t36,
     1 ''X-Rays_from_sealed_tube '')')
      write(scfsfl,'(''_cell_measurement_reflns_used '',t36,a20)')
     1  valout(real(nmcell),-1.0,1)
      if(ctemp.gt.0.0) write(scfsfl,
     1  '(''_cell_measurement_temperature '',t36,a20)')
     2  valout(ctemp,-0.1,1)
      if(ctthx.gt.0.0) write(scfsfl,
     1  '(''_cell_measurement_theta_max '',t36,a20)')
     2  valout(0.5*ctthx,-0.01,1)
      if(ctthn.gt.0.0) write(scfsfl,
     1  '(''_cell_measurement_theta_min '',t36,a20)')
     2  valout(0.5*ctthn,-0.01,1)
      write(scfsfl,'(''_cell_measurement_wavelength '',t36,a20)')
     1  valout(cwave,-0.00001,1)
      write(scfsfl,'(''_cell_volume '',t36,a20)')
     1  valout(volum,evolum,3)
C
C==== Output the CIF chemical items
C
      write(scfsfl,'(/)')
      if (lform(2)) write(scfsfl,'(
     1  ''_chemical_formula_weight_meas '',t36,a20)')
     2  valout(forwg,eforwg,3)
      if (mptemp.gt.1.0) write(scfsfl,'(''_chemical_melting_point '',
     1  t36,a20)') valout(mptemp,-0.01,1)
      k=trim(chemfm)
      if(k.gt.40) then
        write(scfsfl,'(''_chemical_formula_structural '',/,
     1  '''''''',a,'''''''')') chemfm(1:78)
      else if(k.gt.1) then
        write(scfsfl,'(''_chemical_formula_structural '',t36,
     1  '''''''',a,'''''''')') chemfm(1:k)
      endif
C---- Sum formula
      ipnt = -1
      line(1:80) = ' '
      do 2708 k = 1, natmtp
        ipnt = ipnt + 3
        line(ipnt-1:ipnt) = namtyp(k)(1:2)
        if (namtyp(k)(2:2).ne.' ') ipnt = ipnt + 1
        if (nint(fnamtp(k)).ne.1) line(ipnt:)=valout(fnamtp(k),-1.0,1)
        if(nint(fnamtp(k)).gt.9) ipnt = ipnt + 1
        if(nint(fnamtp(k)).gt.99) ipnt = ipnt + 1
 2708 continue
      write(scfsfl,'(''_chemical_formula_sum '',t36,
     1  '''''''',a,'''''''')') line(1:ipnt)
C---- Analytical quantities
      if(lform(1)) then
        ipnt = -1
        line(1:80) = ' '
        do 2709 k = 1, natmtp
          if(peratm(k).gt.0.0) then
            ipnt = ipnt + 3
            line(ipnt-1:ipnt) = namtyp(k)(1:2)
            if (namtyp(k)(2:2).ne.' ') ipnt = ipnt + 1
            line(ipnt:)=valout(peratm(k),eperat(k),1)
            ipnt = trim(line)
          endif
 2709   continue
        write(scfsfl,'(''_chemical_formula_analytical ''
     1    ,t36,'''''''',a,'''''''')') line(1:ipnt)
      endif
      k=trim(cpdsou)
      if(k.gt.40) then
        write(scfsfl,'(''_chemical_compound_source '',/,
     1  '''''''',a,'''''''')') cpdsou(1:k)
      else if(k.gt.1) then
        write(scfsfl,'(''_chemical_compound_source '',t36,
     1  '''''''',a,'''''''')') cpdsou(1:k)
      endif
      k=trim(namecm)
      if (k.gt.40) then
        write(scfsfl,'(''_chemical_name_common '',/,
     1  '''''''',a,'''''''')') namecm(1:k)
      else if (k.gt.1) then
        write(scfsfl,'(''_chemical_name_common '',t36,
     1  '''''''',a,'''''''')') namecm(1:k)
      endif
      k=trim(nameml)
      if (k.gt.40) then
        write(scfsfl,'(''_chemical_name_mineral '',/,
     1  '''''''',a,'''''''')') nameml(1:k)
      else if (k.gt.1) then
        write(scfsfl,'(''_chemical_name_mineral '',t36,
     1  '''''''',a,'''''''')') nameml(1:k)
      endif
C
C==== Output the CIF computing items
C
      write(scfsfl,'(/)')
      write(scfsfl,'(''_computing_data_collection '',t36,a40)') softwr
C
C==== Output the CIF diffrn items
C
      write(scfsfl,'(/)')
      write(scfsfl,'(''_diffrn_measurement_device '',
     1   t36,a40)') device
      if(mstnn.gt.0) then
        if(mstnn.eq.3) then
          write(meas(3)(7:9),'(f3.1)') scnrat
        endif
        write(scfsfl,'(''_diffrn_measurement_method '',
     1   t36,a22)')  meas(mstnn)
      endif
      if(fdeadt.gt.0.0) write(scfsfl,
     1  '(''_diffrn_radiation_detector_dtime '',t36,a20)')
     2  valout(fdeadt,efdead,3)
      write(scfsfl,'(''_diffrn_radiation_source '',t36,
     1 ''X-Rays_from_sealed_tube '')')
      IF(Monoch .EQ. 'b') WRITE(Scfsfl,'(
     1  ''_diffrn_radiation_filter_edge '',t36,a20)')
     2  valout(bedwv,-0.00001,1)
      write(scfsfl,'(''_diffrn_radiation_polarisn_ratio '',
     1  t36,a20)')  valout(ratpol,ertpol,3)
      write(scfsfl,'(''_diffrn_radiation_polarisn_norm '',
     1  t36,a20)')  valout(fmnhed,-0.01,1)
      if(fmhomw.gt.0.0) write(scfsfl,
     1  '(''_diffrn_radiation_inhomogeneity '',t36,a20)')
     2  valout(fmhomw,-0.0001,1)
      write(scfsfl,'(''_diffrn_ambient_temperature '',t36,a20)')
     1  valout(temint,-0.01,1)
      write(scfsfl,'(''_diffrn_reflns_theta_max '',t36,a20)')
     1  valout(rad*asin(stlmax*wavel),-0.01,1)
      write(scfsfl,'(''_diffrn_reflns_theta_min '',t36,a20)')
     1  valout(rad*asin(stlmin*wavel),-0.01,1)
      write(scfsfl,'(''_diffrn_reflns_limit_h_max '',t36,a20)')
     1  valout(real(hmax),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_limit_h_min '',t36,a20)')
     1  valout(real(hmin),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_limit_k_max '',t36,a20)')
     1  valout(real(kmax),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_limit_k_min '',t36,a20)')
     1  valout(real(kmin),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_limit_l_max '',t36,a20)')
     1  valout(real(lmax),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_limit_l_min '',t36,a20)')
     1  valout(real(lmin),-1.0,1)
      write(scfsfl,'(''_diffrn_reflns_number '',t36,a20)')
     1  valout(real(nrefl),-1.0,1)
      if (sinter.gt.0.0) write(scfsfl,'(
     1  ''_diffrn_standards_interval_time '',t36,a20)')
     2  valout(sinter/60.0,-1.0,1)
      write(scfsfl,'(''_diffrn_standards_number '',t36,a20)')
     1  valout(real(numref),-1.0,1)
C
C==== Output the CIF exptl items
C
      write(scfsfl,'(/)')
      k=trim(cryscl)
      if(k.gt.1) write(scfsfl,'(''_exptl_crystal_colour '',t36,
     1  '''''''',a,'''''''')') cryscl(1:k)
      if (mdensi) then
        write(scfsfl,'(''_exptl_crystal_density_meas '',
     1   t36,a20)') valout(densit, edensi,3)
        write(scfsfl,
     1   '(''_exptl_crystal_density_meas_temp '',
     2   t36,a20)') valout(tdensi,-0.01,1)
      endif
      k=trim(crysds)
      if (k.gt.40) then
        write(scfsfl,'(''_exptl_crystal_description '',/,
     1   '''''''',a,'''''''')') crysds(1:k)
      else if(k.gt.1) then
        write(scfsfl,'(''_exptl_crystal_description '',t36,
     1   '''''''',a,'''''''')') crysds(1:k)
      endif
      k=trim(specdt)
      if (k.gt.40) then
        write(scfsfl,'(''_exptl_special_details '',/,
     1   '''''''',a,'''''''')') specdt(1:k)
      else if (k.gt.1) then
        write(scfsfl,'(''_exptl_special_details '',t36,
     1   '''''''',a,'''''''')') specdt(1:k)
      endif
      if (crytyp .eq. 'S') then
        write(scfsfl,'(''_exptl_crystal_size_rad '',
     1   t36,a20)')  valout(dface(1), deface(1), 3)
      endif
C
C==== Output the CIF publ items to indicate persons name, address etc
C
      write(scfsfl,'(/)')
      write(scfsfl,'(''_publ_contact_author '',/,'';'',t36,a,
     1  5(/,t36,a),/,'';'')') family, (addres(i),i=1,5)
      k=trim(email)
      if(k.gt.1) write(scfsfl,'(
     1  ''_publ_contact_author_email '',t36,'''''''',a,'''''''')')
     2  email(1:k)
      k=trim(fax)
      if(k.gt.1) write(scfsfl,'(
     1  ''_publ_contact_author_fax '',t36,'''''''',a,'''''''')')
     2  fax(1:k)
      k = trim(phone)
      if(k.gt.1) write(scfsfl,'(
     1  ''_publ_contact_author_phone '',t36,'''''''',a,'''''''')')
     2  phone(1:k)
C
C==== Output the CIF atom type loop items
C
      if (natmtp .gt. 0) then
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_atom_type_symbol '')')
        write(scfsfl,'(''_atom_type_number_in_cell '')')
        if(lform(1)) write(scfsfl,'(''_atom_type_analytical_mass_% '')')
        do 2719 k = 1, natmtp
          if(lform(1)) then
            write(scfsfl,'(5x,a2,3x,f8.2,3x,a20)')
     1        namtyp(k), zcell*fnamtp(k), valout(peratm(k),eperat(k),1)
          else
            write(scfsfl,'(5x,a2,3x,f8.2)') namtyp(k), zcell*fnamtp(k)
          endif
 2719   continue
      endif
C
C==== Output the CIF  cell measurement loop items
C
      if (numang .gt. 0) then
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_cell_measurement_refln_index_h '')')
        write(scfsfl,'(''_cell_measurement_refln_index_k '')')
        write(scfsfl,'(''_cell_measurement_refln_index_l '')')
        write(scfsfl,'(
     1    ''_cell_measurement_refln_theta '')')
        write(scfsfl,'(
     1    ''_cell_measurement_refln_wl_id_z '')')
        do 2013 i = 1, numang
          write(scfsfl,'(3i5,1x,f10.4,i5)')
     1    hang(i), kang(i), lang(i), 0.5*twoang(i), wang(i)
 2013   continue
      endif
C
C==== Output the CIF diffrn attenuator loop items
C
      if (nfilt .gt. 0) then
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_diffrn_attenuator_code '')')
        write(scfsfl,'(''_diffrn_attenuator_scale '')')
        do 1030 i = 1 , nfilt
          write(scfsfl,'(i2,3x,a20)') nfltcd(i),
     1      valout(attfac(i),ettfac(i),3)
 1030   continue
      endif
C
C==== Output the CIF  diffrn standard loop items
C
      if (numref .gt. 0) then
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_diffrn_standard_refln_index_h '')')
        write(scfsfl,'(''_diffrn_standard_refln_index_k '')')
        write(scfsfl,'(''_diffrn_standard_refln_index_l '')')
        write(scfsfl,'(''_diffrn_standard_refln_code '')')
        do 100 j = 1, numref
          write(scfsfl,'(4i5)') href(j), kref(j), lref(j), j
 100    continue
      endif
C
C==== Output the CIF radiation wavelength loop items
C
      write(scfsfl,'(//''loop_ '')')
      write(scfsfl,'(''_diffrn_radiation_wavelength '')')
      write(scfsfl,'(''_diffrn_radiation_wavelength_wt '')')
      write(scfsfl,'(''_diffrn_radiation_wavelength_id '')')
      write(scfsfl,'(1x,f8.5,'' 1.00000  0'')') wavel
      write(scfsfl,'(1x,2f8.5,''  1'')') flam1, weigh1
      write(scfsfl,'(1x,2f8.5,''  2'')') flam2, weigh2
C
C==== Output the CIF  exptl crystal face loop items
C
      if (crytyp .eq. 'F') then
        write(scfsfl,'(//''loop_ '')')
        write(scfsfl,'(''_exptl_crystal_face_index_h '')')
        write(scfsfl,'(''_exptl_crystal_face_index_k '')')
        write(scfsfl,'(''_exptl_crystal_face_index_l '')')
        write(scfsfl,'(''_exptl_crystal_face_perp_dist '')')
        do 1025 i = 1 , nface
          write(scfsfl,'(3i5,3x,a20)')
     1      hface(i),kface(i),lface(i),valout(dface(i),deface(i),3)
 1025   continue
      endif
C
C==== End of CIF
C
      WRITE(scfsfl,'(//,
     1''#-end-end-end-end-end-end-end-end-end-end-end-end-end-end'')')
C
      return
      end
C
      subroutine set
C
C----- Set up fixed characteristics for CIF output
C
      include 'common.tlb'
C
      extn = '.cif'
C
      return
      end
C
      subroutine dtflnd
C============================================================================C
C      Routine called after the end of a reflection data file to finish off  C
C      output reflection file. This version is for CIF                       C
C============================================================================C
      include 'common.tlb'
      if (hklopn) then
        hklopn = .false.
      endif
C
      return
      end
C
      character*20 function valout (v1,v2,i3)
C
C==== A value (V1) and its esd (V2) are formatted into character form with
C     the necessary number of decimal digits and the esd in parentheses.
C     NUMESD controls the number of digits for the esd.
C     A negative value of V2 means no esd is available and V2 gives instead
C==== the required precision for printing the value.
C
      include 'common.tlb'
      integer   point, numesd, power, powcur, nzerop, valint, jk
      integer   esdint
      integer   i3
      real      v1, v2, value, esd
      logical   esdok, positf
C
C==== Initialisations
C
C     Set up the output string and the current position pointer.
      valout = '                    '
      point  = 1
C
C==== Keep the value and its esd.
C
      value  = v1
      esd    = v2
      numesd = i3
C
C==== Deal with negative VALUE
C
      if (value.lt.0.0) then
        positf = .false.
        value = abs(value)
      else
        positf = .true.
      endif
C
C==== Test to see if esd is OK. If not do fixes
C
      if (esd.gt.0.0) then
C---- Valid esd given
        esdok = .true.
      else if (esd.lt.0.0) then
C---- No esd for this value
        esdok = .false.
        esd = abs(esd)
        numesd = 1
      else
C---- Esd of 0.0 given print with 6 decimals
        esdok = .false.
        esd = 0.00001
      endif
C
C==== Set up the POWER to which VALUE must be raised in order to get
C     only the significant digits output.
C
      power = 0
C---- Make value of esd less than 0.4 in the loop
C**** Loop
 3140 if (esd.lt.0.4) go to 3145
        power = power - 1
        esd = esd/10.0
        go to 3140
 3145 continue
C**** Bottom of loop
C
C---- Now make rounded integer value of esd greater than NUMESD in the loop
      esdint = int(esd + 0.5)
C**** Loop
 3125 if (esdint.ge.numesd) go to 3135
        power = power + 1
        esd = esd*10.0
        esdint = int(esd + 0.5)
        go to 3125
 3135 continue
C**** Bottom of loop
C
C==== Clean up VALUE and take care of the special case of VALUE zero
C
C---- Save the power and number of non-significant zeroes
      powcur =  power
      nzerop = -power
C---- Raise VALUE to POWER, round and make it a "clean" real
      valint = nint((10.0**power)*value)
      value =  real(valint) + 0.5
C---- VALUE of zero has only one zero and is not negative
      if (valint.eq.0) then
        nzerop = min0(nzerop,0) + 1
        positf = .true.
      endif
C
C==== Make  VALUE to be fractional (i.e. less than 1.0)
C==== with POWCUR less or equal to zero
C
C**** Loop
 3170 if ((valint.le.0).and.(powcur.le.0)) go to 3180
        powcur = powcur - 1
        value = value/10.0
        valint = int(value)
      go to 3170
 3180 continue
C**** Bottom of Loop
C
C---- Output the sign
      if (.not. positf) then
        valout(point:point) = '-'
        point = point+1
      endif
C
C==== Output the significant digits of VALUE
C
      do 3185 jk = powcur, power-1
        if (jk.eq.0) then
C----   Put in a decimal point if needed here
          if (jk.eq.powcur) then
C----     Put in 0. instead .
            valout(point:point) = '0'
            point = point + 1
          endif
          valout(point:point) = '.'
          point = point + 1
        endif
C----   Get next digit and update parameter
        value = value*10.0
        valint = int(value)
        value = value - real(valint)
        valout(point:point) = digit(valint)
        point = point + 1
 3185 continue
C
C==== Output the non-significant digits of VALUE
C
      do 3225 jk = 1, nzerop
        valout(point:point) = '0'
        point = point + 1
 3225 continue
C
C==== Output the esd as (xxx) if present.
C
      if (esdok) then
C----   Output the opening parenthesis
        valout(point:point) = '('
        point = point + 1
C----   Save the power
        powcur =  power
C----   Make ESD clean as its already to correct power
        esd =  real(esdint) + 0.5
C----   Make ESD to be fractional (i.e. less than 1.0) in the loop
C**** Loop
 3270 if (esdint.le.0) go to 3280
        powcur = powcur - 1
        esd = esd/10.0
        esdint = int(esd)
      go to 3270
 3280 continue
C**** Bottom of Loop
C
C==== Output the significant digits of ESD
C
      do 3285 jk = powcur, power-1
C----   Get next digit and update parameter
        esd = esd*10.0
        esdint = int(esd)
        esd = esd - real(esdint)
        valout(point:point) = digit(esdint)
        point = point + 1
 3285 continue
C----   Output the non-significant zeroes
        do 3275 jk = 1, -power
          valout(point:point) = '0'
          point = point + 1
 3275   continue
C----   Output the closing parenthesis
        valout(point:point) = ')'
        point = point + 1
      endif
C
      return
      end



------cut here------------------------------------------------------------------


C
C==== SCFS-90 Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
C
      subroutine conout
C
C**** Subroutine to output changes in condition during data collection
C
      include 'common.tlb'
C
C---- Close the Current HKL section if open
C
      if (hklopn) then
C---- EOS for HKL INTEGrated section
        hklopn = .false.
        nseq = nseq + 1
        write (scfsfl,'(''*EOS'',71x,i5)') nseq
      endif
C
      return
      end
C
      subroutine hedout
C
C**** Subroutine to output SCFS header records
C
      include 'common.tlb'
C
C==== Output the TITLE section
C
      nseq = nseq + 1
      write(scfsfl, '(''TITLE'',70x,i5)') nseq
      nseq = nseq + 1
      write(scfsfl, '(''*'',a66,a6,2x,i5)') title, compid, nseq
C
      return
      end
C
      subroutine refout
C
C---- Output a reflection in the SCFS format
C
      include 'common.tlb'
C
C---- Open up the HKL 	INTEgrated section
C
      if(.not. hklopn) then
        hklopn = .true.
        nseq = nseq + 1
        write(scfsfl,'(75x,i5)') nseq
        nseq = nseq + 1
        write(scfsfl,'(''HKL INTEgrated'',61x,i5)') nseq
      endif
C
      nseq = nseq + 2
      write(scfsfl,'(''1INT '',a4,'' 2'',3i5,3f10.0,f9.3,f10.0,i5,/
     1''2INT '',4i2,6f8.3,2f7.3,i5)')
     2dsk, ih, ik, il, peak, backl, backr, psicry, time, nseq-1,
     3ptfilt, ptscl, ptstd, backsn, cnp, cnb, cnl, cnr,
     4swidth, scnrat, dethz, detvt, nseq
C
      return
      end
C
      subroutine basout
C
C**** Subroutine to output to the SCFS the non-intensity information
C
      integer i,j,k
C
      include 'common.tlb'
C
C==== Output the CELL DIMENSIONS  Section
C
      nseq = nseq + 1
      write(scfsfl,'(75x,i5)') nseq
      nseq = nseq + 1
      write(scfsfl,'(''CELL DIMENSION '',60X,i5)') nseq
C---- Cell dimensions
      nseq = nseq + 1
      write(scfsfl,'(1x,''CELL'',a4,1x,6f10.4,2i5)') dsk,
     1  (cell(i), i=1,6), zcell, nseq
C---- Errors in Cell dimensions
      if (ecell(1).gt.0.0) then
        nseq = nseq + 1
        write(scfsfl,'(1x,''ERRS'',a4,1x,6f10.4,5x,i5)') dsk,
     1   (ecell(i), i=1,6), nseq
      endif
C---- Volume of Cell
      nseq = nseq + 1
      write(scfsfl,'(1x,''VOL '',a4,1x,2f10.4,45x,i5)') dsk,
     1  volum, evolum, nseq
C---- PHYS information
      if (mdensi .or. lform(2)) then
        nseq = nseq + 1
        write(line,'(1x,''PHYS'',a4,66x,i5)') dsk, nseq
C---- Measured density
        if (mdensi) write(line(11:40),'(3f10.4)')
     1      densit, edensi, tdensi
C---- Measured Formula weight
        if (lform(2)) write(line(41:60),'(2f10.4)') forwg, eforwg
        write(scfsfl,'(a80)') line
      endif
C---- EOS of cell dimension section
      nseq = nseq + 1
      write(scfsfl,'(''*EOS'',71x,i5)') nseq
C
C==== Output the HKL / ANGL Section
C
      if (numang .gt. 0) then
        nseq = nseq + 1
        write(scfsfl,'(75x,i5)') nseq
        nseq = nseq + 1
        write(scfsfl,'(''HKL'',72X,i5)') nseq
        do 2013 i = 1, numang
C---- ANGL lines
          nseq = nseq + 1
          write(scfsfl,'(1x,''ANGL'',a4,3i5,i1,f10.4,40x,i5)')
     1      dsk, hang(i), kang(i), lang(i), wang(i), twoang(i),
     2      nseq
 2013   continue
C---- EOS of HKL/ ANGL section
        nseq = nseq + 1
        write(scfsfl,'(''*EOS'',71x,i5)') nseq
      endif
C
C==== Output the CONDITIONS Section
C
      nseq = nseq + 1
      write(scfsfl,'(75x,i5)') nseq
      nseq = nseq + 1
      write(scfsfl,'(''CONDITIONS'',65X,i5)') nseq
C
C---- Conditions for measuring cell dimensions
C
      if(nmcell.gt.0) then
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''CELL'',a4,''X'',f10.6,3f10.2,i10,15x,i5)') dsk,
     2    cwave, ctemp, ctthn, ctthx, nmcell, nseq
      endif
C
C---- Conditions used for measuring the intensities
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''INT '',a4,''X'',f10.6,f10.2,2f10.4,i10,15x,i5)') dsk,
     2  wavel, temint, stlmin, stlmax, nrefl, nseq
C
C---- Reflection ranges used in measuring intensities
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''HKL '',a4,a1,6i10,5x,i5)') dsk,
     2  difnam, hmin, hmax, kmin, kmax, lmin, lmax, nseq
C
C---- Details of standard reflections
C
      do 100 j = 1, numref
        if (count(j) .gt. 1.5) then
          avint(j) = avint(j)/count(j)
          sdint(j) = sqrt((sdint(j)-count(j)*avint(j)*avint(j))
     1             /(count(j)-1.0))
        else
          avint(j) = 0.0
          sdint(j) = 0.0
        endif
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''STD '',a4,1x,3i10,2f10.3,15x,i5)') dsk,
     2    href(j), kref(j), lref(j), avint(j), sdint(j), nseq
 100  continue
C
C---- Polyhedral crystal dimensions
C
      if (crytyp .eq. 'F') then
C---- Polyhedron
        do 1025 i = 1 , nface
          nseq = nseq + 1
          write(scfsfl,'(1x,''FACE'',a4,1x,f10.6,3i10,
     1      '' 999999999 999999999'',5x,i5)')
     2      dsk,dface(i),hface(i),kface(i),lface(i),nseq
          nseq = nseq + 1
          write(scfsfl,'(1x,''FARE'',a4,1x,f10.6,3i10,
     1      '' 999999999 999999999'',5x,i5)')
     2      dsk,deface(i),hface(i),kface(i),lface(i),nseq
 1025   continue
      endif
C
C---- Spherical crystal dimensions
C
      if (crytyp .eq. 'S') then
        nseq = nseq + 1
        write(scfsfl,'(1x,''SHAP'',a4,''S'',2f10.6,45x,i5)') dsk,
     1    dface(1), deface(1), nseq
       endif
C
C---- Dead Time Information
C
      if(fdeadt .ne. 0.0) then
        nseq = nseq + 1
        write(scfsfl,'(1x,''DETR'',a4,1x,2f10.2,45x,i5)')
     1  dsk, fdeadt, efdead, nseq
      endif
C
C---- Beta-filter Information
C
      if(monoch .eq. 'B') then
        nseq = nseq + 1
        write(scfsfl,'(1x,''BFIL'',a4,1x,f10.6,55x,i5)')
     1  dsk, bedwv, nseq
      endif
C
C---- Attenuator Factors
C
        do 1030 i = 1 , nfilt
          nseq = nseq + 1
          write(scfsfl,'(1x,''ATTF'',a4,1x,i10,2f10.6,35x,i5)')
     1      dsk,nfltcd(i),attfac(i),ettfac(i),nseq
 1030   continue
C
C---- Monochromator information
C
        nseq = nseq + 1
        write(scfsfl,'(1x,''POLO'',a4,1x,2f10.3,f10.1,f10.3,25x,i5)')
     1  dsk, ratpol, ertpol, fmnhed, fmhomw, nseq
C
C---- Wavelength of spectral line components
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''WAVL'',a4,1x,4f10.6,25x,i5)') dsk,
     2  flam1, weigh1, flam2, weigh2, nseq
C
C---- EOS for CONDITIONS
C
      nseq = nseq + 1
      write(scfsfl,'(''*EOS'',71x,i5)') nseq
C
C==== Output the COMMENTS Section
C
      nseq = nseq + 1
      write(scfsfl,'(75x,i5)') nseq
      nseq = nseq + 1
      write(scfsfl,'(''COMMENTS'',67X,i5)') nseq
C
C---- Experimental description of scan type
C
      if(mstnn.gt.0) then
        if(mstnn.eq.3) then
          write(meas(3)(7:9),'(f3.1)') scnrat
        endif
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''EXPT'',a4,1x,A22,43x,i5)') dsk, meas(mstnn), nseq
      endif
C
C---- Author's name
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''AUTH'',a4,1x,A30,35x,i5)') dsk, family, nseq
C
C---- Trivial name
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''NAME'',a4,1x,A65,i5)') dsk, namecm(1:65), nseq
C
C---- Source of compound
C
      nseq = nseq + 1
      write(scfsfl,
     1  '(1x,''SRCE'',a4,1x,A65,i5)') dsk, cpdsou(1:65), nseq
C
C---- Chemical formula of compound
C
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''FORM'',a4,1x,A65,i5)') dsk, chemfm(1:65), nseq
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''FORM'',a4,1x,A65,i5)') dsk, chemfm(66:130), nseq
        nseq = nseq + 1
        write(scfsfl,
     1    '(1x,''FORM'',a4,1x,A10,55X,i5)') dsk, chemfm(131:140), nseq
C
C---- EOS for COMMENTS
C
      nseq = nseq + 1
      write(scfsfl,'(''*EOS'',71x,i5)') nseq
C
C==== Output the FORMULA  Section
C
      nseq = nseq + 1
      write(scfsfl,'(75x,i5)') nseq
      nseq = nseq + 1
      write(scfsfl,'(''FORMULA'',68X,i5)') nseq
C---- Assumed formula
      j = 0
      k = -9
      line = ' '
      do 1005 i = 1 , natmtp
        j = j + 1
        k = k + 10
        write(line(k:k+9),'(a2,f8.4)') namtyp(j), fnamtp(j)
        if ( (k .eq. 51) .or. (i .eq. natmtp) ) then
          nseq = nseq + 1
          write(scfsfl,'(1x,''FORL'',a4,a60,6x,i5)')
     1     dsk,line(1:60), nseq
          k = -9
          line = ' '
        endif
 1005 continue
C---- Chemical analysis in wt %
      j = 0
      k = -9
      line  = ' '
      lin2 = ' '
      do 1010 i = 1 , natmtp
        j = j + 1
        if ( peratm(j) .gt. 0.0) then
          k = k + 10
          write(line(k:k+9),'(a2,f8.4)') namtyp(j), peratm(j)
          write(lin2(k:k+9),'(a2,f8.4)') namtyp(j), eperat(j)
        endif
        if ( (k.eq.51).or.((i.eq.natmtp).and.(k.gt.0)) ) then
          nseq = nseq + 1
          write(scfsfl,'(1x,''AWPC'',a4,a60,6x,i5)')
     1      dsk,line(1:60), nseq
          nseq = nseq + 1
          write(scfsfl,'(1x,''EWPC'',a4,a60,6x,i5)')
     1      dsk,lin2(1:60), nseq
          k = -9
          line = ' '
          lin2 = ' '
        endif
 1010 continue
C---- EOS of formula section
      nseq = nseq + 1
      write(scfsfl,'(''*EOS'',71x,i5)') nseq
C
C==== END section
C
      nseq = nseq + 1
      write (scfsfl, '(75x,i5)') nseq
      nseq = nseq + 1
      write (scfsfl, '(''END'',72x,i5)') nseq
C
      return
      end
C
      subroutine set
C
C----- Set up fixed characteristics for SCFS output
C
      include 'common.tlb'
C
      extn = '.sfi'
C
      return
      end
C
      subroutine dtflnd
C============================================================================C
C      Routine called after the end of a reflection data file to finish off  C
C      output reflection file.  This version is for SCFS                     C
C============================================================================C
      include 'common.tlb'
C---- EOS line to finish off the HKL INTEGERATED
      if (hklopn) then
        hklopn = .false.
        nseq = nseq + 1
        write (scfsfl,'(''*EOS'',71x,i5)') nseq
      endif
C
      return
      end



------cut here------------------------------------------------------------------


C
C==== CAD4(Version 4) Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
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
C
      include 'common.tlb'
C
C---- In the DIFRAC module for the CAD4 version 4.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
C---- Take care of complicated Enraf-Nonius logic for anti-reflections
C
      if (nfridl.lt.0) then
        ih = -ih
        ik = -ik
        il = -il
      endif
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
      elseif (code6(1:1) .eq. 'I') then
        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
      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) 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*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
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
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
      call refout
      return
      end
      subroutine dtflrd(jj)
C
C===== Subroutine to read a CAD4 diffractometer output file
C
      integer i,jj, mrefl, trim
      real sinsqa
C
      include 'common.tlb'
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
C---- Initialisations
C
      ptscl = 1
      firstr = .false.
      reoren = .false.
      mrefl = 0
C
C---- Get the data file name and open it
C
      open( unit = datafl, 
     1 file = compid(1:trim(compid))//datxtn,
     2 status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     3 )             
CVAX 3 ,readonly)
CUNX 3 )
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(1:5) .eq. '     ') then
C---- Do nothing - we need to ignore entirely blank lines
C
C---- Type 21 line.
C
      elseif (line(1:3) .eq. ' 21') then
        numref = 0
        parclp = 1
        read(line,'(45x,f4.0,3f4.1)')
     1    radius, slit, apmin, apmax
C
C---- Type 22 line.
C
      elseif (line(1:3) .eq. ' 22') then
        parclp = parclp + 2
        read(line,'(25x,2f5.2,f2.0)')
     1    apta, aptb, scnrat
C
C---- Type 23 line.
C
      elseif (line(1:3) .eq. ' 23') then
        parclp = parclp + 4
C
C---- Type 24 line.
C
      elseif (line(1:3) .eq. ' 24') then
        parclp = parclp + 8
C
C---- Type 25 line.
C
      elseif (line(1:3) .eq. ' 25') then
        parclp = parclp + 16
C
C---- Type 26 line.
C
      elseif (line(1:3) .eq. ' 26') then
        read(line,'(34x,2f10.7)')
     1    sinsqa, cosalp
        sinalp = sqrt(sinsqa)
        parclp = parclp + 32
C
C---- Type 0 line. Get the check reflection indices
C
      elseif (line(1:3) .eq. '   ') then
        if (((line(15:15).eq.'H').or.(line(15:15).eq.'N')).and.
     1     ((line(19:19).eq.'I').or.(line(19:19).eq.'F'))) then
          numref = numref + 1
          read(line,'(5X,3I3)')
     1      href(numref), kref(numref), lref(numref)
        endif
C
C---- Type 31 line
C
      elseif (line(1:3) .eq. ' 31') then
        parclp = parclp + 64
        read(line,'(5x,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(1:3) .eq. ' 32') then
        parclp = parclp + 128
        if (parclp .eq. 255) firstr = .true.
        if (parclp .eq. 447) then
          reoren = .true.
          parclp = 255
        endif
        read(line,'(5x,3f9.6,2x,3f9.6)')
     1   ub(3,1),ub(3,2),ub(3,3),flam1,flam2,beta
C
C---- Only allow the reflection lines if control data present
C
      elseif (parclp .eq. 255) then
C
C---- Type 1 line
C
        if (line(1:3) .eq. '  1') then
          read(line,'(3x,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
C
C---- Type 2 line
C
        elseif (line(1:3) .eq. '  2') then
          read (line,'(3x,i6,f8.3,3f9.3,f7.3,f7.0,i3)')
     1      i, theta, phi, omega, kappa, swidth, time, nfridl
          if ( i .eq. mrefl) then
            if (firstr .or. reoren) call initlz
            call refrd(0)
            mrefl = 0
          endif
        endif
      endif
      go to 1000
C
C---- Finish off this data file
C
 2000 continue
      close (unit = datafl, status = 'KEEP')
      return
      end
C
C===== Subroutine to read a CAD4 diffractometer parameter output file
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
C===== On the Enraf-Nonius CAD4, with this programme installed on the
C      diffractometer itself, this subroutine could be made to read the
C      <compid>.CRY and the GONCA0.DAT files for information on the machine.
C      This is not yet implemented in this version but will come soon.
C
      subroutine prflrd
C
      return
      end
      subroutine choos
C
C----- Set up fixed characteristics for CAD4 diffractometer
C
      include 'common.tlb'
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
C
      difnam = 'E'
      device = 'Enraf_Nonius_CAD4_4circle_diffractometer'
      softwr = 'Enraf_Nonius_CAD4_version_4.0           '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = 'pck$:[cristal]address.lst       '
CNORMAL addnam = 'address.lst                     '
      addnam = 'address.lst                     '
C
C---- File extension for data file name
      datxtn = '.dat'
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
C
C---- Options for Quesan
C
      return
      end
      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.tlb'
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
C---- Wavelength
        weigh1 = 0.6666667
        weigh2 = 0.3333333
        weigh3 = 0.0000000
        flam3  = 0.0
        wavel = weigh1*flam1 + weigh2*flam2
C---- Vertical detector opening
        detvt = rad*slit/radius
C
C---- Set Firstr to indicate that INITLZ does not need calling again.
C
        firstr = .false.
      endif
C---- Output the conditions
      call conout
C
      return
      end
      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, In the Press. 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.tlb'
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
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



------cut here------------------------------------------------------------------


C
C==== CAD4(Version 5) Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
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.tlb'
C---- Reflections with Goniometer status C have no intensity measurement
      if ( code6(2:2) .ne. 'C') 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*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.tlb'
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  )             
CVAX 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))//digit(i)//digit(j)//datxtn,
     2  status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     3  )          
CVAX 3  ,readonly)
CUNX 3  )
        write(*,'(/'' Reading file: '',A)')
     1   compid(1:trim(compid))//digit(i)//digit(j)//datxtn
      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.tlb'
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           '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = '[cad4.system]address.lst        '
CNORMAL addnam = 'address.lst                     '
        addnam = '[cad4.system]address.lst        '
C
C---- File extension for data file name
      datxtn = '.dat'
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.tlb'
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)
          endif
 201    continue
        cwave = flam1
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, In the Press. 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.tlb'
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
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.tlb'
      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
C
C---- The R format is only allowed in DATA statements, is very specific
C     to VAX FORTRAN, packs three characters into two bytes and is only
C---- because the Enraf-Nonius software uses it.  Why? Why? Why? Why?
C
      character*6 cption(8), chkl(6), cpsifl, cseq, cnfms
      real        rption(8)
      integer     rhkl(6)
      data rption(1), cption(1) /6rbisect, 'BISECT'/
      data rption(2), cption(2) /6razimut, 'AZIMUT'/
      data rption(3), cption(3) /6rflat  , 'FLAT  '/
      data rption(4), cption(4) /6rneedle, 'NEEDLE'/
      data rption(5), cption(5) /6rvechor, 'VECHOR'/
      data rption(6), cption(6) /6rzigzag, 'ZIGZAG'/
      data rption(7), cption(7) /6rsephkl, 'SEPHKL'/
      data rption(8), cption(8) /6rfile  , 'FILE  '/
      data rhkl(1),   chkl(1)   /3rhkl,    'HKL   '/
      data rhkl(2),   chkl(2)   /3rhlk,    'HLK   '/
      data rhkl(3),   chkl(3)   /3rkhl,    'KHL   '/
      data rhkl(4),   chkl(4)   /3rklh,    'KLH   '/
      data rhkl(5),   chkl(5)   /3rlhk,    'LHK   '/
      data rhkl(6),   chkl(6)   /3rlkh,    'LKH   '/
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.tlb'
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 VAX/VMS
C     used to run the diffractometer. This is because the .PAR and .CRY
C     files are very VAX/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



------cut here------------------------------------------------------------------


C
C==== Niv411 Module of DIFRAC by E. Blanc, University of Lausanne
C==== Version 7th October 1992
C
      subroutine refrd(itype)
C**** ==========================================================================
C     Routine to interpret reflection data from a NICOLET P21
C     All data is output by calling REFOUT
C**** ==========================================================================
      integer itype,i
      real pmt
      logical oka
C
      include 'common.tlb'
C---- Get theta
      theta = twoth/2.0
C
C---- In the DIFRAC module for the Siemens P21 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
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 (nrefl .gt. 0) then
        ptstd = 0
        stdseq = 0
      else
        do 10 i = 1, numref
          if ((ih .eq. href(i)) .and. (ik .eq. kref(i)) .and.
     $        (il .eq. lref(i))) ptstd = i
   10   continue
      endif
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
      time = time*60.0
      if ((time+bastim) .lt. prvtim) bastim = prvtim-time+1.0
      time = time+bastim
      prvtim = time
C---- Set attenuator flag and speed
      if (swidth .ge. 0) then
        ptfilt = 0
      else
        ptfilt = 1
      endif
C
C---- Horizontal and vertical detector(aperture) slit opening
C
C
C---- Scan width
C
      swidth = abs(swidth)
C---- Calculation of transformation coefficients
C---- Scanned background measuring technique
C     Pmt - Total measuring time
      pmt = swidth/spe
C     Bgt - left and right background counting time.
      bgt = 0.5 * ratio * pmt
      cnp = 1.0/spe
      cnb = pmt
      cnl = bgt
      cnr = bgt
      backsn = 1
      nrefl = iabs(nrefl)
C---- Keep one 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
      return
      end
C
      subroutine dtflrd(jj)
C
C===== Subroutine to read a NICOLET diffractometer output file
C
      integer jj, i, trim
C
      include 'common.tlb'
C
      ptscl=1
C
C---- Open input file _REF.DCF)
C
      open( unit=datafl,  
     1  file = compid(1:trim(compid))//datxtn, status = 'OLD', 
     2  access = 'SEQUENTIAL', form = 'UNFORMATTED'
     3  )
CVAX 3  ,readonly)
CUNX 3  )
      if (mstnn .eq. 6) then
        nstep = nstepw
      else
        nstep = 96
      endif
      call initlz
C
C Reads & write all the data file
C
   20 continue
        read(datafl, end = 9999, err = 30) nrefl, ih, ik, il, twoth,
     $  omega, phi, chi, swidth, spe, backl, peak, backr, firl, sir,
     $  time, (profil(i), i = 17, nstep+16)
        call refrd
      go to 20
C
C Error encountered
C
   30 write(*,'(''### Error during reading datafile ###'')')
C
C Finish !!!
C
 9999 close (datafl, status = 'KEEP')
      return
      end
C
C===== Subroutine to read a SYNTEX diffractometer parameter output file
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C===============================================================================
C
C Le diffractometre SYNTEX cree, pour chaque acquisition de donnees, deux
C fichiers, l'un contenant les parametres de la prise de donnee (ouverture
C des detecteurs, ...), et l'autre contenant les intensites mesurees et
C les profils.
C
C Cette subroutine lit les parametres
C
C===============================================================================
      subroutine prflrd

C     Print anything
      character cdum*15
C     File name
      character dcfnme*10
C     Open file type
      character dcopen*10
      character fsd*2
C     User hkl file name
      character hklnme*10
      character lscons*1
C     Anode metal
      character target*2
      character xrysrc*2

C     Loop indices
      integer i, j, trim
      integer ianum
C     Check flag for reference reflections
      integer chknum(12)
      integer dcmode
C     Direction to collect data
      integer idir
C     Exposure time in min
      integer iexp
      integer ifit
      integer iflip
      integer ifrdl
      integer imode
      integer indx
      integer iord(4)
      integer iordr(3)
C     Nb of page in log file
      integer ipage
      integer iprfl
      integer iprint
      integer iprmp
      integer iprof
      integer ipstp
      integer iqus
      integer ircnt
      integer irstp
      integer iscnt
      integer scnhgb
      integer ismpl
      integer ispn
      integer itotl
      integer itype
C     Indices from file flag
      integer iuser
      integer ivect(3)
C     Standard every iwhen refl
      integer iwhen
      integer iwl
      integer iwntme
      integer izm
      integer jtemp
C     Indices of next reflection
      integer ks(3)
      integer ksp(6)
C     Lattice nb
      integer latt
      integer maxfnd
      integer maxp
      integer mlts
      integer msoln
      integer mtape
C     Number of standard
      integer ncheck
      integer ncstp
      integer ngrp
C     Number of reflections left
      integer norefl
      integer nrcnt
C     Integer part of centering reflections
      integer crefl(5, 50)
C     Number of centering reflections
      integer nrf
      integer nsoln
      integer ntrc
C
C     Filter  WARNING- Printed as character
      integer*2 filter
      integer*2 he, hp
      integer*2 iword
      integer*2 lstype
      integer*2 lt
C
C     Used to print systematic absences
      logical print
C
C     Attenuator
      real atnbem
      real axlth(6)
      real baschi
      real basphi
      real baseom
      real basett
      real cfor
C     Chi collision protection limits (NOT READ BY NICOLET.SUB)
      real chibg, chism
C     Chi zero (NOT READ BY NICOLET.SUB)
      real chizro
C     Check value for standard reflections
      real chk(12)
      real chlf
      real cscan
      real ctim
C     Date of data collection
      real date
      real del(4)
      real dev(6)
C     Dispersion (see formula)
      real disp
      real dm
      real dpsi
      real edgmon
C     Low bound used in Wyckoff scans
      real fintm
C     Low & hich correction bounds
      real flcut, fhcut
      real flim1(3), flim2(3)
      real fm3
      real fnl
      real fnono(3)
C     Speed on each circle (NOT READ BY NICOLET.SUB)
      real frat(4)
      real fsig
C     Fastest scan rate
      real fsrt
      real fstt
      real ftpe
      real ofsmic
      real omg1, omg2
C     Omega collision protection limits (NOT READ BY NICOLET.SUB)
      real omgbg, omgsm
C     Omega zero (NOT READ BY NICOLET.SUB)
      real omgzro
      real omspd
      real p3d(10)
      real pcent
C     Phi collision protection limits (NOT READ BY NICOLET.SUB)
      real phibg, phism
      real phis
      real pltim
      real prat
      real pscan
      real psi1, psi2
      real ptim
      real qm
      real ratem
C     Fp part of centering reflections
      real refl(5, 50)
      real refps
      real resmx
C     Range around peak (see formula)
      real rng1, rng2
C     Systematic absences
      real sahkl(12)
      real sahhl(12)
      real sahk(8), sahl(8), sakl(8)
      real sah(6), sak(6), sal(12)
      real sahhlh(2), sahhlr(2), sahblh(2)
      real sclplt
      real slow
C     Slowest scan rate
      real slrt
C     Scan rate
      real srat
C     Estimated time left
      real timest
C     Temporary storage
      real tmp1, tmp2
C     Tolerances (NOT READ BY NICOLET.SUB)
      real tol(4)
C     2 theta collision protection limits (NOT READ BY NICOLET.SUB)
      real tthbg, tthsm
C     2 theta zero (NOT READ BY NICOLET.SUB)
      real tthzro
C     P3 version number
      real versnb
      real xfrat(3)
      real xm
      real ym
      real zm
C
      include 'common.tlb'
C
C****
C     Interpretation of the parameters' file
C****
C
C---- Open input files (_PAR.DCF and _REF.DCF)
C
      open( unit = parmfl, 
     1  file = compid(1:trim(compid))//'_PAR.DCF', status = 'OLD',
     2  form = 'UNFORMATTED')
C
C---- Read title
C
      do 1 i = 1, 60
        title(i:i) = ' '
   1  continue
      read(parmfl) title(1:60), date, versnb
C
C---- Read cell parameters
C
      read(parmfl) (cell(i), i = 1, 6),
     $  (ecell(i), i = 1, 6), ((ub(i, j), j = 1, 3), i = 1, 3),
     $  ((ubinv(i, j), j = 1, 3), i = 1, 3), tthzro, omgzro, chizro,
     $  target, wavel, flam1, flam2, beta, volum, evolum,
     $  edgmon, filter
      if (versnb .ge. 4.12) then
        read(parmfl, err = 90) iprmp, iqus, mstnn, xrysrc, fsd, hp,
     1    he, ispn
      else
        read(parmfl, err = 90) iprmp, iqus, mstnn, fnl, ftpe, ispn
      endif
      go to 91
   90 versnb = 3.43
      mstnn = 0
      ispn = 0
   91 continue
      if (versnb .ge. 3.99) then
        read(parmfl) jtemp, lscons, lstype
      else if (versnb .ne. 3.43) then
        read(parmfl) jtemp
      else
        read(parmfl) iword
        jtemp = iword
      endif
      read(parmfl) fnono, frat, tthsm, tthbg, omgsm,
     $  omgbg, phism, phibg, chism, chibg, imode, lt, basphi,
     $  baschi, basett, baseom, ofsmic
      if ((tthsm .eq. 0.0) .and. (tthbg .eq. 0.0) .and.
     $  (omgsm .eq. 0.0) .and. (omgbg .eq. 0.0)) then
        omgbg = 58.0
        omgsm = 312.0
        tthbg = omgbg + fnono(3)
        tthsm = amod(omgsm - fnono(3), 360.0)
      endif
      read(parmfl) itype, pscan, pltim, iprint,
     $  sclplt, ratem
      read(parmfl) iwl, iord, tol, del, iscnt, cscan,
     $  phis, iprof
      read(parmfl) zm, qm, dm, ptim, xm, ym, prat
      read(parmfl) ivect, ifit, indx, axlth, dev, fm3,
     $  izm, msoln, nsoln
      read(parmfl) omg1, omg2, omspd, ianum
C
C---- Read data collection parameters
C
      read(parmfl) stlmin, stlmax, hmin, kmin, lmin, hmax, kmax, lmax,
     $  ntrc, latt, iuser, nrefl, ks, ipage, iexp, disp, rng1, rng2,
     $  ratio, srat, slrt, fsrt, chlf, cfor, iordr, norefl, timest,
     $  idir
      norefl = iabs(norefl)
      if (versnb .ge. 4.01) then
        read(parmfl, err = 92) ncheck, iwhen, itotl, mtape, ismpl,
     $    iprfl, maxp, dcfnme, dcopen, iwntme, dcmode, hklnme
      else
        read(parmfl, err = 92) ncheck, iwhen, itotl, mtape, ismpl,
     $    iprfl, maxp, dcfnme, dcopen, iwntme, dcmode
      endif
      go to 93
   92 dcfnme = 'SCRATCH'
      dcopen = 'NEW'
      iwntme = 0
   93 continue
      read(parmfl) psi1, psi2, dpsi
      read(parmfl) pcent, ircnt, nrcnt
      read(parmfl) ksp
      read(parmfl) p3d
      read(parmfl) resmx, refps, mlts
      read(parmfl) ncstp, ctim, ngrp, slow
C
C---- Read dead time parameters
C
      read(parmfl) fdeadt, flcut, fhcut, atnbem
      fdeadt = fdeadt * 10**6
      read(parmfl) scnhgb
      read(parmfl) nstepw, ipstp, irstp, fintm,
     $  ifrdl, iflip
C
C---- Read systematic absences
C
      read(parmfl) (sahkl(i), i = 1, 12)
      read(parmfl) (sahk(i), i = 1, 8),
     $  (sahl(i), i = 1, 8), (sakl(i), i = 1, 8)
      read(parmfl) (sah(i), i = 1, 6),
     $  (sak(i), i = 1, 6), (sal(i), i = 1, 12)
      read(parmfl) (sahhl(i), i = 1, 12)
      read(parmfl) (sahhlh(i), i = 1, 2),
     $  (sahblh(i), i = 1, 2)
      read(parmfl) (sahhlr(i), i = 1, 2)
      tmp1 = flim1(1)
      tmp2 = flim1(2)
      read(parmfl, err = 94) maxfnd, fstt, fsig, flim1, flim2, xfrat
      go to 95
   94 if (maxfnd .le. 0) maxfnd = 25
      xfrat(3) = flim1(1)
      flim2(3) = flim1(2)
      flim1(3) = flim1(1)
      flim1(2) = tmp2
      flim1(1) = tmp1
   95 continue
C
C---- Standard reflexions
C
      do 96 i = 1, 10
        href(i) = 0
        kref(i) = 0
        lref(i) = 0
   96 continue
      numref = 0
      do 97 i = 1, 12
        read(parmfl, err = 97) ih, ik, il, chk(i), chknum(i)
        if ((numref .lt. 10) .and. (chknum(i) .gt. 0)) then
          numref = numref+1
          href(numref) = ih
          kref(numref) = ik
          lref(numref) = il
        endif
   97 continue
C
C---- Centering reflexions
C
      nrf = 50
      do 98 i = 1, nrf
        read(parmfl, end = 9999) (crefl(j, i), j = 1, 5),
     $    (refl(j, i), j = 1, 5)
   98 continue
      if (iordr(1) .ne. 0) then
        do 99 i = 1, 3
          iordr(i) = i
   99   continue
      endif
      if (iabs(idir) .ne. 1) idir = 1
      if (edgmon .eq. 0.0) edgmon = 3.1
 9999 close(parmfl, status = 'KEEP')
C
C---- Print parameter file for user
C
      write(*,'('' Do you wish to display the contents of '',
     1 ''the parameter file (Y,N)? '')')
      call fread(000080,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') then
        write(*, '('' Title '', a60, 3X, ''Software version '', F5.2)')
     1    title, versnb
C---- Write cell parameters
        write(*, '(/'' Cell parameters & orientation matrix'')')
        write(*, '(5X, F11.4, 2X, ''a in Angstroems'')') cell(1)
        write(*, '(5X, F11.4, 2X, ''Esd of a'')') ecell(1)
        write(*, '(5X, F11.4, 2X, ''b in Angstroems'')') cell(2)
        write(*, '(5X, F11.4, 2X, ''Esd of b'')') ecell(2)
        write(*, '(5X, F11.4, 2X, ''c in Angstroems'')') cell(3)
        write(*, '(5X, F11.4, 2X, ''Esd of c'')') ecell(3)
        write(*, '(5X, F11.4, 2X, ''alpha in degrees'')') cell(4)
        write(*, '(5X, F11.4, 2X, ''Esd of alpha'')') ecell(4)
        write(*, '(5X, F11.4, 2X, ''beta  in degrees'')') cell(5)
        write(*, '(5X, F11.4, 2X, ''Esd of beta'')') ecell(5)
        write(*, '(5X, F11.4, 2X, ''gamma in degrees'')') cell(6)
        write(*, '(5X, F11.4, 2X, ''Esd of gamma'')') ecell(6)
        write(*, '(5X, F11.4, 2X, ''Volume in Angstroems**3'')')
     1    volum
        write(*, '(5X, F11.4, 2X, ''Esd of volume'')') evolum
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 1)'')')
     1    ub(1, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 2)'')')
     1    ub(1, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 3)'')')
     1    ub(1, 3)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 1)'')')
     1    ub(2, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 2)'')')
     1    ub(2, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 3)'')')
     1    ub(2, 3)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 1)'')')
     1    ub(3, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 2)'')')
     1    ub(3, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 3)'')')
     1    ub(3, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 1)'')') ubinv(1, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 2)'')') ubinv(1, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 3)'')') ubinv(1, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 1)'')') ubinv(2, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 2)'')') ubinv(2, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 3)'')') ubinv(2, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 1)'')') ubinv(3, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 2)'')') ubinv(3, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 3)'')') ubinv(3, 3)
        cdum = '    Constrained'
        if ((lscons .ne. 'C') .and. (lscons .ne. 'c'))
     1    cdum = '  Unconstrained'
        write(*, '(1X, A15, 2X, ''Orientation matrix'')') cdum
C---- Write radiation & filter stuff
        write(*, '(/'' Radiation and filter'')')
        write(*, '(14X, A2, 2X, ''Anode element'')') target
        write(*, '(5X, F11.6, 2X, ''Mean wavelength'')') wavel
        write(*, '(5X, F11.6, 2X, ''Kalpha1 wavelength'')') flam1
        write(*, '(5X, F11.6, 2X, ''Kalpha2 wavelength'')') flam2
        write(*, '(5X, F11.6, 2X, ''Kbeta'')') beta
        if (filter .ne. 0) then
          write(*, '(14X, A2, 2X, ''Absorption filter element'')')
     1      filter
        endif
        write(*, '(14X, I2, 2X, ''Wavelength pointer'')') iwl
C---- Write lattice type & systematic absences
        write(*, '(/'' Lattice type & systematic absences'')')
        cdum = '              ?'
        go to (111, 112, 113, 114, 115, 116, 117, 118) latt+1
  111   cdum = '              P'
        go to 119
  112   cdum = '              A'
        go to 119
  113   cdum = '              B'
        go to 119
  114   cdum = '              C'
        go to 119
  115   cdum = '              F'
        go to 119
  116   cdum = '              I'
        go to 119
  117   cdum = '           Robv'
        go to 119
  118   cdum = '           Rrev'
  119   write(*, '(1X, A15, 2X, ''Lattice type'')') cdum
        print = .false.
        do 130 i = 1, 12
          if (sahkl(i) .ne. 0.0) print = .true.
  130   continue
        if (print) then
          do 135 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hkl'')')
     1        sahkl(i)
  135     continue
        end if
        print = .false.
        do 140 i = 1, 8
          if (sahk(i) .ne. 0.0) print = .true.
  140   continue
        if (print) then
          do 145 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hk'')')
     1        sahk(i)
  145     continue
        end if
        print = .false.
        do 150 i = 1, 8
          if (sahl(i) .ne. 0.0) print = .true.
  150   continue
        if (print) then
          do 155 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hl'')')
     1        sahl(i)
  155     continue
        end if
        print = .false.
        do 160 i = 1, 8
          if (sakl(i) .ne. 0.0) print = .true.
  160   continue
        if (print) then
          do 165 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type kl'')')
     1        sakl(i)
  165     continue
        end if
        print = .false.
        do 170 i = 1, 6
          if (sah(i) .ne. 0.0) print = .true.
  170   continue
        if (print) then
          do 175 i = 1, 6
            write(*, '(5X, F11.2, 2X, ''Systematic absence type h'')')
     1        sah(i)
  175     continue
        end if
        print = .false.
        do 180 i = 1, 6
          if (sak(i) .ne. 0.0) print = .true.
  180   continue
        if (print) then
          do 185 i = 1, 6
            write(*, '(5X, F11.2, 2X, ''Systematic absence type k'')')
     1        sak(i)
  185     continue
        end if
        print = .false.
        do 190 i = 1, 12
          if (sal(i) .ne. 0.0) print = .true.
  190   continue
        if (print) then
          do 195 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type l'')')
     1        sal(i)
  195     continue
        end if
        print = .false.
        do 200 i = 1, 12
          if (sahhl(i) .ne. 0.0) print = .true.
  200   continue
        if (print) then
          do 205 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hhl'')')
     1        sahhl(i)
  205     continue
        end if
        print = .false.
        do 210 i = 1, 2
          if (sahhlh(i) .ne. 0.0) print = .true.
  210   continue
        if (print) then
          do 215 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhl_H'')') sahhlh(i)
  215     continue
        end if
        print = .false.
        do 220 i = 1, 2
          if (sahblh(i) .ne. 0.0) print = .true.
  220   continue
        if (print) then
          do 225 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhbl_H'')') sahblh(i)
  225     continue
        end if
        print = .false.
        do 230 i = 1, 2
          if (sahhlr(i) .ne. 0.0) print = .true.
  230   continue
        if (print) then
          do 235 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhl_R'')') sahhlr(i)
  235     continue
        end if
C---- Write data collection parameters
        write(*, '(/'' Data collection parameters'')')
        cdum = '        Unknown'
        if (mstnn .eq. 1) cdum = '          Omega'
        if (mstnn .eq. 2) cdum = '   Theta:2Theta'
        if (mstnn. eq. 6) cdum = '          Wycoff'
        write(*, '(1X, A15, 2X, ''Data collection type'')') cdum
        cdum = '              ?'
        go to (251, 252, 250, 253, 254, 255) imode+3
  251     cdum = ' Parallel(user)'
          go to 250
  252     cdum = 'Bisecting(user)'
          go to 250
  253     cdum = '      Bisecting'
          go to 250
  254     cdum = '      Bisecting'
          go to 250
  255     cdum = '       Parallel'
  250   continue
        write(*, '(1X, A15, 2X, ''Data collectiom mode '')') cdum
        write(*, '(14x, i2, 2x, ''Low temperature mode '',
     1    ''(0=no low temp)'')') lt
        write(*, '(6X, A10, 2X, ''Data collection file name'')') dcfnme
        write(*, '(10X, A6, 2X, ''Open mode'')') dcopen
C---- Data collection bounds
        write(*, '(/'' Data collection bounds'')')
        write(*, '(5X, F11.2, 2X, ''2 theta min'')') stlmin
        write(*, '(5X, F11.2, 2X, ''2 theta max'')') stlmax
        write(*, '(12X, I4, 2X, ''H min'')') hmin
        write(*, '(12X, I4, 2X, ''K min'')') kmin
        write(*, '(12X, I4, 2X, ''L min'')') lmin
        write(*, '(12X, I4, 2X, ''H max'')') hmax
        write(*, '(12X, I4, 2X, ''K max'')') kmax
        write(*, '(12X, I4, 2X, ''L max'')') lmax
        write(*, '(12x, i4, 2x, ''Number of last reflection '',
     1    ''proceeded'')') nrefl
        write(*, '(12X, I4, 2X, ''H of last reflection proceeded'')')
     1    ks(1)
        write(*, '(12X, I4, 2X, ''K of last reflection proceeded'')')
     1    ks(2)
        write(*, '(12X, I4, 2X, ''L of last reflection proceeded'')')
     1    ks(3)
        write(*, '(8X, I8, 2X, ''Number of reflections proceeded'')')
     1    norefl
        write(*, '(8X, I8, 2X, ''Total exposure time in min'')') iexp
        write(*, '(5X, F11.4, 2X, ''Remaining time in hours'')') timest
        write(*, '(12X, I4, 2X, ''H of first reflection'')') ksp(1)
        write(*, '(12X, I4, 2X, ''K of first reflection'')') ksp(2)
        write(*, '(12X, I4, 2X, ''L of first reflection'')') ksp(3)
        write(*, '(12X, I4, 2X, ''H of last  reflection'')') ksp(4)
        write(*, '(12X, I4, 2X, ''K of last  reflection'')') ksp(5)
        write(*, '(12X, I4, 2X, ''L of last  reflection'')') ksp(6)
        write(*, '(14X, I2, 2X, ''H variation order'')') iordr(1)
        write(*, '(14X, I2, 2X, ''K variation order'')') iordr(2)
        write(*, '(14X, I2, 2X, ''L variation order'')') iordr(3)
        write(*, '(14X, I2, 2X, ''Direction of index variation'')')
     1    idir
        write(*, '(14X, I2, 2X, ''Index restriction code'')') ntrc
        if (iuser .eq. 1) then
          write(*, '(6x, a10, 2x, ''Selected hkl data collection '',
     1      ''file name'')') hklnme
        endif
        if ((dpsi .ne. 0.0) .or. (psi1 .ne. 0.0)) then
          write(*, '(5X, F11.2, 2X, ''Low  psi bound'')') psi1
          write(*, '(5X, F11.2, 2X, ''High psi bound'')') psi2
          write(*, '(5X, F11.4, 2X, ''Delta psi angle'')') dpsi
        endif
C---- Data collection parameters
        write(*, '(/'' Data collection parameters (part 2)'')')
        write(*, '(5X, F11.4, 2X, ''Multiplicative factor for range'')')
     1    disp
        write(*, '(5X, F11.4, 2X, ''Below Kalpha1'')') rng1
        write(*, '(5X, F11.4, 2X, ''Above Kalpha2'')') rng2
        if (srat .eq. 0.0) then
          write(*, '(5X, F11.4, 2X, ''Slowest scan speed'')') slrt
          write(*, '(5X, F11.4, 2X, ''Fastest scan speed'')') fsrt
        else
          write(*, '(5X, F11.4, 2X, ''Fixed scan speed'')') srat
        endif
        write(*, '(5X, F11.6, 2X, ''Background/scan measuring time'')')
     1    ratio
        write(*, '(5X, F11.4, 2X, ''Low  I or |F|'')') chlf
        write(*, '(5X, F11.4, 2X, ''High I or |F|'')') cfor
        if (mstnn .eq. 6) then
          write(*, '(12X, I4, 2X, ''Total steps in Wycoff scan'')')
     1      nstepw
          write(*, '(12X, I4, 2X, ''Calculated steps'')') ipstp
          write(*, '(12X, I4, 2X, ''Extra steps'')') irstp
          write(*, '(5X, F11.2, 2X, ''Intensity'')') fintm
          write(*, '(14X, I2, 2X, ''Collect Friedel pair'')') ifrdl
          write(*, '(14X, I2, 2X, ''Start with hkl or -h-k-l'')') iflip
        endif
        write(*, '(14x, i2, 2x, ''1=profiled prescan, 2=intensity, '',
     1    ''<0=Profiled postscan'')') ispn
C---- Diffractometer parameters
        write(*, '(/'' Diffractometer parameters'')')
        write(*, '(12X, I4, 2X, ''Number of steps'')') ncstp
        write(*, '(12X, I4, 2X, ''Number of group'')') ngrp
        write(*, '(5X, F11.4, 2X, ''Second/step'')') ctim
        write(*, '(5X, F11.4, 2X, ''Min drive speed'')') slow
        write(*, '(5X, F11.4, 2X, ''Dead time in microsec'')') fdeadt
        write(*, '(5X, F11.2, 2X, ''Correction low  bound'')') flcut
        write(*, '(5X, F11.2, 2X, ''Correction high bound'')') fhcut
        if (atnbem .ne. 0.0) then
          write(*, '(5X, F11.6, 2X, ''Attenuation'')') atnbem
        end if
        write(*, '(5x, f11.6, 2x, ''2 theta zero '',
     1    ''(calibration value) '')') tthzro
        write(*, '(5x, f11.6, 2x, ''omega   zero '',
     1    ''(calibration value) '')') omgzro
        write(*, '(5x, f11.6, 2x, ''chi     zero '',
     1    ''(calibration value) '')') chizro
        write(*, '(5X, F11.2, 2X, ''2 Theta collision low  bound'')')
     1    tthsm
        write(*, '(5X, F11.2, 2X, ''2 Theta collision high bound'')')
     1    tthbg
        write(*, '(5X, F11.2, 2X, ''Omega   collision low  bound'')')
     1    omgsm
        write(*, '(5X, F11.2, 2X, ''Omega   collision high bound'')')
     1    omgbg
        write(*, '(5X, F11.2, 2X, ''Phi     collision low  bound'')')
     1    phism
        write(*, '(5X, F11.2, 2X, ''Phi     collision high bound'')')
     1    phibg
        write(*, '(5X, F11.2, 2X, ''Chi     collision low  bound'')')
     1    chism
        write(*, '(5X, F11.2, 2X, ''Chi     collision high bound'')')
     1    chibg
        write(*, '(5X, F11.4, 2X, ''2 Theta fastest speed'')') frat(1)
        write(*, '(5X, F11.4, 2X, ''Omega   fastest speed'')') frat(2)
        write(*, '(5X, F11.4, 2X, ''Phi     fastest speed'')') frat(3)
        write(*, '(5X, F11.4, 2X, ''Chi     fastest speed'')') frat(4)
        write(*, '(5X, F11.4, 2X, ''Microscope offset'')') ofsmic
        write(*, '(5x, f11.2, 2x, ''2 theta home value for '',
     1    ''centering crystal in beam'')') basett
        write(*, '(5x, f11.2, 2x, ''Omega   home value for '',
     1    ''centering crystal in beam'')') baseom
        write(*, '(5x, f11.2, 2x, ''Phi     home value for '',
     1    ''centering crystal in beam'')') basphi
        write(*, '(5x, f11.2, 2x, ''Chi     home value for '',
     1    ''centering crystal in beam'')') baschi
C---- Write standard reflections
        write(*, '(/'' Standard reflections'')')
        write(*, '(14X, I2, 2X, ''Number of check reflections'')')
     1    numref
        write(*, '(12X, I4, 2X, ''Interval between checks'')') iwhen
        write(*, '(12X, I4, 2X, ''Interval between checks in min'')')
     1    iwntme
        write(*, '(5X, F11.4, 2X, ''Recenter if intensity below'')')
     1    pcent
        write(*, '(/''   h   k   l'', 5x, ''low bound'', 3x, ''check''
     1    )')
        do 300 i = 1, numref
          write(*, '(3(1X, I3), 4X, F10.2, 5X, I3)')
     1      href(i), kref(i), lref(i), chk(i), chknum(i)
  300   continue
C---- Write centering reflections
        write(*, '(/'' Centering reflections'')')
        write(*, '(14x, i2, 2x, ''Centering: 1st pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(1)
        write(*, '(14x, i2, 2x, ''Centering: 2nd pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(2)
        write(*, '(14x, i2, 2x, ''Centering: 3rd pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(3)
        write(*, '(14x, i2, 2x, ''Centering: 4th pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(4)
        write(*, '(5X, F11.4, 2X, ''Tolerance on 2 Theta'')') tol(1)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Omega'')') tol(2)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Phi'')') tol(3)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Chi'')') tol(4)
        write(*, '(5X, F11.4, 2X, ''2 Theta scan range'')') del(1)
        write(*, '(5X, F11.4, 2X, ''Omega   scan range'')') del(2)
        write(*, '(5X, F11.4, 2X, ''Phi     scan range'')') del(3)
        write(*, '(5X, F11.4, 2X, ''Chi     scan range'')') del(4)
        write(*, '(5X, F11.4, 2X, ''Centering speed'')') cscan
        write(*, '(5X, F11.4, 2X, ''Phi search speed'')') phis
        write(*, '(10X, I6, 2X, ''Min intensity for phi search'')')
     1    iscnt
        write(*, '(10X, I6, 2X, ''Background intensity'')') scnhgb
        go to (311, 312, 313, 314) ircnt+1
  311     cdum = '           None'
        go to 310
  312     cdum = '       Recenter'
        go to 310
  313     cdum = '  Least squares'
        go to 310
  314     cdum = 'Data Collection'
  310   write(*, '(1X, A15, 2X, ''First operation'')') cdum
        write(*, '(10X, I6, 2X, ''Interval between recenter'')') nrcnt
        write(*, '(14X, I2, 2X, ''Max nb of reflections'')') maxfnd
        write(*, '(5X, F11.6, 2X, ''Search time [sec/step]'')') fstt
        write(*, '(5X, F11.4, 2X, ''Refl significance fact'')') fsig
        write(*, '(5X, F11.2, 2X, ''2 Theta min'')') flim1(1)
        write(*, '(5X, F11.2, 2X, ''Omega   min'')') flim1(2)
        write(*, '(5X, F11.2, 2X, ''Chi     min'')') flim1(3)
        write(*, '(5X, F11.2, 2X, ''2 Theta max'')') flim2(1)
        write(*, '(5X, F11.2, 2X, ''Omega   max'')') flim2(2)
        write(*, '(5X, F11.2, 2X, ''Chi     max'')') flim2(3)
        write(*, '(5X, F11.4, 2X, ''2 Theta max speed'')') xfrat(1)
        write(*, '(5X, F11.4, 2X, ''Omega   max speed'')') xfrat(2)
        write(*, '(5X, F11.4, 2X, ''Chi     max speed'')') xfrat(3)
C---- Orientation matrix calculation
        write(*, '(/'' Orientation matrix calculation'')')
        write(*, '(12X, I4, 2X, ''Number to reflections to be used'')')
     1    ifit
        write(*, '(12X, I4, 2X, ''Max value allowed for indexing'')')
     1    indx
        write(*, '(5X, F11.4, 2X, ''Min axis a value'')') axlth(1)
        write(*, '(5X, F11.4, 2X, ''Max axis a value'')') axlth(2)
        write(*, '(5X, F11.4, 2X, ''Min axis b value'')') axlth(3)
        write(*, '(5X, F11.4, 2X, ''Max axis b value'')') axlth(4)
        write(*, '(5X, F11.4, 2X, ''Min axis c value'')') axlth(5)
        write(*, '(5X, F11.4, 2X, ''Max axis c value'')') axlth(6)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis a value'')')
     1    dev(1)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis a value'')')
     1    dev(2)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis b value'')')
     1    dev(3)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis b value'')')
     1    dev(4)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis c value'')')
     1    dev(5)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis c value'')')
     1    dev(6)
        write(*, '(5X, F11.4, 2X, ''Max index deviation for integer'')')
     1    fm3
        write(*, '(12X, I4, 2X, ''Number of angles'')') izm
        write(*, '(12X, I4, 2X, ''Number solutions'')') msoln
        write(*, '(12X, I4, 2X, ''Number of input solutions'')') nsoln
        write(*, '(14X, I2, 2X, ''Offset for centering reflections'')')
     1    jtemp
        write(*, '(/''  c    h    k    l   wl     2 theta'', 7x,
     1    ''Omega         Phi         Chi           I'')')
        do 350 i = 1, nrf
          write(*, '(5(I3, 2X), 5(F10.4, 2X))')
     1      (crefl(j, i), j = 1, 5), (refl(j, i), j = 1, 5)
  350   continue
C---- Control parameters
        write(*, '(/'' Control parameters'')')
        go to (321, 322, 323) ismpl+1
  321     cdum = '            Off'
        go to 320
  322     cdum = '            Log'
        go to 320
  323     cdum = '          Print'
  320   write(*, '(1X, A15, 2X, ''Log file'')') cdum
        write(*, '(14X, I2, 2X, ''Number of steps to be printed'')')
     1    iprfl
        write(*, '(5X, F11.6, 2X, ''Plot scale'')') sclplt
        go to (331, 332, 333) iprof+1
  331     cdum = '           None'
        go to 330
  332     cdum = '           Line'
        go to 330
  333     cdum = '         Screen'
  330   write(*, '(1X, A15, 2X, ''Profile flag '')') cdum
        if (ismpl .eq. 1) then
          write(*, '(10X, I6, 2X, ''Number of pages in log file '')')
     1      ipage
        endif
C---- Check profiles
        write(*, '(/'' Plot one reflection parameters'')')
        go to (361, 362, 363, 364), itype+1
  361     cdum = '           None'
        go to 360
  362     cdum = '   Omega:2Theta'
        go to 360
  363     cdum = '         2Theta'
        go to 360
  364     cdum = '          Omega'
        go to 360
C 365     cdum = '            Phi'
C       go to 360
C 366     cdum = '            Chi'
  360   write(*, '(1X, A15, 2X, ''Scan type'')') cdum
        go to (341, 342, 343, 344) iprint+1
  341     cdum = '      Peak Plot'
        go to 340
  342     cdum = '         Counts'
        go to 340
  343     cdum = '      Ratemeter'
        go to 340
  344     cdum = '    Powers of 2'
  340   write(*, '(1X, A15, 2X, ''Screen output'')') cdum
        write(*, '(5X, F11.4, 2X, ''Scan range'')') pscan
        write(*, '(5X, F11.6, 2X, ''Sec/step'')') pltim
        write(*, '(5X, F11.4, 2X, ''Drive speed'')') ratem
C---- Photo parameters
        write(*, '(/'' Rotation photo parameters'')')
        write(*, '(5X, F11.4, 2X, ''Exposure time in min'')') ptim
        write(*, '(5X, F11.4, 2X, ''Distance form crystal to film'')')
     1    dm/2.0
        write(*, '(5X, F11.2, 2X, ''2 Theta position'')') zm
        write(*, '(5X, F11.2, 2X, ''Chi     position'')') qm
        write(*, '(5X, F11.2, 2X, ''Phi start'')') xm
        write(*, '(5X, F11.2, 2X, ''Phi end'')') ym
        write(*, '(5X, F11.4, 2X, ''Phi rate'')') prat
        write(*, '(/'' Axial photo parameters'')')
        write(*, '(5X, F11.2, 2X, ''Omega start'')') omg1
        write(*, '(5X, F11.2, 2X, ''Omega end'')') omg2
        write(*, '(5X, F11.4, 2X, ''Omega rate'')') omspd
C        WRITE(formfl, '('' Iprmp      I4  1 '', I6, 74X,
C     $    ''?'')') Iprmp
C        WRITE(formfl, '('' Iqus       I4  1 '', I6, 74X,
C     $    ''?'')') Iqus
C        WRITE(formfl, '('' XrySrc     C   2 '', A2, 78X,
C     $    ''?'')') XrySrc
C        WRITE(formfl, '('' FSD        C   2 '', A2, 78X,
C     $    ''?'')') FSD
C        WRITE(formfl, '('' HP         I2  1 '', I6, 74X,
C     $    ''?'')') HP
C        WRITE(formfl, '('' HE         I2  1 '', I6, 74X,
C     $    ''?'')') HE
C        WRITE(formfl, '('' Fnl        R   1 '', F10.4 ,70X,
C     $    ''?'')') Fnl
C        WRITE(formfl, '('' Ftpe       R   1 '', F10.4, 70X,
C     $    ''?'')') Ftpe
C        WRITE(formfl, '('' LSType     I   2 '', I4, 76X,
C     $    ''?'')') LSType
C      WRITE(formfl, '('' Fnono      R   3 '', 3(F6.2, 2X), 56X,
C     $  ''?'')') Fnono
C      WRITE(formfl, '('' Ivect      I4  3 '', 3(I6, 2X), 56X,
C     $  ''?'')') ivect
C      WRITE(formfl, '('' Ianum      I4  1 '', I6, 74X,
C     $  ''?'')') Ianum
C      WRITE(formfl, '('' Itotl      I4  1 '', I6, 74X,
C     $  ''?'')') Itotl
C      WRITE(formfl, '('' Mtape      I4  1 '', I6, 74X,
C     $  ''?'')') Mtape
C      WRITE(formfl, '('' Maxp       I4  1 '', I6, 74X,
C     $  ''?'')') Maxp
C      WRITE(formfl, '('' P3d        R  10 '', 10(F6.2, 2x),
C     $  ''?'')') p3d
C      WRITE(formfl, '('' Resmx      R   1 '', F10.4, 70X,
C     $  ''?'')') Resmx
C      WRITE(formfl, '('' Refps      R   1 '', F10.4, 70X,
C     $  ''?'')') Refps
C      WRITE(formfl, '('' Mlts       I4  1 '', I6, 74X,
C     $  ''?'')') Mlts
C
C---- End of normal job
C
      endif
      do 10 i = 61, 66
        title(i:i) = ' '
   10 continue
      if (iwl .eq. 2) wavel = flam1
      if (iwl .eq. 3) wavel = flam2
      stlmin = sin(dar*stlmin/2.0)/wavel
      stlmax = sin(dar*stlmax/2.0)/wavel
      return
      end
      subroutine choos
C
C----- Set up fixed characteristics for NICOLET diffractometer
C
      include 'common.tlb'
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 a Nicolet/Syntex diffractometer
C
      difnam = 'S'
      device = 'Siemens_Nicolet_P21_4_circles_diffracto '
      softwr = 'Siemens_Nicolet_P3R3_V4.11              '
C
C---- Tell the programme the name of the address file
      addnam = 'address.lst                     '
C
C---- Set extension for data file name
      datxtn = '_REF.DCF'
C
C---- Set dihedral angle
      fmnhed = 0.0
C
C---- Default wavelength pointer for cell dimension determination, alpha 1
      cdefwv = 1
C
C---- Options for Quesan
      ask(2) = .false.
      ask(3) = .false.
      ask(4) = .false.
C
      return
      end
      subroutine initlz
C
C**** Subroutine for the NICOLET P21 to convert diffractometer parameter
C     information into the standardized form used by these programmes
C     The UB will be converted to the International Tables Volume 4
C     definition.
C****
      integer i,j,k
      real ubdet
C
      include 'common.tlb'
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(i) = 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---- 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)
          endif
 201    continue
        cwave = flam1
C
C---- Horizontal and vertical apertures
C     ouverture : 2.5 mm environ, distance 135 mm environ
C
      detvt = 2.5/135.0
      dethz = 2.5/135.0
C
C---- Output the conditions
      call conout
C
      return
      end
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Siemens-Nicolet P21 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, In the Press. 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 Nicolet P21 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 Nicolet P21 is thus a DSHDF type I diffractometer and
C     hence W = 1 according to DSHDF equation (11). (W and 1 being 3x3
C     matrices).
C
C     The Nicolet P21 uses the Eulerian cradle geometry described in detail
C     in the Nicolet P21
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
      integer  i,j
      real r(3), s(3), sinq, cosq
      real w(9)
      data w/1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
      include 'common.tlb'
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---- *x(-chi)
      sinq = sin(dar*chi)
      cosq = cos(dar*chi)
      s(1) =  r(1)
      s(2) =  r(2)*cosq + r(3)*sinq
      s(3) = -r(2)*sinq + r(3)*cosq
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
      do 115 i = 1, 3
        s(i) = 0.0
        do 115 j = 1, 3
          s(i) = s(i) + w(3*(i-1)+j)*r(j)
  115 continue
C
C---- Crystal based azimuth
      psicry = rad*atan2(s(1),s(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



------cut here------------------------------------------------------------------


C
C==== PW1100 Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
C
      subroutine refrd(itype)
C**** =======================================================================
C     Routine to interpret reflection data from a PW1100. Itype indicates 0
C     for an ordinary reflection and 1 for a reference reflection.
C     All data is output by calling REFOUT
C**** =======================================================================
      integer itype,i
      real pmt
      logical oka
C
      include 'common.tlb'
C
      if (itype .eq. 0) go to 2005
      go to 3005
C---- Ordinary reflection line
 2005 continue
C
C---- In the DIFRAC module for the Philips PW1100 diffractometer, the
C---- UB, UBINV and UBDET variables contain the values used on the
C---- diffractometer even if the determinant terms out to be negative.
C---- So for the psi calculation the diffractometer hkls are used and
C---- are transformed for output.  It has to be done this way because
C---- the PW1100 does not output the setting angles used for each
C---- reflection.
C
C---- Calculate crystal based psi value
      call psic
C---- Transform output h k l for negative UB
      if (trans) then
        ih = -ih
        ik = -ik
        il = -il
      endif
      peak  = values(ptcsum)
      backl = values(ptcb1)
      backr = values(ptcb2)
      time = time + 1.0
      nrefl = nrefl + 1
      ptstd = 0
C---- Calculation of transformation coefficients
C     Pmt - Peak-with-background measuring time
      pmt = values(ptnsc)*values(ptswd)/spe
C     Bgt - left and right background counting time.
C     For bmo = 0 the time is that in values(ptbgt).
C     For bmo = 1 its half the peak scan time
C     For bmo = 2 the time is half that in values(ptbgt).
      if (bmo. eq. 0) then
        bgt = values(ptbgt)
      else if (bmo .eq. 1) then
        bgt = 0.5 * pmt
      else if (bmo .eq. 2) then
        bgt = 0.5 * values(ptbgt)
      endif
      cnp = pmt/values(ptswd)
      cnb = pmt
      cnl = bgt
      cnr = bgt
C---- Scan width
      swidth = values(ptswd)
C---- Scan ratio
      if (abs(theta) .le. sct) then
C---- omega scans
        scnrat = 0.0
      else
C---- omega - 2-theta scans
        scnrat = 2.0
      endif
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
      sintl = sint/wavel
      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
      return
C---- Reference reflection line
 3005 continue
      do 3095 i = 1 , numref
        ih = href(i)
        ik = kref(i)
        il = lref(i)
        peak = refint(i)
        backl = 0.0
        backr = 0.0
        psicry = psirf(i)
        time = time + 1.0
        nrefl = nrefl + 1
        ptfilt = fltref(i)
        ptstd = i
        cnp = cnprf(i)
        cnb = cnbrf(i)
        cnl = 0.0
        cnr = 0.0
C---- Scan width and scan ratio
        swidth = swdrf(i)
        scnrat = scnrf(i)
        avint(i) = avint(i) + peak
        sdint(i) = sdint(i) + peak*peak
        count(i) = count(i) + 1.0
C
        call refout
 3095 continue
      return
      end
C
      subroutine dtflrd(jj)
C
C==== Subroutine to read a PW11OO diffractometer output file
C
C
      integer i,j,k,jj, trim
C
      include 'common.tlb'
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
C---- Initialisations
C
      firstr = .false.
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  )             
CVAX 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))//digit(i)//digit(j)//datxtn,
     2  status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     3  )              
CVAX 3  ,readonly)
CUNX 3  )
        write(*,'(/'' Reading file: '',A)')
     1   compid(1:trim(compid))//digit(i)//digit(j)//datxtn
      endif
C****
C     Interpretation of the .DAT file
C****
 1000 continue
      line = ' '
      read(datafl, '(a80)', end = 2000) line
C---- Look for PAR parameters
      if (line(36:38) .eq. 'SCT') then
        parclp = 1
        read (line, '(7x,f7.4,7x,f5.2,1x,f5.2,7x,f7.2)')
     1   spe, swda, swdb, sct
      elseif (line(32:34) .eq. 'MXN') then
        parclp = parclp + 2
        read (line, '(35x,i3)') mxn
      elseif (line(14:16) .eq. 'BMO') then
        parclp = parclp + 4
        read (line, '(9x,i1,9x,i1)') smo, bmo
        if (bmo. eq. 0) read (line, '(20x,f6.0)') bgt
      elseif (line(4:6) .eq. 'REF') then
        parclp = parclp + 8
        read (line, '(9x,i1,1x,f5.0)') numref, sinter
        read(datafl, '(a80)', end = 2000) line
        do 1005 i = 1 ,numref
          read(datafl, '(a80)', end = 2000) line
          read(line, '(6x,3i4,f8.4,i4)')
     1    href(i), kref(i), lref(i), speref(i), fltref(i)
 1005   continue
      elseif (line(4:6) .eq. 'SUB') then
        parclp = parclp + 16
        do 1010 i = 1 ,3
          read(datafl, '(a80)', end = 2000) line
          read(datafl, '(a80)', end = 2000) line
          read(line, '(3f12.6)') (ub(i,k), k = 1,3)
 1010   continue
      elseif (line(18:20) .eq. 'DET') then
        parclp = parclp + 32
        read(line, '(6x,f8.4,6x,2f6.2)') wavel, dethz, detvt
        if (parclp .eq. 63) then
          firstr = .true.
        else
        endif
      elseif (parclp .eq. 63) then
        if (line(8:10) .eq. 'RFA') then
          read(line, '(18x,f8.0,12x,f8.0,12x,f8.0)')
     1     (refint(i), i = 1 , numref)
          if (firstr) call initlz
          call refrd(1)
        else
C---- Check only for the characters in digit
          ok = .true.
          jj = 0
          do 1020 i = 1, 80
          do 1015 k = 0, 11
            if (line(i:i) .eq. digit(k)) go to 1019
 1015     continue
C---- Make a concession for column 22 when there is a variable scan width
          if ((line(22:22) .eq. '.') .and. (swdb .ne. 0.0)) go to 1020
          ok = .false.
          go to 1020
 1019     if (k .eq. 11) jj = jj+1
 1020     continue
C---- Skip blank lines
          if (jj.eq.80) ok = .false.
          if (ok) then
            if (firstr) call initlz
            read (line, fmt) ih, ik, il, psi, ptfilt,
     1                       (values(i),i=1,lfnorm)
            call refrd(0)
          else
          endif
        endif
      else
      endif
      go to 1000
C
C---- Finish off this data file
C
 2000 continue
      close (unit = datafl, status = 'KEEP')
      return
      end
C
      subroutine prflrd
C
C===== Subroutine to read a PW1100 diffractometer parameter output file
C
C===== On the Philips PW1100 there is no separate parameter output file
C      produced directly by the diffractometer. However we edit the .LOG
C      file to produce a file containing input lines for standalone a LATCON.
C
C      The file has two lines to skip at the top and then lines:
C      h, k, l, omega positive, omega negative
C
      include 'common.tlb'
      integer nrec, i, trim
      real omegp, omegn
      character*80  zharin
      logical zresnt
C
C==== Read the <Compid>.tth file of a PW1100
C
C---- Inquire <Compid>.tth and open it if present
C
      inquire(file=compid(1:trim(compid))//'.tth',
     1 exist=zresnt)
      if(.not.zresnt) then
        write(*,'(''   File: '',A,''.tth not available'')')
     1   compid(1:trim(compid))
      else
        write(*,'(/'' Reading file: '',A,''.tth'')')
     1   compid(1:trim(compid))
        open( unit=parmfl,file=compid(1:trim(compid))//'.tth',
     1  status='OLD', access='SEQUENTIAL', form='FORMATTED',err = 70
     2  )           
CVAX 2  ,readonly) 
CUNX 2  )
C
C---- Header record
C
        read(unit=parmfl,fmt=*)
C
C---- Default wavelength index
C     0 = already defined default, 1 = alpha 1, 2 = alpha 2,
C     anything else = mean alpha
C
        zharin = ' '
        read(unit=parmfl,fmt='(a80)',end=70) zharin
        zharin(78:80) = ' 0 '
        read(unit=zharin,fmt=*) i
        if      (i.eq.1) then
          cdefwv = 1
        else if (i.eq.2) then
          cdefwv = 2
        else if (i.ne.0) then
          cdefwv = 0
        endif
C
C---- Treat  Reflection records
C
        do 2000 nrec = 1, 100
          zharin = ' '
          read(unit=parmfl,fmt='(a80)',end=70) zharin
          zharin(78:80) = ' 0 '
          read(unit=zharin,fmt=*)
     1     hang(nrec), kang(nrec), lang(nrec), omegp, omegn, wang(nrec)
C---- Pick up the two theta angles for "Good" angle reflections
          numang = numang + 1
          twoang(nrec) = abs(omegp) + abs(omegn)
          if (wang(nrec).eq.0) then
            wang(nrec) = cdefwv
          else if (wang(nrec).gt.2) then
            wang(nrec) = 0
          endif
          ctthn = amin1(twoang(nrec),ctthn)
          ctthx = amax1(twoang(nrec),ctthx)
 2000   continue
 70     close(unit=parmfl,status='KEEP')
        nmcell = numang
      endif
C
      return
      end
C
      subroutine choos
C
C----- Set up fixed characteristics for PW1100 diffractometer
C
       include 'common.tlb'
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 a Philips PW1100 diffractometer
      difnam = 'W'
      device = 'Philips_PW1100_4-circle_diffractometer  '
      softwr = 'Philips_PW1100_April_1974               '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = 'pck$:[cristal]address.lst       '
CNORMAL addnam = 'address.lst                     '
      addnam = 'pck$:[cristal]address.lst       '
C
C---- Set extension for data file name
      datxtn = '.dat'
C
C---- Set dihedral angle
      fmnhed = 90.0
C
C---- Default wavelength pointer for cell dimension determination, mean lambda
      cdefwv = 0
C
C---- Options for Quesan
C
      return
      end
C
      subroutine initlz
C
C**** Subroutine to interpret a complete PAR instruction in memory
C     for the PW1100.
C****
      integer i,j,k,iz
      real ubdet, temp2, temp3
      logical oka
C
      include 'common.tlb'
C
C---- Check PAR values
C
      if (spe .le. 0.0)  goto 999
      if (spe .gt. 0.32) goto 999
      if (swda .le. 0.0) goto 999
      if ((smo .lt. 1) .or. (smo .gt. 2)) goto 999
      if ((bmo. lt. 0) .or. (bmo .gt. 2)) goto 999
      if (bgt .lt. 0.0) goto 999
C
C---- Set up the format of the data collection lines and count number of
C---- items not including h,k,l,psi and F
C
      fmt = ' '
      fmt = '(3i4,F4.0,i2'
      lfnorm = 0
C---- Variable scan width?
      if (swdb .ne. 0.0) then
        fmt = fmt(1:index(fmt,' ')-1)//',f6.2'
        lfnorm = lfnorm + 1
        ptswd  = lfnorm
      else
        ptswd  = 7
        values(ptswd) = swda
      endif
C---- Variable number of scans?
      if (mxn .ne. 1) then
        fmt = fmt(1:index(fmt,' ')-1)//',f4.0'
        lfnorm = lfnorm + 1
        ptnsc  = lfnorm
      else
        ptnsc  = 8
        values(ptnsc) = 1.0
      endif
C---- Background counting time on file?
      if (bmo .eq. 2) then
        fmt = fmt(1:index(fmt,' ')-1)//',f6.0'
        lfnorm = lfnorm + 1
        ptbgt  = lfnorm
      elseif (bmo .eq. 1) then
        ptbgt = 9
C---- You don't know the length of the background counting time for bmo 1
C---- at this point in time.
      elseif (bmo .eq. 0) then
        ptbgt = 9
        values(ptbgt) = bgt
      endif
C---- Background 1
      lfnorm = lfnorm + 1
      ptcb1  = lfnorm
C---- Peak with background
      lfnorm = lfnorm + 1
      ptcsum = lfnorm
C---- Background 2
      lfnorm = lfnorm + 1
      ptcb2  = lfnorm
C---- Finish off the format statement
      fmt = fmt(1:index(fmt,' ')-1)//',3f8.0)'
C
C---- Set up the wavelength for output
C
      flam1 = wavel
      flam2 = 0.0
      flam3 = 0.0
      weigh1 = 1.0
      weigh2 = 0.0
      weigh3 = 0.0
      do 201 i = 1, navtab
        if (abs(wavtab(1,i)-wavel).le.0.005) then
          wavel  = wavtab(1,i)
          flam1  = wavtab(2,i)
          flam2  = wavtab(3,i)
          weigh1 = wavtab(4,i)
          weigh2 = wavtab(5,i)
        endif
 201  continue
        if      (cdefwv.eq.1) then
          cwave = flam1
        else if (cdefwv.eq.2) then
          cwave = flam2
        else if (i.ne.0) then
          cwave = wavel
        endif
C
C---- Set up scan type for output: mstnn points into meas
C
      mstnn = 4
      if (sct .lt.  1.0) mstnn = 2
      if (sct. gt. 89.0) mstnn = 1
C---- Time interval in seconds between reference reflections
      sinter = 60.*sinter
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
C---- Scale cell dimensions for wavelength
      do 403 i = 1, 3
       cellub(i) = wavel*cellub(i)
 403  continue
      volub=abs(wavel*wavel*wavel/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.
      temp3 = (cellub(1)+cellub(2)+cellub(3))/3.0
      if((abs(cellub(1)-temp3).lt.0.002*temp3).and.
     1   (abs(cellub(2)-temp3).lt.0.002*temp3).and.
     2   (abs(cellub(3)-temp3).lt.0.002*temp3)) then
        cellub(1) = temp3
        cellub(2) = cellub(1)
        cellub(3) = cellub(1)
      else
        temp2 = 0.5*(cellub(1)+cellub(2))
        if(abs(cellub(1)-cellub(2)).lt.0.002*temp2) then
          cellub(1) = temp2
          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 h k ls*****'')')
      endif
C
C---- Set up fixed values of reference reflections
C
      do 1005 iz = 1 , numref
C---- Get the crystal-based azimuthal angle
        ih = href(iz)
        ik = kref(iz)
        il = lref(iz)
        psi = 0.0
        call psic
        psirf(iz)=psicry
        if(trans) then
C---- Transform the stored reference hkl's
          href(iz) = -href(iz)
          kref(iz) = -kref(iz)
          lref(iz) = -lref(iz)
        endif
C---- Get sin theta needed for scan transition angle and swidth
        if (abs(theta) .le. sct) then
C---- omega scans
          scnrf(iz) = 0.0
        else
C---- omega - 2-theta scans
          scnrf(iz) = 2.0
        endif
        swdrf(iz) = 0.02*real(int(
     1   50.0*(swda + swdb* sint / sqrt(1.0 - sint*sint))))
        cnprf(iz) = 1.0/speref(iz)
        cnbrf(iz) = swdrf(iz)/speref(iz)
C---- Keep on eye on the filter codes that have been used
        oka = .false.
        do 2010 i = 1, nfilt
          if (fltref(iz) .eq. nfltcd(i)) oka = .true.
 2010   continue
        if (.not. oka) then
          nfilt = nfilt + 1
          nfltcd(nfilt) = fltref(iz)
        endif
 1005 continue
C
C---- Set Firstr to indicate that INITLZ does not need calling again.
      firstr = .false.
C
      return
 999  write(*,'('' *****Illegal par values read from diffractometer'',
     1 '' file. Interpretation stopped.***** '')')
      stop
      end
C
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Philips PW1100 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, In the Press. 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 PW1100 corresponds to case (c) of section 3) "Dealing with
C     Manufacturers" design of DSHDF since 2-theta, omega, phi and kappa
C     are not found on the diffractometer output file but the current UB is
C     available (at least in the set up in Geneva where the output to the
C     teletype is captured on the hard disk of an IBM-PC look alike instead
C     of the Philips paper tape or magnetic tape output. The PW1100 is a
C     DSHDF type II diffractometer and hence according to DSHDF equation (11):
C
C                   ( -1   0   0)
C             W =   (  0   1   0)
C                   (  0   0  -1).
C
C                (  1          0           0          )
C     X(chi)  =  (  0          cos(chi)   -sin(chi)   )
C                (  0          sin(chi)    cos(chi)   )
C
C                (  cos(phi)  -sin(phi)    0          )
C     Z(phi)  =  (  sin(phi)   cos(phi)    0          )
C                (  0          0           1          )
C
C                (  cos(omega)  -sin(omega)    0      )
C     Z(omega) = (  sin(omega)   cos(omega)    0      )
C                (  0            0             1      )
C
      integer  i,j
      real x(3), r(3), s(3), d, d2, p, p2, temp
      real sinchb, coschb, sinphb, cosphb
      real sinom, cosom, sinchi, coschi, sinphi, cosphi
      real sindph, cosdph, sinpsi, cospsi, cost
      include 'common.tlb'
C---- Calculate the angles for bisecting position
      do 205 i = 1 ,3
        x(i) = ub(i,1)*real(ih)+ub(i,2)*real(ik)+ub(i,3)*real(il)
 205  continue
C---- Get p and d
      p2 = x(1)*x(1)+x(2)*x(2)
      d2 = p2 + x(3)*x(3)
      p  = sqrt(p2)
      d  = sqrt(d2)
C---- Get sint and cost, theta and two theta
      sint = 0.5*d
      if( sint.gt.1.0 ) then
        write(*,'('' Sint > 1.0 for reflection '',3i5)') ih, ik, il
        return
      endif
      cost = sqrt(1.0 - 0.25*d2)
      theta = rad*atan2( sint, cost)
      twoth = 2.0*theta
C==== Bisecting Position Angle Calculation
      sinchb = -x(3)/d
      coschb =     p/d
      if ( p .eq. 0.0)  then
C---- Set phi to 45 degrees if both X and Y are ZERO
        sinphb = 0.707107
        cosphb = 0.707107
      else
        sinphb = x(1)/p
        cosphb = x(2)/p
      endif
C---- Is the reflection set out of the bisecting position
      if (psi.ne.0.0) then
        sinpsi = sin(dar*psi)
        cospsi = cos(dar*psi)
C     Make omega. Temp will contain tan(omega). Force cos(omega) positive
C     so abs(omega)=< 90 degrees.
        temp = sinpsi*coschb/sinchb
        cosom = 1.0/sqrt(1.0+temp*temp)
        sinom = cosom*temp
C     Make chi. Temp will contain scale factor to ensure that
C     sinchb and sinchi have the same sign
        sinchi = sinchb
        coschi = cosom*coschb*cospsi
        temp = 1.0 / sqrt (sinchi*sinchi + coschi*coschi)
        sinchi = sinchi * temp
        coschi = coschi * temp
C     Make delta phi (dph). Temp will contain scale to ensure that
C     sindph and sinom have opposite signs
        sindph = -sinom
        cosdph =  cosom*coschi
        temp = 1.0 / sqrt (sindph*sindph + cosdph*cosdph)
        sindph = sindph * temp
        cosdph = cosdph * temp
C     Make phi from phi(bisecting) and delta phi
        sinphi = sinphb*cosdph + cosphb*sindph
        cosphi = cosphb*cosdph - sinphb*sindph
      else
C---- Bisecting position
        sinom = 0.0
        cosom = 1.0
        sinchi = sinchb
        coschi = coschb
        sinphi = sinphb
        cosphi = cosphb
      endif
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)
      r(1) =  s(1)*cosphi - s(2)*sinphi
      r(2) =  s(1)*sinphi + s(2)*cosphi
      r(3) =  s(3)
C---- *x(-chi)
      s(1) =  r(1)
      s(2) =  r(2)*coschi - r(3)*sinchi
      s(3) =  r(2)*sinchi + r(3)*coschi
C---- *z(-omega)
      r(1) =  s(1)*cosom - s(2)*sinom
      r(2) =  s(1)*sinom + s(2)*cosom
      r(3) =  s(3)
C---- *W (DSHDF Type II diffractometer)
      r(1) = -r(1)
      r(3) = -r(3)
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==== omega = rad*ATAN2(sinom,cosom) + theta
C==== chi   = rad*ATAN2(sinchi,coschi)
C==== phi   = rad*ATAN2(sinphi,cosphi)
C==== WRITE(12,'(1x,3i5,6f10.2)')
C====1 ih,ik,il, omega, twoth, chi, phi, psi, psicry
C
      return
      end



------cut here------------------------------------------------------------------


C
C==== COMMON Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 7th October 1992
C
C---- SPECIFICATION (DIMENSION,REAL,INTEGER,CHARACTER, ETC) STATATEMENTS
C
      integer scfsfl, nseq, ndtfil, datafl, parmfl, mxn, smo, bmo
      integer numref, href(10), kref(10), lref(10), fltref(10), lfnorm
      integer ptswd, ptnsc, ptbgt, ptcb1, ptcsum, ptcb2, ih, ik, il
      integer ptfilt, hface(50), kface(50), lface(50), nface, intin(10)
      integer nmcell, ipchar, zcell, natmtp, nfilt, nfltcd(20)
      integer ptscl, ptstd, nrefl, parclp, jostat
      integer profil(17:112), nstep, iscan, backsn
      integer hmin, hmax, kmin, kmax, lmin, lmax
      integer stdseq, sdtfil
      integer numang, hang(100), kang(100), lang(100), wang(100)
      integer navtab, mdlmz, mstnn, addrfl, cdefwv
      integer nstepw, nfridl, apc
      integer id1, id2, numstd, qaofl, qanfl, nqaocd
      integer slidc4
C
      real spe, swda, swdb, sct, bgt, speref(10), ub(3,3), dethz, detvt
      real wavel, refint(10), values(9), psi, dface(50), deface(50)
      real realin(10), cell(6), ecell(6), volum, evolum, cwave, ctemp
      real ctthn, ctthx, densit, edensi, tdensi, temint, presnt
      real fnamtp(24), peratm(24), eperat(24), forwg, eforwg
      real attfac(20), ettfac(20), fmnhed, fmn2th, ratpol, ertpol
      real bedwv, fmhomw, fdeadt, efdead, dar, rad, backl, backr
      real psicry, peak, time, cnp, cnb, cnl, cnr, swidth, scnrat
      real x1, x2, x3, sint, psirf(10), cnprf(10), cnbrf(10)
      real swdrf(10), scnrf(10), ubinv(3,3), bastim, prvtim
      real radius, slit, apmin, apmax, apta, aptb, flam1, flam2, flam3
      real sintl, cmn2th, cmix, weigh1, weigh2, weigh3
      real twoth, theta, omega, phi, chi, kappa, range, rate, firl, sir
      real beta, stlmin, stlmax, avint(10), sdint(10), count(10)
      real sinalp, cosalp, attn, gm(3,3), cellub(6), volub
      real twoang(100), wavtab(5,3)
      real csigma, timl, timp, timr
      real mptemp, sinter, ratio
      real stepwd, intnet, signet, stptim
C
      character title*66, compid*6, line*80, digit(0:11)*1
      character dsk*4, fmt*40, crytyp*1, inchar*1, namtyp(24)*2
      character charin(10)*15, alphc(26)*1, alphl(26)*1, lin2*80
      character monoch*1, code6*6, includ*8, exclud*8
      character rmode*2, extn*4, difnam*1, datxtn*8
      character meas(4)*22, device*40, softwr*40
      character family*30, addres(5)*40, email*40, fax*40, phone*40
      character cpdsou*76, namecm*76, nameml*76
      character cryscl*40, crysds*76, specdt*76
      character chemfm*140, addnam*32, qline*80
C
      logical ok, firstr, mdensi,  lform(5), trans, hklopn, condon
      logical reoren, ask(16), statbk, lquan, hspeed
C
C---- COMMON DECLARATIONS
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
      common/difrin/scfsfl, nseq, ndtfil, datafl, parmfl, mxn, smo, bmo
      common/difrin/numref, href, kref, lref, fltref, lfnorm, ih, ik, il
      common/difrin/ptswd, ptnsc, ptbgt, ptcb1, ptcsum, ptcb2
      common/difrin/ptfilt, hface, kface, lface, nface, intin, nmcell
      common/difrin/ipchar, zcell, natmtp, nfilt, nfltcd, ptscl, ptstd
      common/difrin/nrefl, jostat, profil, iscan, nstep, parclp, backsn
      common/difrin/hmin, hmax, kmin, kmax, lmin, lmax
      common/difrin/stdseq, sdtfil
      common/difrin/numang, hang, kang, lang, wang
      common/difrin/navtab, mdlmz, mstnn, addrfl, cdefwv
      common/difrin/nstepw, nfridl, apc
      common/difrin/id1,id2, numstd, qaofl, qanfl, nqaocd
      common/difrin/slidc4
C
      common/difrrl/spe, swda, swdb, sct, bgt, speref, ub, dethz, detvt
      common/difrrl/wavel, refint, values, psi, dface, deface, realin
      common/difrrl/cell, ecell, volum, evolum, cwave, ctemp, ctthn
      common/difrrl/ctthx, densit, edensi, tdensi, temint, presnt
      common/difrrl/fnamtp, peratm, eperat, forwg, eforwg
      common/difrrl/attfac, ettfac, fmnhed, fmn2th, ratpol, ertpol
      common/difrrl/bedwv, fmhomw, fdeadt, efdead, dar, rad
      common/difrrl/backl, backr, psicry, peak, time, cnp, cnb, cnl, cnr
      common/difrrl/swidth, scnrat, x1 ,x2, x3, sint, firl, sir
      common/difrrl/psirf, cnprf, cnbrf, swdrf, scnrf, ubinv
      common/difrrl/bastim, prvtim, radius, slit, apmin, apmax
      common/difrrl/apta, aptb, flam1, flam2, flam3, sintl
      common/difrrl/cmn2th, cmix, weigh1, weigh2, weigh3
      common/difrrl/twoth, theta, omega, chi, phi, kappa, range, rate
      common/difrrl/beta, stlmin, stlmax, avint, sdint, count
      common/difrrl/sinalp, cosalp, attn, gm, cellub, volub
      common/difrrl/twoang, wavtab
      common/difrrl/csigma, timl, timp, timr
      common/difrrl/mptemp, sinter, ratio
      common/difrrl/stepwd,intnet,signet,stptim
C
      common/difrch/title, compid, line, digit, dsk, fmt, crytyp
      common/difrch/inchar, namtyp, charin, alphc, alphl, lin2
      common/difrch/monoch, code6, includ, exclud, rmode, extn
      common/difrch/difnam, datxtn, meas, device, softwr
      common/difrch/family, addres, email, fax, phone
      common/difrch/cpdsou, namecm, nameml
      common/difrch/cryscl, crysds, specdt, chemfm
      common/difrch/addnam, qline
C
      common/difrlg/ok, firstr, mdensi, lform, trans, hklopn, condon
      common/difrlg/reoren, ask, statbk, lquan, hspeed



------cut here------------------------------------------------------------------


Bernardinelli, G.
Laboratoire de Cristallographie         
University of Geneva                    
24, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
bernard@sc2a.unige.ch
[+[41] 22] 781 21 92
[+[41] 22] 702 62 02

Bonhomme, F.
Laboratoire de Cristallographie         
University of Geneva                    
24, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
bonhomme@sc2a.unige.ch
[+[41] 22] 781 21 92
[+[41] 22] 702 63 72

Flack, H.D.                   
Laboratoire de Cristallographie         
University of Geneva                    
24, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
flack@sc2a.unige.ch
[+[41] 22] 781 21 92
[+[41] 22] 702 62 49

Francois, M.
Laboratoire de Cristallographie         
University of Geneva                    
24, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
francois@sc2a.unige.ch
[+[41] 22] 781 21 92
[+[41] 22] 702 67 50

Gladyshevskii, R.E.
Laboratoire de Cristallographie
University of Geneva                    
24, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
gladych@sc2a.unige.ch
[+[41] 22] 781 21 92
[+[41] 22] 702 62 36

Kubel, F.
Department of Inorganic, Analytical and Applied Chemistry
University of Geneva                    
16, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
kubel@sc2a.unige.ch
[+[41] 22] 781 ?? ??
[+[41] 22] 702 60 34

Mao, S.-Y.
Department of Inorganic, Analytical and Applied Chemistry
University of Geneva                    
16, quai Ernest-Ansermet                 
CH-1211 Geneva 4                        
Switzerland                             
mao@sc2a.unige.ch
[+[41] 22] 781 ?? ??
[+[41] 22] 702 60 34


------cut here------------------------------------------------------------------
