C 
C Program GRDWRL to convert gridded data into the GOCAD representation
C
C Version: 5.60
C Date: 2002, May 17
C
C Coded by: Ludek Klimes & Vaclav Bucha
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mails: klimes@seis.karlov.mff.cuni.cz
C              bucha@seis.karlov.mff.cuni.cz
C
C Reference:
C     
C     GOCAD
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 Input colour-map file:
C     COLORS='string'... Name of the file containing the data describing
C             the colour map.
C             Description of file COLORS
C             Default: COLORS='hsv.dat'
C Input/output file:
C     WRL='string'... Name of the file to be supplemented with surfaces
C             or to be copied to the beginning of the output file.
C             If the filename is blank, output file starts from a
C             scratch (mostly not reasonable).
C             The default name of the output file is equal to WRL.
C             It is recommended to specify WRL rather than to use
C             the default name.
C             Default: WRL='out.wrl'
C     WRLOUT='string'... Name of the output file if different from WRL.
C             Default: WRLOUT=WRL
C Data specifying the form of the output file:
C     VRML='string'... Virtual reality scene description language.
C             VRML='VRML1': VRML (Virtual Reality Modeling Language)
C                           version 1.0.  This value is accepted but
C                           no output is generated.
C             VRML='VRML2': VRML97 according to ISO/IEC 14772 standard.
C                           This value is accepted but no output is
C                           generated.
C             VRML='GOCAD': GOCAD description of Voxet (data grid) is
C                           generated.
C             Default: VRML='VRML2'.
C     NAME='string'... String containing the GOCAD name of the gridded
C             values.  Be sure to select different names for all objects
C             within the GOCAD file.
C             The same name is used for the corresponding colour scale.
C             Used only if VRML='GOCAD'.  Obligatory parameter, must be
C             specified and cannot be blank if VRML='GOCAD'.
C     PROPERTY='string'... String containing the GOCAD property name of
C             the gridded quantity.
C             Used only if VRML='GOCAD'.
C             Default: PROPERTY=NAME
C Data specifying the name of the file with gridded values:
C     GRD='string'... String with the name of the input ascii file
C             containing the gridded values.  The file is used to
C             determine the minimum and maximum grid values for colour
C             mapping.
C             Undefined grid values are allowed but dangerous.
C             No default, GRD must be specified and cannot be blank.
C     IN='string'... String with the name of the input binary file
C             containing the gridded values.  The file should contain
C             just the 4 byte big-endian (workstation-like) IEEE reals,
C             even on a little-endian computer.  The length of the file
C             is thus exactly 4*N1*N2*N3 bytes.  The gridded data should
C             be converted into this form by programs
C             ascbin.for and
C             swap.for.
C             Note that the file given by parameter IN is not read nor
C             checked by this program, only its name is written to the
C             output file.  The user is thus responsible for generating
C             the file with gridded values in the correct form.
C             No default, IN must be specified and cannot be blank.
C Data specifying the parameters of the input grid:
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     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 Data specifying the colour scale:
C     VADD=real, VMUL=real, VPER=real, VREF=real, CREF=real, CREF1=real,
C     CREF2=real, CREF3=real, etc... Refer to file
C             colors.for.
C     TRANSP=real... Transparency of the displayed sections through the
C             gridded data.  Values from 0 to 1.
C             Initial position of the sections corresponds to the six
C             sides of the voxet.
C             Default: TRANSP=0.
C     NODATA=real... Value used as the 'property no data value' in
C             GOCAD.  This parameter should not be specified if there
C             are undefined grid values.  Initial transparency of the
C             'no data' values is set to 0.80.
C             Default: NODATA=undefined value used in
C                      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     External functions and subroutines:
      EXTERNAL LENGTH,UARRAY,RSEP1,RSEP3T,RSEP3I,RSEP3R,ERROR,FORM2
      INTEGER  LENGTH
      REAL     UARRAY
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FGRD,FBIN,FIN,FOUT
      INTEGER LU1,LU2
      PARAMETER (LU1=1,LU2=2)
C
C     Other variables:
      CHARACTER*27 FORMAT
      CHARACTER*5  VRML
      CHARACTER*255 NAME,TEXT
      INTEGER N1,N2,N3,I0,I
      REAL O1,O2,O3,D1,D2,D3,TRANSP
      REAL OUT(3),OUTMIN(1),OUTMAX(1),R,G,B,AUX,AUXA(1)
C
C     Undefined grid value:
      REAL UNDEF
      UNDEF=UARRAY()
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+GRDWRL: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       GRDWRL-01
        CALL ERROR('GRDWRL-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
      CALL RSEP1(LU1,FSEP)
      WRITE(*,'(A)') '+GRDWRL: Working...            '
C
C     Reading input and output filenames:
      CALL RSEP3T('GRD',FGRD ,' ')
      IF(FGRD.EQ.' ') THEN
C       GRDWRL-02
        CALL ERROR('GRDWRL-02: No ascii grid file specified')
C       Ascii file with gridded data must be specified.
C       There is no default filename.
      END IF
      CALL RSEP3T('IN',FBIN ,' ')
      IF(FBIN.EQ.' ') THEN
C       GRDWRL-03
        CALL ERROR('GRDWRL-03: No binary grid file specified')
C       Binary file with gridded data must be specified.
C       There is no default filename.
      END IF
      CALL RSEP3T('WRL'   ,FIN  ,'out.wrl')
      CALL RSEP3T('WRLOUT',FOUT ,FIN      )
      CALL RSEP3T('VRML'  ,VRML ,'VRML2'  )
      CALL LOWER(VRML)
C
C     Reading grid dimensions:
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
C
C     Reading input parameters for surface appearance:
      CALL RSEP3R('TRANSP',TRANSP,0.00)
C
C     Opening the output file and writing its beginning:
      CALL WRL1(LU1,LU2,FIN,FOUT,VRML,1)
C
C     Writing the prolog for the voxet (part 1):
      IF (VRML.EQ.'vrml1') THEN
C       GRDWRL-51
        CALL WARN('GRDWRL-51: Nothing to do')
C       If VRML='VRML1', no output file with gridded values is created.
C       This program generates output only if VRML='GOCAD'.
      ELSE IF (VRML.EQ.'vrml2') THEN
C       GRDWRL-52
        CALL WARN('GRDWRL-52: Nothing to do')
C       If VRML='VRML2', no output file with gridded values is created.
C       This program generates output only if VRML='GOCAD'.
      ELSE IF (VRML.EQ.'gocad') THEN
        CALL RSEP3T('NAME',NAME,' ')
C       Subroutine WRL has already checked that NAME is not blank.
        WRITE(LU2,'(A)')
     *   'GOCAD Voxet 1.0'
        WRITE(LU2,'(2A)')
     *   'HDR name:',NAME(1:LENGTH(NAME))
        CALL RSEP3T('PROPERTY',TEXT,NAME)
        I=LENGTH(TEXT)
        WRITE(LU2,'(2A)')
     *   'HDR *property:',TEXT(1:I)
     *  ,'HDR *grid3d1:',TEXT(1:I)
        WRITE(LU2,'(2A)')
     *   'HDR ','*cage:false'
     *  ,'HDR ','*shaded:true'
     *  ,'HDR ','*precise:true'
     *  ,'HDR ','*smoothed:true'
        FORMAT='(3A,I0,A,I0,A,I0)'
        FORMAT( 6: 6)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N1)-0.5)))
        FORMAT(11:11)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N2)-0.5)))
        FORMAT(16:16)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N3)-0.5)))
        WRITE(LU2,FORMAT)
     *   'HDR *',TEXT(1:I),'*sections: 6  1 1 0  1 1 ',N1-1
     *                                ,'  2 1 0  2 1 ',N2-1
     *                                ,'  3 1 0  3 1 ',N3-1
        FORMAT='(A,'
        OUT(1)=O1
        OUT(2)=O2
        OUT(3)=O3
        CALL FORM2(3,OUT,OUT,FORMAT(4:27))
        WRITE(LU2,FORMAT)
     *   'AXIS_O ',O1,' ',O2,' ',O3
        OUT(1)=D1
        OUT(2)=D2
        OUT(3)=D3
        CALL FORM2(3,OUT,OUT,FORMAT(4:27))
        FORMAT(25:27)=')  '
        WRITE(LU2,FORMAT)
     *   'AXIS_U ',D1,' ',0.,' ',0.
     *  ,'AXIS_V ',0.,' ',D2,' ',0.
     *  ,'AXIS_W ',0.,' ',0.,' ',D3
        FORMAT='(A,I0,A,I0,A,I0)'
        FORMAT( 5: 5)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N1)+0.5)))
        FORMAT(10:10)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N2)+0.5)))
        FORMAT(15:15)=CHAR(ICHAR('1')+INT(ALOG10(FLOAT(N3)+0.5)))
        WRITE(LU2,FORMAT)
     *   'AXIS_MIN ',0   ,' ',0   ,' ',0
     *  ,'AXIS_MAX ',N1-1,' ',N2-1,' ',N3-1
     *  ,'AXIS_N   ',N1  ,' ',N2,  ' ',N3
        WRITE(LU2,'(3A)')
     *   'PROPERTY 1 "',TEXT(1:I),'"'
     *  ,'PROPERTY_CLASS 1 "',TEXT(1:I),'"'
     *  ,'PROPERTY_CLASS_HEADER ',TEXT(1:I),' {'
C       The output file now waits for the colour scale.
      ELSE
C       GRDWRL-04
        CALL ERROR('GRDWRL-04: Invalid string in VRML')
C       Valid string specifying the form of the output file is:
C       VRML='GOCAD', 'VRML1' or 'VRML2'.  The output is created only
C       if VRML='GOCAD'.  Default value is 'VRML2'.
      END IF
C
C     Determining the minimum and maximum grid values:
      IF(N1*N2*N3.GT.MRAM) THEN
C       GRDWRL-05
        CALL ERROR('GRDWRL-05: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM) to allocate the grid values.
C       If possible, increase dimension MRAM in include file
C       ram.inc.
      END IF
      CALL RARRAY(LU1,FGRD,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,RAM)
      I0=0
      DO 10 I=1,N1*N2*N3
        IF(RAM(I).NE.UNDEF) THEN
          IF(I0.EQ.0) THEN
            OUTMIN(1)=RAM(I)
            OUTMAX(1)=RAM(I)
            I0=1
          ELSE
            OUTMIN(1)=AMIN1(OUTMIN(1),RAM(I))
            OUTMAX(1)=AMAX1(OUTMAX(1),RAM(I))
          END IF
        END IF
   10 CONTINUE
      IF(I0.EQ.0) THEN
C       GRDWRL-06
        CALL ERROR('GRDWRL-06: All grid values undefined')
C       At least one grid value should be defined.
      END IF
C
C     Determining the colour map:
      IF(VRML.EQ.'gocad') THEN
        CALL COLOR1(LU1,MRAM,IRAM,RAM,1,OUTMIN,OUTMAX)
        WRITE(LU2,'(2A)')
     *   ' *colormap:',TEXT(1:LENGTH(TEXT))
        WRITE(LU2,'(A)')
     *   ' *colormap*nodata:true'
     *  ,' *colormap*ndtransparency: 0.80'
        FORMAT='(A,'
        CALL FORM2(1,OUTMIN,OUTMAX,FORMAT(4:11))
        FORMAT(9:11)=')  '
        IF(OUTMAX(1).GT.OUTMIN(1)) THEN
          WRITE(LU2,FORMAT)
     *     ' *low_clip: ',OUTMIN(1)
     *    ,' *high_clip:',OUTMAX(1)
        ELSE
          WRITE(LU2,FORMAT)
     *     ' *low_clip: ',OUTMIN(1)
     *    ,' *high_clip:',OUTMIN(1)+1.
        END IF
        WRITE(LU2,'(4A)')
     *   ' *colormap*',TEXT(1:LENGTH(TEXT)),'*colors: ',CHAR(92)
        AUX=(OUTMAX(1)-OUTMIN(1))/255.
        DO 31 I=0,255
          AUXA(1)=OUTMIN(1)+FLOAT(I)*AUX
          CALL COLOR2(MRAM,IRAM,RAM,1,AUXA,R,G,B)
          IF (I.LT.255) THEN
            WRITE(LU2,'(I5,3(1X,F4.2),2A)')
     *       I,R,G,B,' ',CHAR(92)
          ELSE
            WRITE(LU2,'(I5,3(1X,F4.2),2A)')
     *       I,R,G,B
          END IF
   31   CONTINUE
        IF(TRANSP.GT.0.) THEN
          WRITE(LU2,'(2A)')
     *     ' *colormap*alphas: ',CHAR(92)
          DO 32 I=0,255
            IF (I.LT.255) THEN
              WRITE(LU2,'(I5,1X,F4.2,2A)')
     *         I,TRANSP,' ',CHAR(92)
            ELSE
              WRITE(LU2,'(I5,1X,F4.2,2A)')
     *         I,TRANSP
            END IF
   32     CONTINUE
        END IF
      END IF
C
C     Writing the prolog for the voxet (part 2):
      IF (VRML.EQ.'gocad') THEN
        WRITE(LU2,'(A)')
     *   '}'
     *  ,'PROP_FORMAT 1 RAW'
     *  ,'PROP_ETYPE 1 IEEE'
     *  ,'PROP_ESIZE 1 4'
        CALL RSEP3R('NODATA',AUX,UNDEF)
        WRITE(LU2,'(A,G15.9)')
     *   'PROP_NO_DATA_VALUE 1 ',AUX
        WRITE(LU2,'(2A)')
     *   'PROP_FILE 1 ',FBIN
      END IF
C
C     Writing the trailor for the GOCAD object:
      IF (VRML.EQ.'gocad') THEN
        WRITE(LU2,'(A)') 'END'
      END IF
      CLOSE(LU2)
      WRITE(*,'(A)') '+GRDWRL: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'colors.for'
C     colors.for
      INCLUDE 'wrl.for'
C     wrl.for
C
C=======================================================================
C