C
C Program GRDISO for identification of points at isosurfaces
C of 3-D gridded values.  If one of two neigbouring gridpoints displays
C the value greater than VALUE, and if the second one of the two
C gridpoints displays the value less or eaqual to VALUE, the coordinates
C of the point which lies between the two gridpoints and corresponds
C to the value VALUE are computed and stored in the output file.
C
C Version: 6.70
C Date: 2012, November 7
C
C Coded by Petr Bulant
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.......................................................................
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 input file:
C     GRD='string'... Name of the input ASCII file with the grid values.
C             Default: GRD='grd.out'
C     For general description of the files with gridded data refer
C     to file forms.htm.
C Parameters defining the basic regular rectangular grid:
C     N1=positive integer... Number of gridpoints of the basic grid
C             along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints of the basic grid
C             along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints of the basic grid
C             along the X3 axis.
C             Default: N3=1
C     O1=real... X1 coordinate of the origin of the grid.
C             Default: O1=0.
C     O2=real... X2 coordinate of the origin of the grid.
C             Default: O2=0.
C     O3=real... X3 coordinate of the origin of the grid.
C             Default: O3=0.
C     D1=real... Grid spacing along the X1 axis.
C             Default: D1=1.
C     D2=real... Grid spacing along the X2 axis.
C             Default: D2=1.
C     D3=real... Grid spacing along the X3 axis.
C             Default: D3=1.
C Data specifying dimensions of the grid, on which the calculation
C of points at isosurfaces will be performed:
C     N1NEW=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1NEW=N1
C     N2NEW=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2NEW=N2
C     N3NEW=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3NEW=N3
C     NO1=positive integer... Index of the first gridpoint along
C             the X1 axis.
C             Default: NO1=1
C     NO2=positive integer... Index of the first gridpoint along
C             the X2 axis.
C             Default: NO2=1
C     NO3=positive integer... Index of the first gridpoint along
C             the X3 axis.
C             Default: NO3=1
C     ND1=positive integer... Multiplication factor of the grid interval
C             along the X1 axis.
C             Default: ND1=1
C     ND2=positive integer... Multiplication factor of the grid interval
C             along the X2 axis.
C             Default: ND2=1
C     ND3=positive integer... Multiplication factor of the grid interval
C             along the X3 axis.
C             Default: ND3=1
C     The grid for calculation should be always a subgrid of the
C     original grid.  Two gridpoints located at the different sides of
C     the isosurface are searched for on the grid for calculation.
C     The coordinates of the point at isosurface are then calculated on
C     the original grid.
C Value corresponding to the points to be calculated:
C     VALUE=real... Value corresponding to the isosurface being
C             searched for.
C             Default: VALUE=0.
C Mode of the calculation:
C     MODE=integer
C             MODE=0... The points corresponding to VALUE are
C               searched for at all gridlegs.
C             MODE=1,2,3... The points are searched for only
C               at the gridlegs parallel with X1, X2 or X3 axis.
C               All the points are written to the file PTS.
C             MODE=-1,-2,-3... The points are searched for only
C               at the gridlegs parallel with X1, X2 or X3 axis.
C               The points with values increasing with
C               the corresponding gridlegs are written to the file PTS,
C               other points to the file PTS1.
C             Default: MODE=0
C Names of the output files:
C     PTS='string'... Name of the output file with the coordinates
C             of the points. The points are the intersections of
C             the gridlegs with the isosurface of the value VALUE.
C             The file is not generated if PTS=' '.
C             Description of file PTS.
C             Default: PTS='pts.out'
C     PTS1='string'... Name of the second output file with the points
C             in case of negative MODE.
C             The file is not generated if PTS1=' '.
C             Description of file PTS1.
C             Default: PTS1='pts1.out'
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 Output file PTS or PTS1 with the points at isosurface:
C (1) /
C (2) For each point data (2.1):
C (2.1) 'NNNNNN',X1,X2,X3,/
C     'NNNNNN'...  Name of the point - six-digit integer index of the
C             point.
C      X1,X2,X3... Coordinates of the point.
C (3) /
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL GIGLEG,GICP,
     *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM1
C     GIGLEG,GICP... This file.
C     ERROR... File error.for.
C     RSEP1,RSEP3T,RSEP3R,RSEP3I... File sep.for.
C     FORM1... File forms.for.
C
C Common block /GIC/:
      INCLUDE 'grdiso.inc'
C grdiso.inc
C
C.......................................................................
C
      EXTERNAL UARRAY
      REAL UARRAY
C
C     Auxiliary storage locations:
      CHARACTER*80 FSEP,FGRD,FPTS,FPTS1
      INTEGER LU,LU1,LENG,MODE
      PARAMETER (LU=1,LU1=2,LENG=6)
      REAL UNDEF,VALUE
      INTEGER I1,I2,I3,I4,J2,JCOOR,IP1,IP2,IP3,NPTS,
     *  NLEG2,NLEG3
      REAL XX(3),YY(3),ZZ(3),OUTMIN,OUTMAX,VAL1,VAL2,VAL3
      CHARACTER*20 FORMAT,TEXTP
C
      UNDEF=UARRAY()
C
C.......................................................................
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDISO: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.EQ.' ') THEN
C       GRDISO-01
        CALL ERROR('GRDISO-01: SEP file not given')
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.
      ENDIF
C
C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',DD(1),1.)
      CALL RSEP3R('D2',DD(2),1.)
      CALL RSEP3R('D3',DD(3),1.)
      CALL RSEP3I('N1NEW',N1N,N1)
      CALL RSEP3I('N2NEW',N2N,N2)
      CALL RSEP3I('N3NEW',N3N,N3)
      IF (N1.EQ.N1N.AND.N2.EQ.N2N.AND.N3.EQ.N3N) THEN
        NO1=1
        NO2=1
        NO3=1
        ND1=1
        ND2=1
        ND3=1
      ELSE
        CALL RSEP3I('NO1',NO1,1)
        CALL RSEP3I('NO2',NO2,1)
        CALL RSEP3I('NO3',NO3,1)
        CALL RSEP3I('ND1',ND1,1)
        CALL RSEP3I('ND2',ND2,1)
        CALL RSEP3I('ND3',ND3,1)
        IF ((NO1+ND1*(N1N-1).GT.N1).OR.
     *      (NO2+ND2*(N2N-1).GT.N2).OR.
     *      (NO3+ND3*(N3N-1).GT.N3)) THEN
C         GRDISO-02
          CALL ERROR('GRDISO-02: Wrong grid for calculation')
C         The grid for calculation must be a subgrid
C         of the original grid.
        ENDIF
      ENDIF
C
C     Mode of the calculation and the isovalue:
      CALL RSEP3I('MODE',MODE,0)
      CALL RSEP3R('VALUE',VALUE,0.)
C
C     Recalling the input and output filenames:
      CALL RSEP3T('GRD',FGRD,'grd.out')
      CALL RSEP3T('PTS',FPTS,'pts.out')
      IF (MODE.LT.0) CALL RSEP3T('PTS1',FPTS1,'pts1.out')
C
C
C     Reading input grid values:
      IF (N1*N2*N3.GT.MRAM) THEN
C       GRDISO-03
        CALL ERROR('GRDISO-03: Small array RAM.')
C       Try to enlarge the dimension MRAM in the file
C       ram.inc.
      ENDIF
      CALL RARRAY(LU,FGRD,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,RAM)
      NPTS=0
C
C     Opening output files, initializing:
      WRITE(*,'(A)') '+GRDISO: Working ... '
      IF (FPTS.NE.' ') THEN
        OPEN(LU,FILE=FPTS,FORM='FORMATTED')
        WRITE(LU,'(A)') '/'
      ENDIF
      IF ((MODE.LT.0).AND.(FPTS1.NE.' ')) THEN
        OPEN(LU1,FILE=FPTS1,FORM='FORMATTED')
        WRITE(LU1,'(A)') '/'
      ENDIF
C
      N1N2=   N1   * N2
C
      NLEG1=(N1N-1)*N2N*N3N
      NLEG2=N1N*(N2N-1)*N3N
      NLEG3=N1N*N2N*(N3N-1)
      NLEG12=NLEG1+NLEG2
      NLEG=NLEG12+NLEG3
C
      NN11=    N1N-1
      NN1N2=   N1N   * N2N
      NN11N2= (N1N-1)* N2N
      NN1N21=  N1N   *(N2N-1)
C     DIP - number of points in RAM between neighboring points
C           in calculation grid
      DIP(1)=1
      DIP(2)=N1
      DIP(3)=N1*N2
C     I1 - index of the gridleg being processed
      I2=1
      I3=NLEG
      I4=IABS(MODE)
      IF     (I4.EQ.1) THEN
        I3=NLEG1
      ELSEIF (I4.EQ.2) THEN
        I2=NLEG1
        I3=NLEG12
      ELSEIF (I4.EQ.3) THEN
        I2=NLEG12
      ENDIF
C     Loop along all gridlegs, recording points at isosurfaces:
      DO 29, I1=I2,I3
        CALL GIGLEG(I1,IP1,IP2)
        VAL1=RAM(IP1)
        VAL2=RAM(IP2)
        IF (((VAL1.LE.VALUE).AND.(VAL2.GT.VALUE)).OR.
     *      ((VAL1.GT.VALUE).AND.(VAL2.LE.VALUE))) THEN
C         The points of the gridleg are at different sides
C         of the isosurface - searching for two neighbouring
C         gridpoints of the original grid, which are also
C         at the different sides of the isosurface:
C         JCOOR - index of the coordinate axis
          JCOOR=3
          IF (I1.LE.NLEG12) JCOOR=2
          IF (I1.LE.NLEG1)  JCOOR=1
C         J2 - number of gridsteps of the original grid
C              between IP1 and IP2
  10      CONTINUE
            CALL GICP(IP1,XX)
            CALL GICP(IP2,YY)
            J2=IABS(NINT((XX(JCOOR)-YY(JCOOR))/DD(JCOOR)))
            J2=J2/2
            IF (J2.NE.0) THEN
              IP3=IP1+J2*DIP(JCOOR)
              VAL3=RAM(IP3)
              IF (((VAL1.LE.VALUE).AND.(VAL3.LE.VALUE)).OR.
     *            ((VAL1.GT.VALUE).AND.(VAL3.GT.VALUE))) THEN
                IP1=IP3
                VAL1=VAL3
              ELSE
                IP2=IP3
                VAL2=VAL3
              ENDIF
              GOTO 10
            ENDIF
C         End of the search for neighbouring gridpoints.
C         Computing the coordinates ZZ of the point at the isosurface:
          ZZ(1)=XX(1)
          ZZ(2)=XX(2)
          ZZ(3)=XX(3)
          IF ((VAL2-VAL1).NE.0.) ZZ(JCOOR)=
     *      XX(JCOOR)+(VALUE-VAL1)*(YY(JCOOR)-XX(JCOOR))/(VAL2-VAL1)
C         Writing the point:
C         Name of the point:
          NPTS=NPTS+1
          DO 204, I2=0,LENG-1
            TEXTP(LENG-I2:LENG-I2)=
     *      CHAR(ICHAR('0')+MOD(NPTS,10**(I2+1))/10**I2)
 204      CONTINUE
C         Setting the output format:
          FORMAT='(3A,03(F00.0,1X),A)'
          OUTMIN=0.
          OUTMAX=0.
          DO 214, I2=1,3
            IF (OUTMIN.GT.ZZ(I2)) OUTMIN=ZZ(I2)
            IF (OUTMAX.LT.ZZ(I2)) OUTMAX=ZZ(I2)
 214      CONTINUE
          CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15))
          FORMAT(14:17)= '1X),'
C         Writing:
          IF ((MODE.LT.0).AND.(VAL1.GT.VAL2).AND.(FPTS1.NE.' ')) THEN
            WRITE(LU1,FORMAT) '''',TEXTP(1:LENG),''' ',ZZ,'/'
          ELSEIF (FPTS.NE.' ') THEN
            WRITE(LU,FORMAT) '''',TEXTP(1:LENG),''' ',ZZ,'/'
          ENDIF
        ENDIF
  29  CONTINUE
      IF (FPTS.NE.' ') THEN
        WRITE(LU,'(A)') '/'
        CLOSE(LU)
      ENDIF
      IF ((MODE.LT.0).AND.(FPTS1.NE.' ')) THEN
        WRITE(LU1,'(A)') '/'
        CLOSE(LU1)
      ENDIF
C
      WRITE(*,'(A)')
     *'+GRDISO: Done.                                        '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE GIGLEG(ILEG,IPOIN1,IPOIN2)
C
C-----------------------------------------------------------------------
C
      INTEGER ILEG,IPOIN1,IPOIN2
C Input:
C     ILEG... Index of gridleg.
C Output:
C     IPOIN1,IPOIN2... Indices of gridpoints forming the gridleg.
C.......................................................................
C Common block /GIC/:
      INCLUDE 'grdiso.inc'
C grdiso.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,I2,I3,J1,J2,J3,J21,J31,JCOOR,ILEG1
C.......................................................................
C
      ILEG1=ILEG-1
      IF (ILEG.LE.NLEG1) THEN
        I31=ILEG1 / NN11N2
        I21=(ILEG1 - I31*NN11N2) / NN11
        I1=ILEG - I31*NN11N2 - I21*NN11
        JCOOR=1
      ELSEIF (ILEG.LE.NLEG12) THEN
        I31=(ILEG1-NLEG1) / NN1N21
        I21=((ILEG1-NLEG1) - I31*NN1N21) / N1N
        I1=(ILEG-NLEG1) - I31*NN1N21 - I21*N1N
        JCOOR=2
      ELSEIF (ILEG.LE.NLEG) THEN
        I31=(ILEG1-NLEG12) / NN1N2
        I21=((ILEG1-NLEG12) - I31*NN1N2) / N1N
        I1=(ILEG-NLEG12) - I31*NN1N2 - I21*N1N
        JCOOR=3
      ELSE
C       GRDISO-04
        CALL ERROR ('GRDISO-04: Wrong ILEG.')
C       This error should not appear. Contact the author.
      ENDIF
      I2=I21+1
      I3=I31+1
      J1=NO1+(I1-1)*ND1
      J2=NO2+(I2-1)*ND2
      J3=NO3+(I3-1)*ND3
      J21=J2-1
      J31=J3-1
      IPOIN1=J31*N1N2+J21*N1+J1
      IPOIN2=IPOIN1+DIP(JCOOR)
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE GICP(IPOIN,XX)
C
C-----------------------------------------------------------------------
C
      INTEGER IPOIN
      REAL XX(3)
C     IPOIN... Index of a point.
C     XX...   Coordinates of the point.
C.......................................................................
C Common block /GIC/:
      INCLUDE 'grdiso.inc'
C grdiso.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I11,IPOIN1
C.......................................................................
C
      IPOIN1=IPOIN-1
      I31= IPOIN1 / N1N2
      I21=(IPOIN1 - I31*N1N2) / N1
      I11= IPOIN  - I31*N1N2 - I21*N1 - 1
      XX(1)=O1+FLOAT(I11)*DD(1)
      XX(2)=O2+FLOAT(I21)*DD(2)
      XX(3)=O3+FLOAT(I31)*DD(3)
      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
C
C=======================================================================
C