C                  INTERNATIONAL AVS CENTER
C        (This disclaimer must remain at the top of all files)
C
C WARRANTY DISCLAIMER
C
C This module and the files associated with it are distributed free of charge.
C It is placed in the public domain and permission is granted for anyone to use,
C duplicate, modify, and redistribute it unless otherwise noted.  Some modules
C may be copyrighted.  You agree to abide by the conditions also included in
C the AVS Licensing Agreement, version 1.0, located in the main module
C directory located at the International AVS Center ftp site and to include
C the AVS Licensing Agreement when you distribute any files downloaded from
C that site.
C
C The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module provide absolutely
C NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
C the quality and performance of this software is with the user.  IN NO EVENT
C WILL The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module BE LIABLE TO
C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.
C
C This AVS module and associated files are public domain software unless
C otherwise noted.  Permission is hereby granted to do whatever you like with
C it, subject to the conditions that may exist in copyrighted materials. Should
C you wish to make a contribution toward the improvement, modification, or
C general performance of this module, please send us your comments:  why you
C liked or disliked it, how you use it, and most important, how it helps your
C work. We will receive your comments at avs@ncsc.org.
C
C Please send AVS module bug reports to avs@ncsc.org.
C
C/* NOTE:  THIS MODULE AND SOURCE CODE IS FOR USE
C   WITH THE AVS SOFTWARE ENVIRONMENT ONLY */
C---------------------------------------------------------------------r
C Read EJ65
C    Read the EJ65 Binary Files containing 3D array data,
C    and generate an AVS Field output port.
C    This is a template read module in Fortran,
C    based on the standard example "test_fld2_f.f"
C    The EJ65 format has two files, the mesh and the flow field,
C    describing a single block mesh about wing-body
C    aircraft configuration.
C
C    This module performs the following operations:
C    1. Reads an Input file for geometry coordinates
C       The double prec. coords are converted to reals
C    2. Reads an Input file for flow domain variables
C       The flow field is cell-centered, and hence has
C       one less point in each dimension.
C    3. Sends the data to the AVS Output port,
C       in a 3D array.
C
C    It will do these things every time, even if you
C    have not changed the mesh file, for instance.
C
C----------------------------------------------------------------------
C Author: Ian J. Curington, Advanced Visual Systems, Inc.
C                           Staines, UK
C         On request of Dr. D. Wood, BTRL (generalized template)
C         Converted to BAe EJ65 Specific Format, for J. Davies
C
C   British Aerospace Aerodynamics, Filton, UK
C
C Revision History:
C  7 July 1992   ianc - Original
C  8 July 1992   ianc - Added Processing dials
C 16 April 1993  ianc - converted read_hologram.f to read_ej65
C 29 April 1993  ianc - added freeze button, port to SGI
C                       added convex option
C                       added interp to node mode
C
C 11 MAY 1993    J.Giles - edited by J Giles, BAe Filton
C                     -  changes to array dimensions  
C                     - mesh data only generated for testing 
C 12 May 1993    ianc - field mesh structure revised,
C                     - fixed declaration bugs.
C                     - debugged flow data also.
C----------------------------------------------------------------------
C
C	"@(#)read_ej65.f	7.1 AVS 92/02/25"
C			Copyright (c) 1990,1991 by
C			Advanced Visual Systems Inc.
C			All Rights Reserved
C	
C	This software comprises unpublished confidential information of
C	Advanced Visual Systems Inc. and may not be used, copied or made
C	available to anyone, except in accordance with the license
C	under which it is furnished.
C	
C     ************************************************************
C     AVS module description routine - defines the module,
C     its output port, and its parameters
C     ************************************************************
C
      subroutine AVSinit_modules
C   IAC CODE CHANGE :       include 'avs/avs.inc'
	INCLUDE '/usr/avs/include/avs.inc'
      external read_ej65_compute
      integer iparm
C     
C     Set the module name
C     
      call AVSset_module_name('Read EJ65', 'data')
C
      call AVSset_module_flags( single_arg_data )
C
C     Define the output port: an AVS two-dimensional field of 
C     reals, with a two-dimensional mapping space
C
      oport = AVScreate_output_port('field',
     *     'field 3D 7-vector real 3-space irregular')
C
C     This tells AVS to automatically free the output data after
C     each call.  This way, we can reallocate it each time without
C     difficulty.
C
      call AVSautofree_output(oport)
C     
C     Add the parameters, AVS will automatically add the
C     user interface controls:
C
      iparm = AVSadd_parameter('EJ65 Mesh File',
     *  'string', ' ',' ',' ')
      call AVSconnect_widget(iparm,'browser')
C
      iparm = AVSadd_parameter('EJ65 FlowField File',
     *  'string', ' ',' ',' ')
      call AVSconnect_widget(iparm,'browser')
C
      iparm = AVSadd_parameter('Convex Non-IEEE', 'boolean',0,0,1)
C
      iparm = AVSadd_parameter('Freeze', 'boolean',1,0,1)
C     
      iparm = AVSadd_parameter('Interp to Nodes', 'boolean',0,0,1)
C     
C     Set the compute procedure
C     
      call AVSset_compute_proc(read_ej65_compute)
      
      return
      end
      
C ************************************************************
C read_ej65_compute function - creates the field based on the
C state of the parameters (described above).
C
C This module generates a 3D field, so the compute function expects
C an output field "pointer" value which it uses as an argument to the
C field accessor and allocation functions. The accessor functions
C expect to be told that this is an output field since it will handle
C output field pointers slightly differently than input field pointers.
C
C Through AVSdata_alloc and the accessor functions, this function will
C fill the data and points arrays with its computation.
C
C The parameters are
C passed to this subroutine by the AVS flow exective based on
C the settings of the on-screen control widgets.
C ************************************************************

      integer function read_ej65_compute(output_field,
     *    mesh_filename, flow_filename, convex, freeze, interp )

C   IAC CODE CHANGE :       include 'avs/avs.inc'
	INCLUDE '/usr/avs/include/avs.inc'

C The input (none), output and parameter arguments
C   The output field is received as an integer
C
      integer output_field
      integer convex, freeze, interp
      character*128 mesh_filename 
      character*128 flow_filename 
C
C
C Internal declarations
C
C====================== USER CONFIGURABLE DIMENSIONS ==========
C  NOTE: Set these parameters to the largest i-j mesh size
C        you expect to use, but not much bigger, 
C        because they take up static memory.
C
      parameter (IMAX = 600)
      parameter (JMAX = 600)
C
C==============================================================
C
      external read_ej65_compute2
      external read_ej65_compute3
      integer dims, iresult, read_ej65_compute2
      integer read_ej65_compute3
      integer ofield, pfield
      integer ocoords, pcoords
      integer itemp
      real field_data
      dimension field_data(1), dims(3)
      dimension coords(1)
      character*64 field_descriptor
      character*128 errbuf
      integer ni, nj, nk
      double precision mesh (3, IMAX, JMAX )
C change J Giles, 11 May 1993 
C     real   w( 7, IMAX, JMAX ) 
      real   w( 7, IMAX+1, JMAX+1)   
      save   dims, ofield, pfield, ni, nj, nk
      save   ofield, pfield, ocoords, pcoords
      save   coords, mesh, w
C
C define the logical unit for the input file
C
      integer IN_UNIT
      save    IN_UNIT
      data    IN_UNIT  / 18 /
C
C Open the Input file, attatch to logical unit, and read
C the first header line for dimensions to follow
C
      if ( mesh_filename(1:1) .eq. ' ' .or. mesh_filename .eq. '') then
         read_ej65_compute = 0
         return
      endif
      if ( flow_filename(1:1) .eq. ' ' .or. flow_filename .eq. '') then
         read_ej65_compute = 0
         return
      endif
C
C Check if freeze is on, if so, exit, without processing anything.
C
      if ( freeze .eq. 1 ) then
         read_ej65_compute = 0
         return
      endif

C
C Open the file first, containing the mesh coordinates
C
      call open_ej65  ( mesh_filename, IN_UNIT )
C
C Get header dimensions from start of mesh file
C
      call get_header ( IN_UNIT, ni, nj, nk, convex )
C
C Determine what kind of field is being requested and provide that 
C information in the descriptor used to AVSdata_alloc. Once set by 
C AVSdata_alloc the individual field scalar descriptors (uniform, type,
C etc) cannot be changed
C
       field_descriptor = 'field 3D 7-vector real 3-space irregular'
C     
C  Setup dimensions array used in AVSdata_alloc
C     AVSdata_alloc( field_descriptor, dimensions_array)
C
      dims(1) = ni
      dims(2) = nj
      dims(3) = nk
C
C  AVSdata_alloc will automatically initialize the field variables
C  covering the dimensions, type, etc.       
C
      output_field = AVSdata_alloc(field_descriptor, dims)
      if (output_field .eq. 0) then 
         write(errbuf,90) field_descriptor
 90      format('Error allocating field ',A)
         call AVSerror(errbuf)
         read_ej65_compute = 0
         return
      endif
C
C  Example of accessors for field scalar values - itemp is value
C
C      itemp = AVSfield_get_int(output_field, AVS_field_ndim)
C      itemp = AVSfield_get_int(output_field, AVS_field_nspace)
C      itemp = AVSfield_get_int(output_field, AVS_field_veclen)
C      itemp = AVSfield_get_int(output_field, AVS_field_type)
C      itemp = AVSfield_get_int(output_field, AVS_field_size)
C      itemp = AVSfield_get_int(output_field, AVS_field_uniform)      
C
C  The internal arrays in the field have been allocated in 
C  AVSdata_alloc so all we have to do is to retrieve the offsets to 
C  the data array and the coordinates array. This works just like the
C  AVSptr_offset function.
C
      iresult = AVSfield_data_offset(output_field, field_data, ofield)
      iresult = AVSfield_points_offset(output_field, coords, ocoords)
C
C  Now call a second function which directly manipulates the field
C  arrays - passing the first actual element of the array 
C  allocations allows the second routine
C  to receive them as adjustable arrays.
C  The first one will access the mesh file, and fill the coord array
C
      write(*,*) 'Read EJ65: Reading Mesh File...'
C      
      iresult=read_ej65_compute2(field_data (ofield +1),
     *                           coords     (ocoords+1),
     *                           mesh,
     *                           ni, nj, nk, convex,
     *                           IN_UNIT )
C
C     Close the mesh file here
      close ( IN_UNIT )
C
C  Open the Flowfield solution file,
C
      call open_ej65  ( flow_filename, IN_UNIT )
C
C  Now call a third function which directly manipulates the field
C  arrays - passing the first actual element of the array 
C  allocations allows the third routine
C  to receive them as adjustable arrays.
C  This one will access the FLOWFIELD file, and fill the field_data array
C 
C
C      
      iresult=read_ej65_compute3(field_data (ofield +1),
     *                           coords     (ocoords+1),
     *                           w,
     *                           ni, nj, nk, convex, 
     *                           IN_UNIT )
 
C      Close the flow file here
       close ( IN_UNIT )
C
C remove interp function, J Giles 11 May 1993.
C
C If interpolation to nodal positions is needed, pass to interp routine

      if ( interp .eq. 1 ) then
          call interp_to_nodes(  field_data, coords, ni, nj, nk )
      endif
C
C Set the field attribute label strings, for downstream annotation
C
      call AVSfield_set_labels (output_field,
     * 'density#x-mom#y-mom#z-mom#total-energy#pressure#dt/vol',
     * '#' )
C
C Report results of operation
C
      write(*,*) 'Read EJ65: AVS Field Complete.'
      if ( iresult .ne. 1 ) 
     &    write(*,*) 'Read EJ65: error exit.'
      read_ej65_compute = iresult
      return


      end

C *******************************************************************
C This function fills the Output port COORDINATE section with values
C from the MESH file, already assumed openned, and positioned
C past the header.
C *******************************************************************
      integer function read_ej65_compute2(field_data, coords,
     *     mesh,
     *     ni, nj, nk, convex, 
     *     IN_UNIT )
      
C   IAC CODE CHANGE :       include 'avs/avs.inc'
	INCLUDE '/usr/avs/include/avs.inc'
C
C     Pointers to the field and the mapping coordinates
C
C
      integer      ni, nj, nk
      integer      i, j, k
      integer      ii, jj, kk,m
      real         field_data (7, ni, nj, nk)
      real         coords (ni, nj, nk, 3)
      integer      convex
      integer      IN_UNIT
      integer      pp
      real         correction_factor
      real         temp_a, temp_b, temp_c
C
C  array size changed, J Giles, 11 May 1993 
      double precision mesh (3, ni+1, nj+1 )
C       double precision mesh (3, ni, nj ) 
C
C     Correct for Convex Non-IEEE floating point format, if needed
      if ( convex .eq. 1 ) then
          correction_factor = 0.25
      else
          correction_factor = 1.0
      endif
C
C     Load in the field data
C     loop over each k slice, bring it in as doubles, then
C     convert to real in the output field structure
C     the inner loop, over x-y-z is unrolled for speed
C
      do k=1,nk
          read ( IN_UNIT ) mesh
C          read ( IN_UNIT, * )
C     &       (((mesh(m,ii,jj),m=1,3),ii=1,ni+1),jj=1,nj+1)
          do j=1,nj
              do i=1,ni
                  coords(i,j,k,1) = mesh (1,i,j) * correction_factor
                  coords(i,j,k,2) = mesh (2,i,j) * correction_factor
                  coords(i,j,k,3) = mesh (3,i,j) * correction_factor
C
C                  temp_a = coords(i,j,k,1)
C                  temp_b = coords(i,j,k,2)
C                  temp_c = coords(i,j,k,3)
C                  if ( i .eq. 1 .and. j .eq. 1 .and. k .eq. 1 ) 
C     &             write(*,*) temp_a, temp_b, temp_c


              enddo
          enddo
      enddo

C
C     Return a 1 to indicate success
C
      read_ej65_compute2 = 1

      return
      end


C *******************************************************************
C This function fills the Output port DATA section with values
C from the FLOWFIELD file, already assumed openned, and positioned
C at the BEGINNING of the file.
C *******************************************************************
      integer function read_ej65_compute3(field_data, coords,
     *     w,
     *     ni, nj, nk, convex, 
     *     IN_UNIT )
      
C   IAC CODE CHANGE :       include 'avs/avs.inc'
	INCLUDE '/usr/avs/include/avs.inc'
C
C     Pointers to the field and the mapping coordinates
C
      integer      ni, nj, nk
      integer      i, j, k
      real         field_data (7, ni, nj, nk)
      real         coords (ni, nj, nk, 3)
      integer      convex
      integer      IN_UNIT
      integer      vector, veclen
C
      double precision title(20), header(32)
      double precision cfl, hm, qfil, vis2, vis4, bc, mach
      double precision alpha, cl, cd, cd0, cm
      integer local_ni, local_nj, local_nk
      real         correction_factor
C
C  Local Temp Storage for field, one slice at a time
C  Note critical dimension adjustment!
C
C  change array dimensions, J Giles 11 May 1993
      real  w( 7, ni+1, nj+1 )
C     real  w( 7, ni-1, nj-1 ) 
C
C     vector length is fixed by the file and field formats:
      veclen = 7
C
C     Correct for Convex Non-IEEE floating point format, if needed
      if ( convex .eq. 1 ) then
          correction_factor = 0.25
      else
          correction_factor = 1.0
      endif
C 
C     Read the Header Section of the Flow Field
C
      read ( IN_UNIT ) title
      read ( IN_UNIT ) header
C
C     Extract some key parameters, for later use (?)
C
      cfl   = header ( 1) * correction_factor
      hm    = header ( 2) * correction_factor
      qfil  = header ( 3) * correction_factor
      vis2  = header ( 4) * correction_factor
      vis4  = header ( 5) * correction_factor
      bc    = header ( 6) * correction_factor
      mach  = header ( 7) * correction_factor
      alpha = header ( 8) * correction_factor
      cl    = header ( 9) * correction_factor
      cd    = header (10) * correction_factor
      cd0   = header (11) * correction_factor
      cm    = header (12) * correction_factor
C
      local_ni = header (16) * correction_factor
      local_nj = header (17) * correction_factor
      local_nk = header (18) * correction_factor
C
       write(*,*) ' read in header ok for flow data'
C     The mesh is one size bigger in each dimension than
C     the FlowField file, since the flow data is assumed
C     to be centered in the middle of each cell,
C     so the local dims should be one less than the mesh size,
C     lets check it to make sure.
C
C change test, J Giles, 11 May1993
C
C     if ( ni .ne. local_ni+1 .or. 
C    *     nj .ne. local_nj+1 .or.
C    *     nk .ne. local_nk+1 ) then
C 
      if ( ni .ne. local_ni-1 .or.
     *     nj .ne. local_nj-1 .or.
     *     nk .ne. local_nk-1 ) then
          write(*,*) 'Read EJ65: mesh and flow dimensions do not match'
          write(*,*) ' Mesh Res: ',ni,nj,nk
          write(*,*) ' Flow Res: ',local_ni,local_nj,local_nk
          read_ej65_compute3 = 0
          return
      endif
       write(*,*) ' dims ok for flow data'
C
C     Pre-Zero the output field array, so we can average,
C     or partially fill, or whatever, and not leave
C     garbage values.
C
      do k=1, nk
          do j=1, nj
              do i=1, ni
                  do vector=1, veclen
                      field_data(vector,i,j,k) = 0.0
                  enddo
              enddo
          enddo
      enddo
C
C     Load in the field data
C     loop over K slices, convert from cell center to corners
C     Each slice is read into a temporary array, w, then
C     copied with the correct dimension adjustments to
C     the AVS field, missing out the last row and column
C
C changes, J Giles, BAe, 11 May 1993  
C
      do k=2, nk    
C     do k=2, local_nk
C
C         Get One slice from file, full speed i/o binary access
          read ( IN_UNIT ) w
C
C         Copy over to field data structure
C changes, J Giles
          do j=1, nj    
C         do j=1, local_nj
              do i=1, ni
C             do i=1, local_ni
                  do vector=1, veclen
                      field_data(vector,i,j,k) = w(vector,i,j)
     &                    * correction_factor
                  enddo
              enddo
          enddo
      enddo
C
       write(*,*) ' ok for flow data, good exit'
C
C
C     Return a 1 to indicate success
C
      read_ej65_compute3 = 1

      return
      end

C *******************************************************************
C routine to Interpolate the flowfield onto the mesh node points
C *******************************************************************
      subroutine interp_to_nodes(  field_data, coords,
     *     ni, nj, nk )
C
C   IAC CODE CHANGE :       include 'avs/avs.inc'
	INCLUDE '/usr/avs/include/avs.inc'
C
C     Pointers to the field and the mapping coordinates
C
      integer      ni, nj, nk
      integer      i, j, k
      real         field_data (7, ni, nj, nk)
      real         coords (ni, nj, nk, 3)
      integer      vector, veclen
      real         fraction
C
      veclen = 7
      fraction = 1.0 / 8.0
C
C  Average interior data from cell centers to node corners
C
      do k=nk, 2, -1
          do j=nj, 2, -1
              do i=ni, 2, -1
                  do vector=1, veclen
                      field_data(vector,  i,  j,  k) = 
     &                field_data(vector,  i,  j,  k) * fraction +
     &                field_data(vector,i-1,  j,  k) * fraction +
     &                field_data(vector,  i,j-1,  k) * fraction +
     &                field_data(vector,i-1,j-1,  k) * fraction +
     &                field_data(vector,  i,  j,k-1) * fraction +
     &                field_data(vector,i-1,  j,k-1) * fraction +
     &                field_data(vector,  i,j-1,k-1) * fraction +
     &                field_data(vector,i-1,j-1,k-1) * fraction 
                  enddo
              enddo
          enddo
      enddo


      return
      end


C *******************************************************************
C This function OPENS the mesh or solution file
C *******************************************************************
      subroutine open_ej65 ( filename, IN_UNIT )

      character*128 filename 
      integer      IN_UNIT
      integer      ier
C
C  Open the BINARY Geometry or FlowField File for reading
C
C============= USER CONFIGURATION AREA: =========================
C  Note: Keywords here may vary on different machines!
C
      open( unit   = IN_UNIT,
     &      file   = filename,
     &      status = 'OLD',
     &      form   = 'UNFORMATTED',
     &      iostat = ier )
C
C
C================================================================
C
C
      rewind ( IN_UNIT )
C
C
      return
      end
      
C *******************************************************************
C This function extracts key header parameters from the mesh file
C *******************************************************************
      subroutine get_header ( IN_UNIT, ni, nj, nk, convex )

      integer      IN_UNIT
      integer      ni, nj, nk
      integer      ier
      integer      convex
C
      double precision gtitle(20), ghead(32)
      integer ktip, nls, nus,jb, nwake
      integer itl, itu, nm, nw
      double precision xmin, xmax, ymin, ymax, zmin, zmax
      double precision correction_factor
C
C     Correct for Convex Non-IEEE floating point format, if needed
      if ( convex .eq. 1 ) then
          correction_factor = 0.25
      else
          correction_factor = 1.0
      endif
C
C Pull in the header into local arrays:
C
C      read ( IN_UNIT, '(20a4)' ) gtitle
C      read ( IN_UNIT, * ) (ghead(i),i=1,32)
C
      read ( IN_UNIT ) gtitle
      read ( IN_UNIT ) ghead
C
C Extract the dimensions, and other derived header quantities:
C Only array sizes are passed back as yet, others may be
C used later.
C
      ni   = ghead ( 1)* correction_factor
      nj   = ghead ( 2)* correction_factor
      nk   = ghead ( 3)* correction_factor
      ktip = ghead ( 4)* correction_factor
      nls  = ghead ( 5)* correction_factor
      nus  = ghead ( 6)* correction_factor
      jb   = ghead ( 7)* correction_factor
      nwake= (ni - (nls + nus) + 1 ) / 2
      itl  = nwake + 1
      itu  = ni - nwake
      nm   = itu - itl + 2
      nw   = itl - 1
C


      xmin = ghead (11)* correction_factor
      xmax = ghead (12)* correction_factor
      ymin = ghead (13)* correction_factor
      ymax = ghead (14)* correction_factor
      zmin = ghead (15)* correction_factor
      zmax = ghead (16)* correction_factor
C
      return
      end