C 
C One-purpose program TCROT to change the material parameters
C in the "Twisted Crystal" model in such a way, that the plane
C of polarization vectors in the new "Oblique Twisted Crystal" model
C is slanted.
C
C Version: 8.00
C Date: 2022, April 24
C
C Coded by Petr Bulant according to Ludek Klimes's program MODMOD
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Name of the file with the original model:
C     MODEL='string'... String containing the name of the input data
C             file specifying the original model.  The program is
C             intended to be used only for the "Twisted Crystal" model,
C             specification of MODEL='tc-mod.dat' is strongly
C             recommended.
C             Default: MODEL='tc-mod.dat'
C Rotation vector:
C     TCE1=real, TCE2=real, TCE3=real ... Three components of
C              the rotation vector describing the rotation of the plane
C              of polarization vectors in the model.
C              The plane is rotated along the axis given by the rotation
C              vector, the angle of the rotation in radians is given by
C              the Cartesian length of the vector.
C              Default: TCE1=TCE2=TCE3=0. (no rotation).
C Name of the output rotated model file:
C     MODOUT='string'... Name of the output file describing the new
C             rotated model.  File MODOUT is a copy of file MODEL,
C             with the functional values replaced by new ones.
C             Default: MODOUT='tco-mod.out'
C Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             
C             forms.for.
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             
C             forms.for.
C Form of the input/output files:
C     FORM='string'... Form of the input/output files: either
C             'FORMATTED' or 'UNFORMATTED' (case insensitive).
C             Default: FORM='FORMATTED'
C     FORMW='string'... Form of the output file. If both FORM and FORMW
C             are specified, FORMW is used for output files.
C             Default: FORMW=FORM
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C.......................................................................
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,MODEL1,NEWMOD,NEWVAL
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FILE1,FILE2
      INTEGER LU1,LU2
      PARAMETER (LU1=1,LU2=2)
C
C     Auxiliary storage locations:
      CHARACTER*251 LINE
      INTEGER NSRF,NCB,N,I1,I2
      CHARACTER*3 TSRF(1)
      REAL REA11,IMA11,REA21,IMA21,REA31,IMA31,
     *     REA12,IMA12,REA22,IMA22,REA32,IMA32,
     *     REA13,IMA13,REA23,IMA23,REA33,IMA33,REA(3,3),IMA(3,3),
     * REE11,IME11,REE21,IME21,REE31,IME31,
     * REE12,IME12,REE22,IME22,REE32,IME32,
     * REE13,IME13,REE23,IME23,REE33,IME33,REE(3,3),IME(3,3),
     * REB(3,3),IMB(3,3),REET(3,3),IMET(3,3),
     * AE,AE2,COSAE,SINAE,AUX1,AUX2,E1,E2,E3,VS
      EQUIVALENCE (REA11,REA(1,1)),(IMA11,IMA(1,1)),(REA21,REA(2,1)),
     *            (IMA21,IMA(2,1)),(REA31,REA(3,1)),(IMA31,IMA(3,1)),
     *            (REA12,REA(1,2)),(IMA12,IMA(1,2)),(REA22,REA(2,2)),
     *            (IMA22,IMA(2,2)),(REA32,REA(3,2)),(IMA32,IMA(3,2)),
     *            (REA13,REA(1,3)),(IMA13,IMA(1,3)),(REA23,REA(2,3)),
     *            (IMA23,IMA(2,3)),(REA33,REA(3,3)),(IMA33,IMA(3,3)),
     *            (REE11,REE(1,1)),(IME11,IME(1,1)),(REE21,REE(2,1)),
     *            (IME21,IME(2,1)),(REE31,REE(3,1)),(IME31,IME(3,1)),
     *            (REE12,REE(1,2)),(IME12,IME(1,2)),(REE22,REE(2,2)),
     *            (IME22,IME(2,2)),(REE32,REE(3,2)),(IME32,IME(3,2)),
     *            (REE13,REE(1,3)),(IME13,IME(1,3)),(REE23,REE(2,3)),
     *            (IME23,IME(2,3)),(REE33,REE(3,3)),(IME33,IME(3,3))
      DATA REE,IME/18*0./
      DATA TSRF/'   '/
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+TCROT: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
      IF(FILE1.EQ.' ') THEN
C       TCROT-01
        CALL ERROR('TCROT-01: No input file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      END IF
      WRITE(*,'(A)') '+TCROT: Working...            '
      CALL RSEP1(LU1,FILE1)
C
C     Checking input data MODEL:
      CALL RSEP3T('MODEL',FILE1,'tc-mod.dat')
C
C     Reading input data MODEL for the model to be updated:
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
C     Rotation matrix REE:
      CALL RSEP3R('TCE1',E1,0.)
      CALL RSEP3R('TCE2',E2,0.)
      CALL RSEP3R('TCE3',E3,0.)
      AE2=E1*E1+E2*E2+E3*E3
      IF (AE2.EQ.0.) THEN
        REE11=1.
        REE22=1.
        REE33=1.
      ELSE
        AE=SQRT(AE2)
        COSAE=COS(AE)
        AUX1=(1-COSAE)/AE2
        SINAE=SIN(AE)
        AUX2=SINAE/AE
        REE11=COSAE+E1*E1*AUX1
        REE21=      E2*E1*AUX1                -E3*AUX2
        REE31=      E3*E1*AUX1        +E2*AUX2
        REE12=      E1*E2*AUX1                +E3*AUX2
        REE22=COSAE+E2*E2*AUX1
        REE32=      E3*E2*AUX1-E1*AUX2
        REE13=      E1*E3*AUX1        -E2*AUX2
        REE23=      E2*E3*AUX1+E1*AUX2
        REE33=COSAE+E3*E3*AUX1
      ENDIF
      CALL TCGTRA(REE,IME,REET,IMET)
C
C
C     Rewriting the beginning of the file with the model:
C     Lines to the end of surfaces:
      CALL RSEP3T('MODOUT',FILE2,'tco-mod.out' )
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      OPEN(LU2,FILE=FILE2)
      CALL NEWMOD(LU1,LU2,NSRF,NCB)
      CALL NEWVAL(LU1,LU2,1,NSRF,1,TSRF)
C     Lines to the end of S-wave velocity:
      N=-4
      CALL NEWLIN(LU1,LU2,N)
      CALL NEWLIN(LU1,LU2,N)
      CALL NEWLIN(LU1,LU2,N)
      CALL NEWLIN(LU1,LU2,N)
      CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*) VS
C
C     Reading A55:
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,*) (RAM(I1),I1=1,1001)
      READ(LU1,*) (RAM(I1),I1=1*1001+1,1*1001+1001)
C     Reading A44:
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,*) (RAM(I1),I1=1,1001)
      READ(LU1,*) (RAM(I1),I1=4*1001+1,4*1001+1001)
C     Reading A45:
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,'(A)') LINE
      READ(LU1,*) (RAM(I1),I1=1,1001)
      READ(LU1,*) (RAM(I1),I1=2*1001+1,2*1001+1001)
C
C     Rotating the model
C     (calculating new values of A55, A45, A35, A44, A34, A33):
      DO 10, I2=1,1001
        IMA11=0.
        IMA21=0.
        IMA31=0.
        IMA12=0.
        IMA22=0.
        IMA32=0.
        IMA13=0.
        IMA23=0.
        IMA33=0.
        REA11=RAM(1*1001+I2)
        REA21=RAM(2*1001+I2)
        REA31=0.
        REA12=RAM(2*1001+I2)
        REA22=RAM(4*1001+I2)
        REA32=0.
        REA13=0.
        REA23=0.
        REA33=(1.73205*VS)**2
        CALL TCGMAT(REE,IME,REA,IMA,REB,IMB)
        CALL TCGMAT(REB,IMB,REET,IMET,REA,IMA)
        RAM(1*1001+I2)=REA11
        RAM(2*1001+I2)=REA21
        RAM(3*1001+I2)=REA31
        RAM(4*1001+I2)=REA22
        RAM(5*1001+I2)=REA32
        RAM(6*1001+I2)=REA33
  10  CONTINUE
C
C     Writing:
      WRITE(LU2,*) ' '
      WRITE(LU2,*) '''A55'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''A45'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(2*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''A35'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(3*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''A44'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(4*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''A34'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(5*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''A33'' 1'
      WRITE(LU2,*) ' 3 0 0 0 /'
      WRITE(LU2,*) '1001'
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(1))
      CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,1001,RAM(6*1001+1))
      WRITE(LU2,*) ' '
C
      WRITE(LU2,*) '''END OF THE DATA FOR THE MODEL'' /'
C
      CLOSE(LU1)
      CLOSE(LU2)
      WRITE(*,'(A)') '+TCROT: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE NEWMOD(LU1,LU2,NSRF,NCB)
      INTEGER LU1,LU2,NSRF,NCB
C
C Subroutines and external functions required:
      EXTERNAL NEWLIN
C
C-----------------------------------------------------------------------
C
      CHARACTER*1 TEXTM
      INTEGER I,J,K,N,NEXPV,NEXPQ,NSB
      REAL AUX
C
C.......................................................................
C
C     Model description:
      N=0
   11 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=11) TEXTM
C
C     Model indices:
      N=0
   12 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
        NEXPV=1
        NEXPQ=1
      READ(LU2,*,END=12) I,NEXPV,NEXPQ,I
C
C     Model boundaries:
      N=0
   13 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=13) (AUX,I=1,6)
C
C     Number of surfaces:
      N=0
   14 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=14) NSRF
C
C     Number of simple blocks:
      N=0
   20 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=20) NSB
C
C     Indices of surfaces bounding simple blocks:
      DO 22 J=1,NSB
        N=0
   21   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
        READ(LU2,*,END=21) (K,I=1,99)
   22 CONTINUE
C
C     Number of complex blocks:
      N=0
   30 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=30) NCB
C
C     Indices of simple blocks forming complex blocks:
      DO 32 J=1,NCB
        N=0
   31   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
        READ(LU2,*,END=31) (K,I=1,99)
   32 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE NEWVAL(LU1,LU2,ICLASS,NGROUP,NFUNCT,TFUNCT)
      INTEGER LU1,LU2,ICLASS,NGROUP,NFUNCT
      CHARACTER*(*) TFUNCT(NFUNCT)
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL ERROR,WARRAY,VAL2,NEWLIN,LENGTH
      INTEGER LENGTH
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
      CHARACTER*3 TEXT
      CHARACTER*120 LINE
      CHARACTER*40 FORMAT
      LOGICAL WHAT
      INTEGER NVAR,IVAR(3),NX(3),MX,IGROPU,IFUNCT,JFUNCT,IADR
      INTEGER I1,I2,I3,I,N
      REAL GROUP,POWERW,COOR(3),F(10,47),POWER(47),AUX
C
C.......................................................................
C
C     Flag if the physical meaning of the functions is included in the
C     input data:
      WHAT=.FALSE.
      DO 10 I=1,NFUNCT
        IF(TFUNCT(I).NE.'   ') WHAT=.TRUE.
   10 CONTINUE
C
C     Loop for groups of functions:
      N=0
   11 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
        GROUP=1.
      READ(LU2,*,END=11) TEXT,GROUP
      DO 90 IGROPU=1,NGROUP
C
C       Loop for functions of the current group:
        DO 80 IFUNCT=1,NFUNCT
C
C         Physical meaning of the function:
          IF(WHAT) THEN
            N=0
   21       CONTINUE
              CALL NEWLIN(LU1,LU2,N)
              GROUP=1.
            READ(LU2,*,END=21) TEXT,GROUP
            DO 22 I=1,NFUNCT
              IF(TFUNCT(I).EQ.TEXT) THEN
                JFUNCT=I
                GO TO 23
              END IF
   22       CONTINUE
            GO TO 89
   23       CONTINUE
          ELSE
            JFUNCT=IFUNCT
          END IF
C
C         Initial address of the function parameters:
          I2=IPAR(ICLASS-1)+IGROPU
          DO 25 I1=IPAR(I2-1)+1,IPAR(I2-1)+NFUNCT
            IADR=IPAR(I1-1)
            IF(IPAR(IADR+1).EQ.JFUNCT) THEN
              GO TO 26
            END IF
   25     CONTINUE
C           TCROT-02
            CALL ERROR('TCROT-02: Function not found')
C           Function specified in data MODIN has not been specified in
C           data MODEL.
   26     CONTINUE
C
C         Reading spline grid:
          DO 31 I=1,3
            IVAR(I)=0
            NX(I)=1
   31     CONTINUE
          N=0
   32     CONTINUE
            CALL NEWLIN(LU1,LU2,N)
            IVAR(1)=0
            IVAR(2)=0
            IVAR(3)=0
            POWERW=1.
          READ(LU2,*,END=32) (IVAR(I),I=1,3),AUX,POWERW
          NVAR=3
          I2=0
   41     CONTINUE
            I2=I2+1
            IF(IVAR(I2).LE.0) THEN
              NVAR=NVAR-1
              DO 42 I1=I2,NVAR
                IVAR(I1)=IVAR(I1+1)
   42         CONTINUE
              I2=I2-1
            END IF
          IF(I2.LT.NVAR) GO TO 41
          IF(NVAR.GT.0) THEN
            N=0
   44       CONTINUE
              CALL NEWLIN(LU1,LU2,N)
            READ(LU2,*,END=44) (NX(I),I=1,NVAR)
          END IF
          MX=MAX0(NX(1),NX(2),NX(3))
          RAM(     1)=0.
          RAM(  MX+1)=0.
          RAM(2*MX+1)=0.
          IF(4*MX.GT.MRAM) THEN
C           TCROT-03
            CALL ERROR('TCROT-03: Small array RAM')
          END IF
          DO 46 I2=1,NVAR
            IF(NX(I2).GT.0) THEN
              N=0
   45         CONTINUE
                CALL NEWLIN(LU1,LU2,N)
              READ(LU2,*,END=45)
     *                         (RAM(I1),I1=(I2-1)*MX+1,(I2-1)*MX+NX(I2))
            ELSE
              NX(I2)=1
            END IF
   46     CONTINUE
          READ(LU1,*) (AUX,I=1,NX(1)*NX(2)*NX(3))
C
C         Changing coordinate indices to 1,2,3:
          DO 53 I2=3,5
            IF(IPAR(IADR+I2).LE.0) THEN
              IPAR(IADR+I2)=0
            ELSE
              DO 51 I1=1,NVAR
                IF(IPAR(IADR+I2).EQ.IVAR(I1)) THEN
                  IPAR(IADR+I2)=I1
                  GO TO 52
                END IF
   51         CONTINUE
C               TCROT-04
                CALL ERROR('TCROT-04: Wrong independent variable')
C               Function in data MODEL depends on different variables
C               than the corresponding function in data MODIN.
   52         CONTINUE
            END IF
   53     CONTINUE
C
C         Calculating and writing grid values of the given function:
          DO 63 I3=1,NX(3)
            IF(NX(1).NE.1.AND.NX(2).NE.1.AND.NX(3).NE.1) THEN
C             Separating 2-D slices of 3-D grid by a blank line
              WRITE(LU2,*)
            END IF
            COOR(3)=RAM(2*MX+I3)
            DO 62 I2=1,NX(2)
              COOR(2)=RAM(MX+I2)
              DO 61 I1=1,NX(1)
                COOR(1)=RAM(I1)
                CALL VAL2(ICLASS,IGROPU,NFUNCT,COOR,F,POWER)
                AUX=GROUP*POWERW/POWER(JFUNCT)
                RAM(3*MX+I1)=F(1,JFUNCT)
                IF(WHAT) THEN
                  IF(AUX.NE.1.) THEN
                    IF(RAM(3*MX+I1).GE.0.) THEN
                      RAM(3*MX+I1)=RAM(3*MX+I1)**AUX
                    ELSE
                      FORMAT='(A,I2,A,I2,A,'
                      CALL FORM2(3,COOR,COOR,FORMAT(14:37))
C                     
C                     TCROT-05
                      WRITE(LINE,FORMAT)
     *                  'TCROT-05: Negative value. Block',IGROPU,
     *                  ', function',JFUNCT,
     *                  ', coordinates ',COOR(1),' ',COOR(2),' ',COOR(3)
                      CALL ERROR(LINE(1:LENGTH(LINE)))
C                     Function with negative values is interpolated
C                     while its non-unit power should be written.
C                     Such an operation is not permitted.
                    END IF
                  END IF
                END IF
   61         CONTINUE
              CALL WARRAY(LU2,' ',.FALSE.,0.,.FALSE.,0.,
     *                                                NX(1),RAM(3*MX+1))
   62       CONTINUE
   63     CONTINUE
   80   CONTINUE
C       End of loop for functions
C
        N=0
   81   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
          GROUP=1.
        READ(LU2,*,END=81) TEXT,GROUP
   89   CONTINUE
   90 CONTINUE
C     End of loop for groups of functions
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE NEWLIN(LU1,LU2,N)
      INTEGER LU1,LU2,N
C
C Subroutines and external functions required:
      EXTERNAL LENGTH
      INTEGER  LENGTH
C
C-----------------------------------------------------------------------
C
      CHARACTER*251 LINE
      INTEGER I
C
C.......................................................................
C
C     Returning from the position after the end of file:
      IF(N.GT.0) THEN
        BACKSPACE(LU2)
      END IF
C
C     Copying one more line:
      READ (LU1,'(A)') LINE
      WRITE(LU2,'(A)') LINE(1:LENGTH(LINE))
      N=N+1
C
C     Rewinding to the position before reading:
      DO 10 I=1,N
        BACKSPACE(LU2)
   10 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TCGMAT(A,B,C,D,E,F)
C
C----------------------------------------------------------------------
C Subroutine to compute the product of two 3x3 complex matrices
C (E+iF)=(A+iB)*(C+iD)
      REAL A(3,3),B(3,3),C(3,3),D(3,3),E(3,3),F(3,3)
C Input:
C     A,B,C,D ... Real and imaginary parts of the two matrices.
C Output:
C     E,F ... Real and imaginary parts of the resulting matrix.
C
C.......................................................................
      E(1,1)=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1)+
     *       A(1,3)*C(3,1)-B(1,3)*D(3,1)
      E(2,1)=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1)+
     *       A(2,3)*C(3,1)-B(2,3)*D(3,1)
      E(3,1)=A(3,1)*C(1,1)-B(3,1)*D(1,1)+A(3,2)*C(2,1)-B(3,2)*D(2,1)+
     *       A(3,3)*C(3,1)-B(3,3)*D(3,1)
      E(1,2)=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2)+
     *       A(1,3)*C(3,2)-B(1,3)*D(3,2)
      E(2,2)=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2)+
     *       A(2,3)*C(3,2)-B(2,3)*D(3,2)
      E(3,2)=A(3,1)*C(1,2)-B(3,1)*D(1,2)+A(3,2)*C(2,2)-B(3,2)*D(2,2)+
     *       A(3,3)*C(3,2)-B(3,3)*D(3,2)
      E(1,3)=A(1,1)*C(1,3)-B(1,1)*D(1,3)+A(1,2)*C(2,3)-B(1,2)*D(2,3)+
     *       A(1,3)*C(3,3)-B(1,3)*D(3,3)
      E(2,3)=A(2,1)*C(1,3)-B(2,1)*D(1,3)+A(2,2)*C(2,3)-B(2,2)*D(2,3)+
     *       A(2,3)*C(3,3)-B(2,3)*D(3,3)
      E(3,3)=A(3,1)*C(1,3)-B(3,1)*D(1,3)+A(3,2)*C(2,3)-B(3,2)*D(2,3)+
     *       A(3,3)*C(3,3)-B(3,3)*D(3,3)
C
      F(1,1)=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1)+
     *       A(1,3)*D(3,1)+B(1,3)*C(3,1)
      F(2,1)=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1)+
     *       A(2,3)*D(3,1)+B(2,3)*C(3,1)
      F(3,1)=A(3,1)*D(1,1)+B(3,1)*C(1,1)+A(3,2)*D(2,1)+B(3,2)*C(2,1)+
     *       A(3,3)*D(3,1)+B(3,3)*C(3,1)
      F(1,2)=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2)+
     *       A(1,3)*D(3,2)+B(1,3)*C(3,2)
      F(2,2)=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2)+
     *       A(2,3)*D(3,2)+B(2,3)*C(3,2)
      F(3,2)=A(3,1)*D(1,2)+B(3,1)*C(1,2)+A(3,2)*D(2,2)+B(3,2)*C(2,2)+
     *       A(3,3)*D(3,2)+B(3,3)*C(3,2)
      F(1,3)=A(1,1)*D(1,3)+B(1,1)*C(1,3)+A(1,2)*D(2,3)+B(1,2)*C(2,3)+
     *       A(1,3)*D(3,3)+B(1,3)*C(3,3)
      F(2,3)=A(2,1)*D(1,3)+B(2,1)*C(1,3)+A(2,2)*D(2,3)+B(2,2)*C(2,3)+
     *       A(2,3)*D(3,3)+B(2,3)*C(3,3)
      F(3,3)=A(3,1)*D(1,3)+B(3,1)*C(1,3)+A(3,2)*D(2,3)+B(3,2)*C(2,3)+
     *       A(3,3)*D(3,3)+B(3,3)*C(3,3)
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TCGTRA(A,B,C,D)
C
C----------------------------------------------------------------------
C Subroutine to compute the transpose matrix to the 3x3 complex matrix
C (C+iD)=(A+iB)^T
      REAL A(3,3),B(3,3),C(3,3),D(3,3)
C Input:
C     A,B ... Real and imaginary parts of the input matrix.
C Output:
C     C,D ... Real and imaginary parts of the transposed matrix.
C
C.......................................................................
      C(1,1)=A(1,1)
      C(2,1)=A(1,2)
      C(3,1)=A(1,3)
      C(1,2)=A(2,1)
      C(2,2)=A(2,2)
      C(3,2)=A(2,3)
      C(1,3)=A(3,1)
      C(2,3)=A(3,2)
      C(3,3)=A(3,3)
C
      D(1,1)=B(1,1)
      D(2,1)=B(1,2)
      D(3,1)=B(1,3)
      D(1,2)=B(2,1)
      D(2,2)=B(2,2)
      D(3,2)=B(2,3)
      D(1,3)=B(3,1)
      D(2,3)=B(3,2)
      D(3,3)=B(3,3)
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
C
C=======================================================================
C