*     ****************
      SUBROUTINE INPUT
*     ****************
*
*     This subroutine reads and prints the input data and initializes
*     other data.
*
      CHARACTER ATEXT*4,NAME(10)*2,TEXT*80
      REAL F(21,10),FLAG(6),S(21,10)
      INTEGER IFLAG(6)
      INCLUDE 'COMMONS.FOR'
      COMMON/CHARS/ATEXT(30),TEXT
      COMMON/SPACE/JUMP,A,B,C,COSALF,COSBET,COSGAM,SINALF,SINBET,SINGAM
      COMMON/IO/EPS,RELAX,RELAXB,RELAXC,SUMOBS,SUMW,SUMWBS,SW,IPONS,
     *          LMAT,LOUT,NUM
      COMMON/DATA/CUT,NREFL,NH(3),IALFA,ISTAP,NCHAN,ICODE(999),
     *       HKL(3,999),SJJ(999),YOBS(15000),WEIGHT(15000),IORD(15000)
      COMMON/PARAMS/AINSTR(10),ACRYST(22), ATOM(14,30),
     *              LINSTR(10),LCRYST(22),LATOM(14,30),
     *              PINSTR(10),PCRYST(22),PATOM(14,30)
      COMMON/CONSTR/ML,LN(4),COEF(4,10),IL(4,9),JL(4,9),IQ(4,9),JQ(4,9)
      SAVE
      DATA RAD/1.7453292520E-2/

      CALL TIMES(1)
*
*     Read and print problem title. Subroutine TIMES gives the
*     elapsed CPU-time.
*
      WRITE(6,1001)
      READ(5,'(A)') TEXT
      CALL TIMES(2)
      WRITE(6,1002) TEXT
*
*     Find out the space group of crystal structure.
*     Read incident angle to specify measurement geometry
*     and define which radiation was used, neutrons or X-rays.
*
      CALL GROUP(ANGLE,XRAYS)
*
*     Read and print problem parameters.
*
      CALL READER(5,TEXT,0,FLAG,1,IFLAG,2,.TRUE.)
      EPS=FLAG(1)
      K=IFLAG(1)
      LIM=IFLAG(2)
      CALL READER(5,TEXT,0,FLAG,0,IFLAG,6,.TRUE.)
      ICENT=IFLAG(1)
      IRL=IFLAG(2)
      KL=IFLAG(3)
      KM=IFLAG(4)
      N=IFLAG(5)
      MAGROT=IFLAG(6)
      CALL READER(5,TEXT,0,FLAG,1,IFLAG,0,.FALSE.)
      SLABDA=FLAG(1)
      IRL=IRL+1
      WRITE(6,1003) EPS
      HALPER=.FALSE.
      UNIAX=.FALSE.
      CUBIC=.FALSE.
      IF(K.EQ.1) THEN
         WRITE(6,1004)
         HALPER=.TRUE.
      ELSE IF(K.EQ.2) THEN
         WRITE(6,1005)
         UNIAX=.TRUE.
      ELSE IF(K.EQ.3) THEN
         WRITE(6,1006)
         CUBIC=.TRUE.
      END IF
*
*     At this point SLABDA equals to wavelength in Angstroms.
*
      WRITE(6,1007) SLABDA
*
*     Read equivalent position vector and matrices.
*
      DO 30 I=1,3
      EQUIV(4,I,1)=0
      DO 30 J=1,3
      IJ=(I-1)*3+J
      IF(I.EQ.J) THEN
         EQUIV(J,I,1)=1
      ELSE
         EQUIV(J,I,1)=0
      END IF
      IF(I.EQ.J.OR.MAGROT.GT.0) THEN
         DO 20 K=1,MAGROT
         IF(I.EQ.J) RM(IJ,K,1)=1
         IF(I.NE.J) RM(IJ,K,1)=0
   20    CONTINUE
      END IF
   30 CONTINUE
*                                   Line 5 and 6.
      DO 50 IR=2,IRL
      CALL READER(5,TEXT,0,EQUIV(1,1,IR),12,IFLAG,0,.TRUE.)
      DO 40 K=1,MAGROT
   40 CALL READER(5,TEXT,0,RM(1,K,IR),9,IFLAG,0,.TRUE.)
   50 CONTINUE
*
*     Read scattering lengths and dispersion corrections, if any.
*     Lines 7-9.
*
      IF(KL.GT.0) CALL READER(5,TEXT,0,SCAT,KL,IFLAG,0,.TRUE.)
      IF(KM.GT.0 .AND. XRAYS) THEN
         CALL READER(5,TEXT,0,DISPRE,KM,IFLAG,0,.TRUE.)
         CALL READER(5,TEXT,0,DISPIM,KM,IFLAG,0,.TRUE.)
      END IF
*
*     Read atom names and scattering curves.
*
      DO 60 I=1,KM
      CALL READER(5,NAME(I),2,FLAG,0,IFLAG,0,.TRUE.)
      TYPE=.TRUE.
      DO 60 J=1,21
      IF(TYPE) THEN
         CALL READER(5,TEXT,0,S(J,I),1,IFLAG,0,.TRUE.)
         CALL READER(5,TEXT,0,F(J,I),1,IFLAG,0,.FALSE.)
         IF(S(J,I).EQ.-100) THEN
            TYPE=.FALSE.
            S(J,I)=S(J-1,I)
            X=S(J,I)
            F(J,I)=F(J-1,I)
            Y=F(J,I)
         END IF
      ELSE
         S(J,I)=X
         F(J,I)=Y
      END IF
   60 CONTINUE
*
*     Read number of cycles, relaxation factors and output parameters.
*     Line 12.
*
      CALL READER(5,TEXT,0,FLAG,0,IFLAG,1,.TRUE.)
      NCYCLE=IFLAG(1)
      CALL READER(5,TEXT,0,FLAG,3,IFLAG,3,.FALSE.)
      RELAXC=FLAG(1)
      RELAXB=FLAG(2)
      RELAX=FLAG(3)
      LOUT=IFLAG(1)
      IPONS=IFLAG(2)
      LMAT=IFLAG(3)
*
*     From this on SLABDA equals to wavelength**2/4.
*
      SLABDA=SLABDA*SLABDA/4
*
*     Print all input data.
*
      WRITE(6,1008) NCYCLE
      IF(LOUT.EQ.1) WRITE(6,1009)
      IF(LMAT.EQ.1) WRITE(6,1010)
      IF((IPONS.EQ.1).OR.(IPONS.EQ.3)) WRITE(6,1011)
      IF((IPONS.EQ.2).OR.(IPONS.EQ.3)) WRITE(6,1012)
      WRITE(6,1013) RELAXC,RELAXB,RELAX
      IF(ICENT.EQ.2) THEN
         WRITE(6,1014)
      ELSE
         WRITE(6,1015)
      END IF
      WRITE(6,1016)
      DO 80 IR=2,IRL
      WRITE(6,1017) ((EQUIV(J,I,IR),J=1,4),I=1,3)
      DO 70 K=1,MAGROT
   70 WRITE(6,1018) (RM(IJ,K,IR),IJ=1,9)
   80 CONTINUE
      IF(KL.GT.0) WRITE(6,1020) (I,SCAT(I),I=1,KL)
      IF(KM.GT.0) WRITE(6,1021)
      DO 100 I=1,KM
      WRITE(6,1022) NAME(I)
      DO 90 J=1,20
      WRITE(6,1023) S(J,I),F(J,I)
      IF((S(J,I).EQ.S(J+1,I)).AND.(F(J,I).EQ.F(J+1,I))) GO TO 100
   90 CONTINUE
  100 CONTINUE
      WRITE(6,1024)
*
*     Read scale and temp factors.
*
      CALL READER(5,TEXT,0,FLAG,2,IFLAG,0,.TRUE.)
      SCALE=FLAG(1)
      PCRYST(7)=FLAG(2)
      PINSTR(2)=1/SCALE
*
*     Read and print data for each atom.
*
      DO 120 I=1,N
      CALL READER(5,ATEXT(I),4,FLAG,0,IFLAG,3,.TRUE.)
      NTYP(I)=IFLAG(1)
      MTYP(I)=IFLAG(2)
      MEQ(I)=IFLAG(3)
      CALL READER(5,TEXT,0,PATOM(1,I),8,IFLAG,0,.FALSE.)
      CALL READER(5,TEXT,0,PATOM(9,I),6,IFLAG,0,.TRUE.)
      WRITE(6,1025) ATEXT(I),NTYP(I),MTYP(I),MEQ(I),(PATOM(J,I),J=1,14)
  120 CONTINUE
      WRITE(6,1026) SCALE,PCRYST(7)
*
*     Read and print halfwidth and zeropoint parameters.
*
      CALL READER(5,TEXT,0,PINSTR(3),8,IFLAG,0,.TRUE.)
      WRITE(6,1027) PINSTR(9),PINSTR(10)
      WRITE(6,1028) (PINSTR(I),I=3,8)
*
*     IFLAG(1) .GE. 0  <===>  a, b, c, alpha, beta, gamma are given.
*     IFLAG(1) .LT. 0  <===>  The coefficients of products of Miller
*                             indices are given.
*     In former case the conversion to latter case is performed.
*
      CALL READER(5,TEXT,0,FLAG,0,IFLAG,1,.TRUE.)
      CALL READER(5,TEXT,0,PCRYST,6,IFLAG,0,.TRUE.)
      IF(IFLAG(1) .LT. 0) THEN
         CALL PROFIL
      ELSE
         A=PCRYST(1)
         B=PCRYST(2)
         C=PCRYST(3)
         COSALF=COS(RAD*PCRYST(4))
         COSBET=COS(RAD*PCRYST(5))
         COSGAM=COS(RAD*PCRYST(6))
         SINALF=SQRT(1-COSALF**2)
         SINBET=SQRT(1-COSBET**2)
         SINGAM=SQRT(1-COSGAM**2)
         VOL=A*B*C*SQRT(1-COSALF*COSALF-COSBET*COSBET-COSGAM*COSGAM+
     *       2*COSALF*COSBET*COSGAM)
         AS=B*C*SINALF/VOL
         BS=C*A*SINBET/VOL
         CS=A*B*SINGAM/VOL
         COSAS=(COSBET*COSGAM-COSALF)/(SINBET*SINGAM)
         COSBS=(COSGAM*COSALF-COSBET)/(SINGAM*SINALF)
         COSCS=(COSALF*COSBET-COSGAM)/(SINALF*SINBET)
         PCRYST(1)=AS*AS
         PCRYST(2)=BS*BS
         PCRYST(3)=CS*CS
         PCRYST(4)=2*BS*CS*COSAS
         PCRYST(5)=2*AS*CS*COSBS
         PCRYST(6)=2*AS*BS*COSCS
*
*     CONVERT B(I,J) TO BETA(I,J)
*
         DO 130 I=1,N
         PATOM( 9,I)=PATOM( 9,I)*AS*AS/4
         PATOM(10,I)=PATOM(10,I)*BS*BS/4
         PATOM(11,I)=PATOM(11,I)*CS*CS/4
         PATOM(12,I)=PATOM(12,I)*AS*BS/4
         PATOM(13,I)=PATOM(13,I)*AS*CS/4
  130    PATOM(14,I)=PATOM(14,I)*BS*CS/4
      END IF
      WRITE(6,1029) (PCRYST(I),I=1,6)
*
*     Get preparation data.
*
      CALL PREPAR
      CALL GEOMET(ANGLE)
      J=0
      SW=0
      SUMOBS=0
      SUMWBS=0
      DO 10 I=1,NCHAN
      IF(ABS(YOBS(I)).LT.0.001) THEN
         J=J+1
         YOBS(I)=0
      ELSE
         SW=SW+WEIGHT(I)
         SUMOBS=SUMOBS+ABS(YOBS(I))
         SUMWBS=SUMWBS+WEIGHT(I)*YOBS(I)*YOBS(I)
      END IF
   10 CONTINUE
      NUM=NCHAN-J
*
*     Read and write asymmetry parameter and texture parameters.
*
      CALL READER(5,TEXT,0,PINSTR,1,IFLAG,0,.TRUE.)
      CALL READER(5,TEXT,0,PCRYST(8),15,IFLAG,0,.TRUE.)
      WRITE(6,1030) PINSTR(1),(PCRYST(J),J=8,22)
      DO 160 I=1,N
      DO 160 J=1,14
      ATOM(J,I)=0
  160 LATOM(J,I)=0
      WRITE(6,1031) NREFL,CUT
      MAXCON=0
*
*     Read the number of refinement parameters.
*
      CALL READER(5,TEXT,0,FLAG,0,IFLAG,1,.TRUE.)
      MAXS=IFLAG(1)
      IF(MAXS.EQ.0) GO TO 230
      WRITE(6,1032)
*
*     Read and print codewords for atomic parameters.
*
      DO 180 I=1,N
      CALL READER(5,TEXT,0,ATOM(1,I),8,IFLAG,0,.TRUE.)
      CALL READER(5,TEXT,0,ATOM(9,I),6,IFLAG,0,.TRUE.)
      WRITE(6,1033) ATEXT(I),(ATOM(J,I),J=1,14)
      DO 170 J=1,14
  170 ATOM(J,I)=CODING(ATOM(J,I),LATOM(J,I))
  180 CONTINUE
*
*     Read and print codewords
*        - for scale, temperature factor
*        - zero point and halfwidth parameters
*        - for cell constants
*        - for asymmetry parameter and texture parameters.
*
      CALL READER(5,TEXT,0,AINSTR(2),1,IFLAG,0,.TRUE.)
      CALL READER(5,TEXT,0,ACRYST(7),1,IFLAG,0,.FALSE.)
      WRITE(6,1034) AINSTR(2)
      WRITE(6,1035) ACRYST(7)
      CALL READER(5,TEXT,0,AINSTR(3),8,IFLAG,0,.TRUE.)
      WRITE(6,1036) AINSTR(9),AINSTR(10)
      WRITE(6,1037) (AINSTR(J),J=3,8)
      CALL READER(5,TEXT,0,ACRYST,6,IFLAG,0,.TRUE.)
      WRITE(6,1038) (ACRYST(J),J=1,6)
      CALL READER(5,TEXT,0,AINSTR,1,IFLAG,0,.TRUE.)
      CALL READER(5,TEXT,0,ACRYST(8),15,IFLAG,0,.TRUE.)
      WRITE(6,1039) AINSTR(1)
      WRITE(6,1040) (ACRYST(J),J=8,22)
      DO 190 I=1,10
  190 AINSTR(I)=CODING(AINSTR(I),LINSTR(I))
      DO 200 I=1,22
  200 ACRYST(I)=CODING(ACRYST(I),LCRYST(I))
*
*     STEADY = TRUE means cell constants to be varied. So COORD and
*     PROFIL must be called at the beginning of each cycle (see ITER).
*
      STEADY = LCRYST(1) + LCRYST(2) + LCRYST(3) +
     *         LCRYST(4) + LCRYST(5) + LCRYST(6) .EQ. 0
*
*     Read constraints if any.
*
      CALL READER(5,TEXT,0,FLAG,0,IFLAG,2,.TRUE.)
      ML=IFLAG(1)
      MQ=IFLAG(2)
      MAXCON=ML+MQ
      K=0
      DO 220 I=1,MAXCON
  210 K=K+1
*
*     Read coefficient and indices for constraint functions.
*
      CALL READER(5,TEXT,0,COEF(I,K),1,IFLAG,4,.TRUE.)
*
*                             Linear or quadratic indices.
      IF(IFLAG(1).GT.0) THEN
         IL(I,K)=IFLAG(1)
         JL(I,K)=IFLAG(2)
*                             Quadratic indices.
         IQ(I,K)=IFLAG(3)
         JQ(I,K)=IFLAG(4)
         IF(K.LT.10) GO TO 210
      END IF
*                             Last term.
      COEF(I,10)=COEF(I,K)
      LN(I)=K-1
      K=0
  220 CONTINUE
  230 WRITE(6,1041) MAXS
      IF(MAXCON.EQ.0) THEN
         WRITE(6,1042)
         GO TO 260
      END IF
      WRITE(6,1043) MAXCON
      DO 240 I=1,ML
      LNI=LN(I)
      WRITE(6,1044) (COEF(I,J),IL(I,J),JL(I,J),J=1,LNI)
  240 WRITE(6,1045) COEF(I,10)
      DO 250 K=1,MQ
      I=K+ML
      LNI=LN(I)
      WRITE(6,1046) (COEF(I,J),IL(I,J),JL(I,J),IQ(I,J),JQ(I,J),J=1,LNI)
  250 WRITE(6,1045) COEF(I,10)
*
*     Calculate for each reflection alpha**2.
*
  260 DO 310 I=1,NREFL
      DO 280 J=1,3
  280 HKL(J,I)=HKL(J,I)/NH(J)
      X=PCRYST(1)*HKL(1,I)*HKL(1,I)+PCRYST(2)*HKL(2,I)*HKL(2,I)+
     *  PCRYST(3)*HKL(3,I)*HKL(3,I)+PCRYST(4)*HKL(2,I)*HKL(3,I)+
     *  PCRYST(5)*HKL(1,I)*HKL(3,I)+PCRYST(6)*HKL(1,I)*HKL(2,I)
      IF(UNIAX.AND.(ICODE(I).NE.1)) THEN
         Y=(PCRYST(3)*HKL(3,I)+PCRYST(4)*HKL(2,I)/2+
     *                         PCRYST(5)*HKL(1,I)/2)/NH(3)**2
         FAC(I)=(Y*NH(3))**2/(X*PCRYST(3))
      END IF
*
*     COORD is called define cartesian coordinate orientation.
*
      CALL COORD(HKL(1,I),I)
*
      X=SQRT(X)/2
      DO 300 J=1,KM
      DO 290 K=1,20
      K1=MAX(K-1,1)
      IF(X.GT.S(K,J)) GO TO 290
      FF(I,J)=(S(K1,J)-X)/(S(K1,J)-S(K,J))*(F(K,J)-F(K1,J))+F(K1,J)
      GO TO 300
  290 CONTINUE
      FF(I,J)=F(20,J)
  300 CONTINUE
  310 CONTINUE
      CALL TIMES(1)
      RETURN
 1001 FORMAT('1NUCLEAR AND MAGNETIC STRUCTURE REFINEMENT PROGRAM FOR ',
     1       'NEUTRON AND X-RAY POWDER PROFILES - FORTRAN 77 VERSION ',
     2       '- OCTOBER 1987')
 1002 FORMAT('0******    ',A80,' ******',/)
 1003 FORMAT(' Convergence criteria: EPS =',F6.1)
 1004 FORMAT(' MAGNETIC INTENSITIES CALCULATION ACCORDING TO HALPERN',
     1       ' AND JOHNSON.')
 1005 FORMAT(' UNIAXIAL CONFIGURATIONAL SPIN SYMMETRY.')
 1006 FORMAT(' CUBIC CONFIGURATIONAL SPIN SYMMETRY.')
 1007 FORMAT('0WAVELENGTH =',F11.7)
 1008 FORMAT(' NUMBER OF CYCLES=',I6)
 1009 FORMAT(' PRINTED OUTPUT OF OBS + CALC INTENSITIES ON LAST CYCLE')
 1010 FORMAT(' PRINTED OUTPUT OF CORRELATION  -  MATRIX ON LAST CYCLE')
 1011 FORMAT(' PUNCHED OUTPUT OF STRUCTURE FACTOR TAPE  ON LAST CYCLE')
 1012 FORMAT(' PUNCHED OUTPUT OF OBS + CALC INTENSITIES ON LAST CYCLE')
 1013 FORMAT('0RELAXATION FACTORS',/,' FOR COORDINATES    =',F5.2,/,
     1      ' FOR TEMP FACTORS   =',F5.2,/,' FOR SCALING FACTORS=',F5.2)
 1014 FORMAT('0THE STRUCTURE IS CENTROSYMMETRIC.')
 1015 FORMAT('0THE STRUCTURE IS NON CENTROSYMMETRIC.')
 1016 FORMAT('0EQUIVALENT POSITIONS=')
 1017 FORMAT(5X,3F5.1,F13.2,/,:,' R=  ',3F5.1,4X,'T=',F7.2)
 1018 FORMAT(5X,3F5.1,/,' M=  ',3F5.1,/,5X,3F5.1,/)
 1020 FORMAT('0SCATTERING LENGTHS:',(T22,'B(',I1,') =',F7.3))
 1021 FORMAT('0FORMFACTORS')
 1022 FORMAT(/,10X,A2,/,' SINTH/LAMBDA      F')
 1023 FORMAT(1X,F8.4,F13.4)
 1024 FORMAT('0***INITIAL PARAMETERS***',/,' ATOM NTP MTP MRT',3X,'X',
     1        7X,'Y',7X,'Z',7X,'B',7X,'N',7X,'KX',6X,'KY',6X,'KZ',6X,
     2        'B11',5X,'B22',5X,'B33',5X,'B12',5X,'B13',5X,'B23')
 1025 FORMAT(1X,A4,I3,2I4,14F8.5)
 1026 FORMAT('0OVERALL SCALE FACTOR=',F12.6,/,' OVERALL TEMP. FACTOR=',
     *       F12.6)
 1027 FORMAT('0ZEROPOINT=',F23.6,F12.6)
 1028 FORMAT('0HALFWIDTH PARAMETERS:',/,' GAUSSIAN  =',10X,4F12.4,/,
     *       ' LORENTZIAN=',10X,2F12.4)
 1029 FORMAT('0CELL CONSTANTS=',6X,6F12.7)
 1030 FORMAT('0ASYMMETRY PARAMETER= ',F12.6,/,
     *       ' TEXTURE PARAMETERS = ',8F12.6,/,22X,8F12.6)
 1031 FORMAT('0For this analysis ',I4,' reflections',
     *       ' are used with CUT =',F5.1,'.') 
 1032 FORMAT('0***CODING OF VARIABLES***',/,' ATOM',16X,'X',7X,'Y',7X,
     1       'Z',7X,'B',7X,'N',6X,'KX',6X,'KY',6X,'KZ',6X,'B11',5X,
     2       'B22',5X,'B33',5X,'B12',5X,'B13',5X,'B23')
 1033 FORMAT(1X,A4,10X,14F8.3)
 1034 FORMAT('0OVERALL SCALE FACTOR=',F8.2)
 1035 FORMAT(' OVERALL TEMP. FACTOR=',F8.2)
 1036 FORMAT('0ZEROPOINT=',F19.2,F8.2)
 1037 FORMAT('0HALFWIDTH PARAMETERS:',/,' GAUSSIAN  =',10X,4F8.2,/,
     *       ' LORENTZIAN=',10X,2F8.2)
 1038 FORMAT('0CELL CONSTANTS=',6X,6F8.2)
 1039 FORMAT('0ASYMMETRY PARAMETER=',F9.2)
 1040 FORMAT(' TEXTURE PARAMETERS = ',8F8.2,/,22X,8F8.2,/)
 1041 FORMAT(/,'0TOTAL NUMBER OF VARIABLES =',I3)
 1042 FORMAT(' NO CONSTRAINT FUNCTIONS.')
 1043 FORMAT(1X,I2,' CONSTRAINT FUNCTIONS =',/)
 1044 FORMAT(8(F6.2,' * X(',I1,',',I2,')'))
 1045 FORMAT(10X,'=',E12.5,/)
 1046 FORMAT(5(F6.2,' * X(',I1,',',I2,') * X(',I1,',',I2,')'))
      END
