C 
C Program GRDTE to compute the values of a real or complex function,
C described in terms of the Taylor expansions of its amplitude and
C phase, on a given grid
C
C Version: 5.40
C Date: 2000, May 9
C
C Coded by Karel Zacek
C     Department of Geophysics, Charles University Prague
C     E-mail: zacek@karel.troja.mff.cuni.cz
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 Names of the output files:
C     GRDTER='string'
C     GRDTEI='string'... Names of the output files with the calculated
C             grid values of the given complex-valued function.
C             File GRDTER will contain the real part and file GRDTEI
C             will contain the imaginary part of the function.
C             If the string is blank, the corresponding file is not
C             generated.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Default: GRDTER=' ', GRDTEI=' '
C Data specifying the dimensions of the grid for discretization:
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Default: O1=0.  O2=0.  O3=0.
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=0.
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=0.
C     D3=real... Grid interval along the X3 axis.
C             Default: D3=0.
C Data specifying the origins of individual discretized functions:
C     N10, N20, N30, O10, O20, O30, D10, D20 and D30... Specify the set
C             of origins X0j for individual shifted functions.  Their
C             meaning is analogous to parameters N1, N2, N3, O1, O2, O3,
C             D1, D2 and D3.
C Additional data specific to this program:
C     COEFFICIENT=real
C             where COEFFICIENT stands for any of the following 80 names
C               AR0, AI0, TR0, TI0,
C               AR1, AR2, AR3, AI1, AI2, AI3,
C               TR1, TR2, TR3, TI1, TI2, TI3,
C               AR11, AR12, AR22, AR13, AR23, AR33,
C               AI11, AI12, AI22, AI13, AI23, AI33,
C               TR11, TR12, TR22, TR13, TR23, TR33,
C               TI11, TI12, TI22, TI13, TI23, TI33,
C               AR111, AR112, AR122, AR222, AR113,
C               AR123, AR223, AR133, AR233, AR333,
C               AI111, AI112, AI122, AI222, AI113,
C               AI123, AI223, AI133, AI233, AI333,
C               TR111, TR112, TR122, TR222, TR113,
C               TR123, TR223, TR133, TR233, TR333,
C               TI111, TI112, TI122, TI222, TI113,
C               TI123, TI223, TI133, TI233, TI333.
C             The calculated complex-valued function has the form of
C                F(Xm) = A(Xm) exp(T(n)) ,
C             where the complex-valued amplidtude is
C                A(Xm) = A0 + Aj*Yj + Ajk*Yj*Yk/2 + Ajkl*Yj*Yk*Yl/6
C             and the complex-valued phase is
C                T(Xm) = T0 + Tj*Yj + Tjk*Yj*Yk/2 + Tjkl*Yj*Yk*Yl/6 ,
C             with
C                Yj = Xj - X0j .
C             Here
C               A0   = AR0   + i*AI0   ,
C               Aj   = ARj   + i*AIj   ,
C               Ajk  = ARjk  + i*AIjk  ,
C               Ajkl = ARjkl + i*AIjkl ,
C               T0   = TR0   + i*TI0   ,
C               Tj   = TRj   + i*TIj   ,
C               Tjk  = TRjk  + i*TIjk  ,
C               Tjkl = TRjkl + i*TIjkl .
C             The origins X0j of individual discretized functions are
C             specified by input parameters N10, N20, N30, O10, O20,
C             O30, D10, D20 and D30.  The individual functions are
C             written as individual snapshots.  It means that they may
C             be processed or displayed like N4=N10*N20*N30 spatial
C             grids.
C             Defaults: COEFFICIENT=0
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
      INTEGER N1,N2,N3,N10,N20,N30,LU1,LU2,PT,PTM
      REAL    D1,D2,D3,O1,O2,O3,D10,D20,D30,O10,O20,O30
      INTEGER IN1,IN2,IN3,IN10,IN20,IN30,I,J,K,L
      REAL    X(3),X0(3)
      REAL    FR,FI,AR,AI,TR,TI,AR0,AI0,TR0,TI0,DX,DXI,DXJ,DXK
      REAL    AR1(3),AR2(3,3),AR3(3,3,3),TR1(3),TR2(3,3),TR3(3,3,3)
      REAL    AI1(3),AI2(3,3),AI3(3,3,3),TI1(3),TI2(3,3),TI3(3,3,3)
      CHARACTER*80 FILE,OUTR,OUTI
      PARAMETER (LU1=1,LU2=2)
C
C-----------------------------------------------------------------------
C
C     Reading main input data:
      WRITE(*,'(A)')'+GRDTE: Enter input filename: '
      FILE=' '
      READ(*,*) FILE
      IF(FILE.EQ.' ') THEN
C       GRDTE-01
        CALL ERROR('GRDTE-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)')'+GRDTE: Reading, calculating... '
C
C     Reading grid dimensions:
      CALL RSEP1(LU,FILE)
      CALL RSEP3T('GRDTER',OUTR,' ')
      CALL RSEP3T('GRDTEI',OUTI,' ')
      IF (OUTR.EQ.' ' .AND. OUTI.EQ.' ') THEN
C       GRDFFT-03
        CALL ERROR('GRDFFT-03: No output files specified.')
      ENDIF
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('N10',N10,1)
      CALL RSEP3I('N20',N20,1)
      CALL RSEP3I('N30',N30,1)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3R('D10',D10,1.)
      CALL RSEP3R('D20',D20,1.)
      CALL RSEP3R('D30',D30,1.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('O10',O10,0.)
      CALL RSEP3R('O20',O20,0.)
      CALL RSEP3R('O30',O30,0.)
C
      N1N2N3=N1*N2*N3
      IF(2*N1N2N3.GT.MRAM) THEN
C       GRDTE-02
        CALL ERROR('GRDTE-02: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM) to allocate both N1*N2*N3 real and
C       N1*N2*N3 imaginary output grid values.  If possible, increase
C       dimension MRAM in include file ram.inc.
      END IF
      IF(OUTR.NE.' ') THEN
        OPEN (LU1,FILE=OUTR,FORM='FORMATTED')
      END IF
      IF(OUTI.NE.' ') THEN
        OPEN (LU2,FILE=OUTI,FORM='FORMATTED')
      END IF
      PTM=N10*N20*N30
      PT=0
C
C     Reading function parameters
      CALL RSEP3R('AR0',AR0,0.)
      CALL RSEP3R('AI0',AI0,0.)
      CALL RSEP3R('TR0',TR0,0.)
      CALL RSEP3R('TI0',TI0,0.)
      CALL RSEP3R('AR1',AR1(1),0.)
      CALL RSEP3R('AI1',AI1(1),0.)
      CALL RSEP3R('TR1',TR1(1),0.)
      CALL RSEP3R('TI1',TI1(1),0.)
      CALL RSEP3R('AR2',AR1(2),0.)
      CALL RSEP3R('AI2',AI1(2),0.)
      CALL RSEP3R('TR2',TR1(2),0.)
      CALL RSEP3R('TI2',TI1(2),0.)
      CALL RSEP3R('AR3',AR1(3),0.)
      CALL RSEP3R('AI3',AI1(3),0.)
      CALL RSEP3R('TR3',TR1(3),0.)
      CALL RSEP3R('TI3',TI1(3),0.)
      CALL RSEP3R('AR11',AR2(1,1),0.)
      CALL RSEP3R('AR12',AR2(1,2),0.)
      CALL RSEP3R('AR22',AR2(2,2),0.)
      CALL RSEP3R('AR13',AR2(1,3),0.)
      CALL RSEP3R('AR23',AR2(2,3),0.)
      CALL RSEP3R('AR33',AR2(3,3),0.)
      CALL RSEP3R('AI11',AI2(1,1),0.)
      CALL RSEP3R('AI12',AI2(1,2),0.)
      CALL RSEP3R('AI22',AI2(2,2),0.)
      CALL RSEP3R('AI13',AI2(1,3),0.)
      CALL RSEP3R('AI23',AI2(2,3),0.)
      CALL RSEP3R('AI33',AI2(3,3),0.)
      CALL RSEP3R('TR11',TR2(1,1),0.)
      CALL RSEP3R('TR12',TR2(1,2),0.)
      CALL RSEP3R('TR22',TR2(2,2),0.)
      CALL RSEP3R('TR13',TR2(1,3),0.)
      CALL RSEP3R('TR23',TR2(2,3),0.)
      CALL RSEP3R('TR33',TR2(3,3),0.)
      CALL RSEP3R('TI11',TI2(1,1),0.)
      CALL RSEP3R('TI12',TI2(1,2),0.)
      CALL RSEP3R('TI22',TI2(2,2),0.)
      CALL RSEP3R('TI13',TI2(1,3),0.)
      CALL RSEP3R('TI23',TI2(2,3),0.)
      CALL RSEP3R('TI33',TI2(3,3),0.)
      CALL RSEP3R('AR111',AR3(1,1,1),0.)
      CALL RSEP3R('AR112',AR3(1,1,2),0.)
      CALL RSEP3R('AR122',AR3(1,2,2),0.)
      CALL RSEP3R('AR222',AR3(2,2,2),0.)
      CALL RSEP3R('AR113',AR3(1,1,3),0.)
      CALL RSEP3R('AR123',AR3(1,2,3),0.)
      CALL RSEP3R('AR223',AR3(2,2,3),0.)
      CALL RSEP3R('AR133',AR3(1,3,3),0.)
      CALL RSEP3R('AR233',AR3(2,3,3),0.)
      CALL RSEP3R('AR333',AR3(3,3,3),0.)
      CALL RSEP3R('AI111',AI3(1,1,1),0.)
      CALL RSEP3R('AI112',AI3(1,1,2),0.)
      CALL RSEP3R('AI122',AI3(1,2,2),0.)
      CALL RSEP3R('AI222',AI3(2,2,2),0.)
      CALL RSEP3R('AI113',AI3(1,1,3),0.)
      CALL RSEP3R('AI123',AI3(1,2,3),0.)
      CALL RSEP3R('AI223',AI3(2,2,3),0.)
      CALL RSEP3R('AI133',AI3(1,3,3),0.)
      CALL RSEP3R('AI233',AI3(2,3,3),0.)
      CALL RSEP3R('AI333',AI3(3,3,3),0.)
      CALL RSEP3R('TR111',TR3(1,1,1),0.)
      CALL RSEP3R('TR112',TR3(1,1,2),0.)
      CALL RSEP3R('TR122',TR3(1,2,2),0.)
      CALL RSEP3R('TR222',TR3(2,2,2),0.)
      CALL RSEP3R('TR113',TR3(1,1,3),0.)
      CALL RSEP3R('TR123',TR3(1,2,3),0.)
      CALL RSEP3R('TR223',TR3(2,2,3),0.)
      CALL RSEP3R('TR133',TR3(1,3,3),0.)
      CALL RSEP3R('TR233',TR3(2,3,3),0.)
      CALL RSEP3R('TR333',TR3(3,3,3),0.)
      CALL RSEP3R('TI111',TI3(1,1,1),0.)
      CALL RSEP3R('TI112',TI3(1,1,2),0.)
      CALL RSEP3R('TI122',TI3(1,2,2),0.)
      CALL RSEP3R('TI222',TI3(2,2,2),0.)
      CALL RSEP3R('TI113',TI3(1,1,3),0.)
      CALL RSEP3R('TI123',TI3(1,2,3),0.)
      CALL RSEP3R('TI223',TI3(2,2,3),0.)
      CALL RSEP3R('TI133',TI3(1,3,3),0.)
      CALL RSEP3R('TI233',TI3(2,3,3),0.)
      CALL RSEP3R('TI333',TI3(3,3,3),0.)
C
C
C     Loop over Xi0:
      DO 120 IN30=1,N30
        X0(3)=O30+D30*(IN30-1)
        DO 110 IN20=1,N20
          X0(2)=O20+D20*(IN20-1)
          DO 100 IN10=1,N10
            X0(1)=O10+D10*(IN10-1)
            L=0
C
C           Loop over Xi:
            DO 90 IN3=1,N3
              X(3)=O3+D3*(IN3-1)
              DO 80 IN2=1,N2
                X(2)=O2+D2*(IN2-1)
                DO 70 IN1=1,N1
                  X(1)=O1+D1*(IN1-1)
                  AR=AR0
                  AI=AI0
                  TR=TR0
                  TI=TI0
C                 Calculating the functional value:
                  DO 10 I=1,3
                    DX=X(I)-X0(I)
                    AR=AR+AR1(I)*DX
                    AI=AI+AI1(I)*DX
                    TR=TR+TR1(I)*DX
                    TI=TI+TI1(I)*DX
   10             CONTINUE
                  DO 30 J=1,3
                    DXJ=(X(J)-X0(J))/2.
                    DO 20 I=1,J
                      DXI=X(I)-X0(I)
                      DX=DXI*DXJ
                      AR=AR+AR2(I,J)*DX
                      AI=AI+AI2(I,J)*DX
                      TR=TR+TR2(I,J)*DX
                      TI=TI+TI2(I,J)*DX
   20               CONTINUE
   30             CONTINUE
                  DO 60 K=1,3
                    DXK=(X(K)-X0(K))/6.
                    DO 50 J=1,K
                      DXJ=X(J)-X0(J)
                      DO 40 I=1,J
                        DXI=X(I)-X0(I)
                        DX=DXI*DXJ*DXK
                        AR=AR+AR3(I,J,K)*DX
                        AI=AI+AI3(I,J,K)*DX
                        TR=TR+TR3(I,J,K)*DX
                        TI=TI+TI3(I,J,K)*DX
   40                 CONTINUE
   50               CONTINUE
   60             CONTINUE
                  FR=AR*COS(TI)*EXP(TR)-AI*SIN(TI)*EXP(TR)
                  FI=AR*SIN(TI)*EXP(TR)+AI*COS(TI)*EXP(TR)
                  L=L+1
                  RAM(L)=FR
                  RAM(N1N2N3+L)=FI
C
C                 End of calculation
C
   70           CONTINUE
   80         CONTINUE
   90       CONTINUE
C
C     Writing output grid values
      IF(OUTR.NE.' ') THEN
        CALL WARRAY(LU1,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                             N1N2N3,RAM)
      END IF
      IF(OUTI.NE.' ') THEN
        CALL WARRAY(LU2,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                             N1N2N3,RAM(N1N2N3+1))
      END IF
C
C     End of loop over Xi
C
          PT=PT+1
          WRITE(*,'(''+GRDTE: '',I16,'' loops over Xi0 of'',I9)') PT,PTM
  100     CONTINUE
  110   CONTINUE
  120 CONTINUE
C
C     End of loop over Xi0
C
      IF(OUTR.NE.' ') THEN
        CLOSE(LU1)
      END IF
      IF(OUTI.NE.' ') THEN
        CLOSE(LU2)
      END IF
C
      WRITE(*,'(A)')'+GRDTE: Done.'
      STOP
      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
C
C=======================================================================
C