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	"@(#)phoenics_read.f"
C
C    This module reads in PHOENICS results files (sequential access)
C    and creates labelled AVS fields from the information contained
C    in them.  If a BFC case is called for, a file browser widget is
C    called to read in the sequential access XYZ file.
C
C    Substantial portions of this code were originally written as a
C    filter (translator) to convert PHOENICS results files into 
C    MPGS readable files on a Cray system.  The interpolation routine
C    comes from this source.
C
C    Written by Michael Rangitsch,   July, 1991
C
C    version 2.0                     September 23, 1991
C       this version allows the choice of interpolation of the grid to
C       the cell center locations (see version 1.0) with the velocities
C       assumed to be at the cell centers and to an interpolated version
C       of the velocity and scalar field values.  the values are inter-
C       polated to the nodal locations using the areas (or projected
C       areas) for the velocity components and the cell volumes for the
C       scalar fields.  the velocity at a node is considered to be 
C       given by multiplying the velocity at a face by the area of the
C       face and summing up these values for the four faces adjacent to
C       a nodal point and then dividing by the total area.  The boundary
C       values are found by extrapolating the two values adjacent to the
C       boundary.   Scalar values (cell centered) are found similarly by 
C       using the volumes of the four cells adjacent to a nodal point.
C	
C    version 1.0                     September 1,1991
C       this version interpolates the grid information to the
C       cell center, i.e. the data values are assumed to be cell
C       centered (even the velocities).  The u1 velocity at the
C       high x plane, the v1 velocity at the high y plane and the
C       w1 velocity at the high z plane are all assumed to be the
C       same velocity as at the previous plane
C
C       known problems:    BFC velocities are not handled correctly.
C                          Cartesian and Cylindrical velocities are correct,
C                          although there are small problems at R=0 in
C                          cylindrical grids
C
C
C     ************************************************************
C     AVS module description routine - defines the module,
C     its output port, and its parameters
C     ************************************************************
C
      subroutine AVSinit_modules
      include '/usr/avs/include/avs.inc'
      external phoenics_field_compute
      integer iparm
C     
C     Set the module name
C     
      call AVSset_module_name('phoenics interp read', 'data')
C
C     Define the output port: an AVS three-dimensional field of 
C     reals, with a three-dimensional mapping space
C
      oport = AVScreate_output_port('field',
     *     'field 3D real 3-space')

C
C     Define the geometry output port (not functional yet!)
C
C     oport1 = AVScreate_output_port('geometry','geom')
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:
C
C        the first is the file browser for the phi file
C        the second is the file browser for the xyz file
C        the third is the toggle switch to produce an interpolated field
C
      iparm1 = AVSadd_parameter('PHI file browser','string',
     *    ' ',' ',' ')
      call AVSconnect_widget(iparm1,'browser')
C
      iparm2 = AVSadd_parameter('XYZ file browser','string',
     *    ' ',' ',' ')
      call AVSconnect_widget(iparm2,'browser')
C
c      iparm3 = AVSadd_parameter('Interpolate values','boolean',0,1,1)
c      call AVSconnect_widget(iparm3,'toggle')
C     
C     Set the compute procedure
C     
      call AVSset_compute_proc(phoenics_field_compute)
      
      return
      end
      
C     ************************************************************
C     phoenics_field_compute function - creates the field based on the
C     state of the parameters (described above).
C
C     ************************************************************
      integer function phoenics_field_compute(pfield,nxl,nyl,nzl,nval,
     *    coordflag,nspace,pcoords,phiname,xyzname)
C     ,interp)
      
      include '/usr/avs/include/avs.inc'
C
C     Pointers to the field and the mapping coordinates
C
      integer pfield, pcoords,coordflag,ofield,ocoords,gfield
      integer nspace,interp,pax,oax,pay,oay,paz,oaz,pvols,ovols

      external phoenics_field_compute2,phoenics_cart_read,
     *    phoenics_bfc_read,phoenics_cyl_conv,phoenics_data_read,
     *    phoenics_bfc_reorder,phoenics_cart_areas,phoenics_cart_vols,
     *    phoenics_cyl_areas,phoenics_cyl_vols,phoenics_bfc_areas,
     *    phoenics_bfc_vols,uvel_extrap,vvel_extrap,wvel_extrap
     *    uvel_interp,vvel_interp,wvel_interp,oth_interp,field_reorder
      integer phoenics_field_compute2,phoenics_cart_read,
     *    phoenics_bfc_read,phoenics_cyl_conv,phoenics_data_read,
     *    phoenics_bfc_reorder,phoenics_cart_areas,phoenics_cart_vols,
     *    phoenics_cyl_areas,phoenics_cyl_vols,phoenics_bfc_areas,
     *    phoenics_bfc_vols,uvel_extrap,vvel_extrap,wvel_extrap
     *    uvel_interp,vvel_interp,wvel_interp,oth_interp,field_reorder
      integer iresult
C
      parameter (nphimax = 50, nxmax=100, nymax=100, nzmax=100)
      integer stored,l
      real field(1),coords(1)
      real ax(0:nxmax,nymax,nzmax),ay(nxmax,0:nymax,nzmax),
     *   az(nxmax,nymax,0:nzmax),vols(nxmax,nymax,nzmax)
c     real ax(1),ay(1),az(1),vols(1)
      character*(*) phiname,xyzname
C 
C    declare the arrays for PHOENICS input (the header arrays)
C
      parameter (iphi = 9)
      character*40 title
      character*4  name(nphimax+2),name_stored(nphimax+2)
      character*1  cartes,onephs,bfc,xcycle,store(nphimax+2)
      character*250 labels
      real xc(0:nxmax),yc(0:nymax),zc(0:nzmax),pcor(nzmax)
      integer nx,ny,nz,nphi,den1,den2,epor,npor,hpor,vpor,lenrec
      real rinner
      integer nv,nx1,ny1,nz1
C
C   if no file name is specified, return
C
      print*,' in compute function'
      if(phiname .eq. '') then
       phoenics_field_compute = 0
       return
      endif
C
C   if the file name is specified, open it and read the header
C
      open(unit=iphi,file=phiname,status='OLD',iostat=ier)
      rewind(iphi)
C
      READ(iphi,5) TITLE
 5    format(a40)
      READ(iphi,6) CARTES,ONEPHS,BFC,XCYCLE
 6    format(1x,4(A1))
      READ(iphi,*) NX,NY,NZ,NPHI,DEN1,DEN2,EPOR,NPOR,HPOR,VPOR,LENREC
      READ(iphi,*) RINNER
c     print*, rinner
      rinner = max(rinner,1.e-10)
      READ(iphi,7) (NAME(I),I=1,19)
      READ(iphi,7) (NAME(I),I=20,38)
      READ(iphi,7) (NAME(I),I=39,50)
 7    format(1x,19(A4))
      READ(iphi,9) (XC(I),I=1,NX)
      READ(iphi,9) (YC(I),I=1,NY)
      READ(iphi,9) (ZC(I),I=1,NZ)
      READ(iphi,9) (PCOR(I),I=1,NZ)
 9    format(6(1pe13.6))
      READ(iphi,8) (STORE(I),I=1,NPHI)
 8    format(1x,50(A1)) 
C
C  check that an xyz file has been chosen for bfc cases
C
      if(bfc .eq. 'T') then
        if(xyzname .eq. '') then
          phoenics_field_compute = 0
          return
        endif
      endif
C
C  now step through name array and store array and build stored name
C  array
C
      nstored = 1
      do i = 1,nphi
         if(store(i) .eq. 'T') then
             name_stored(nstored) = name(i)
             nstored = nstored + 1
         endif
      enddo
C
C  set up the proper sizes to allow for space allocation for the field 
C  array
C
      ncoord=3
      nv=nstored-1
      nval = nv
C
C   add room for another velocity component if the problem is 2 dimensional
C   if we have a cylindrical or cartesian case, add the appropriate velocity
C   component, if we have a bfc case at the appropriate cartesian 
C   component.  If we are dealing with 2 phases, add the components for 
C   both phases
C
      if(nx .eq. 1) then
       if(onephs .eq. 'T') then
          nphi = nphi + 1
          nval = nv + 1
          store(nphi) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval) = 'UCRT'
          else
            name_stored(nval) = 'U1  '
          endif
        else
          nphi = nphi + 2
          nval = nv + 2
          store(nphi) = 'T'
          store(nphi-1) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval-1) = 'U2CR'
            name_stored(nval) = 'UCRT'
          else
            name_stored(nval-1) = 'U2  '
            name_stored(nval) = 'U1  '
          endif
        endif
       endif
C
      if(ny .eq. 1) then
       if(onephs .eq. 'T') then
          nphi = nphi + 1
          nval = nv + 1
          store(nphi) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval) = 'VCRT'
          else
            name_stored(nval) = 'V1  '
          endif
        else
          nphi = nphi + 2
          nval = nv + 2
          store(nphi) = 'T'
          store(nphi-1) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval-1) = 'V2CR'
            name_stored(nval) = 'VCRT'
          else
            name_stored(nval-1) = 'V2  '
            name_stored(nval) = 'V1  '
          endif
        endif
       endif
C
      if(nz .eq. 1) then
       if(onephs .eq. 'T') then
          nphi = nphi + 1
          nval = nv + 1
          store(nphi) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval) = 'WCRT'
          else
            name_stored(nval) = 'W1  '
          endif
        else
          nphi = nphi + 2
          nval = nv + 2
          store(nphi) = 'T'
          store(nphi-1) = 'T'
          if(bfc .eq. 'T') then
            name_stored(nval-1) = 'W2CR'
            name_stored(nval) = 'WCRT'
          else
            name_stored(nval-1) = 'W2  '
            name_stored(nval) = 'W1  '
          endif
        endif
       endif
C
C   now set the indices for the data read statements
C           
c      nxf = 1
c      nyf = 1
c      nzf = 1
c      nxl = nx
c      nyl = ny
c      nzl = nz
C
C  if interpolation is called for, then increase the data array size
C  by 1 in each dimension and start the read counters at 2 to allow
C  for interpolation to the first node position
C
c      if(interp .eq. 1) then
        nxl = nx + 1
        nyl = ny + 1
        nzl = nz + 1
        nxf = 2
        nyf = 2
        nzf = 2
c      endif
c     print*,nxl,nyl,nzl,pfield,ofield
C
C  allocate space for the field array, making sure to zero out the
C  data values so the add
C  first free the data field if interpolation is called
C  for
C
c      if(interp .eq. 1) then
c         call AVSptr_free('field',pfield,ofield)
c         call AVSptr_free('field',pcoords,ocoords)
c         print*, 'freeing ',ires
c      endif
c        
      iresult = AVSptr_alloc('field', nxl*nyl*nzl*nval, 4, 1, field,
     *     pfield, ofield)
      print*, 'allocating',iresult,' nv ',nv,' nval ',nval
C
C  call the read function with the proper indices to allow for 
C  interpolation space
C
      iresult = phoenics_data_read(field(ofield+1),nxf,nxl,nyf,nyl,
     *                                 nzf,nzl,nv,nval,iphi)
C
C   if interpolation is called for allocate space for the face area
C   arrays and the cell volume array
C
c      if(interp .eq. 1) then
c       iresult = AVSptr_alloc(' ', nxl*ny*nz, 4, 1, ax, 0, oax)
c       iresult = AVSptr_alloc(' ', nx*nyl*nz, 4, 1, ay, 0, oay)
c       iresult = AVSptr_alloc(' ', nx*ny*nzl, 4, 1, az, 0, oaz)
c       iresult = AVSptr_alloc(' ', nx*ny*nz, 4, 1, vols, 0 ,ovols)
c      endif
C
C     Since we are allocating a field from scratch, we must allocate
C     space for mapping coordinates, for non bfc cases we must allow
C     non-even grid spacing, for bfc cases we must allow for storage of
C     general x,y,z triplet location information
C
C
      if(bfc .eq. 'F') then
C
C  cartesian case
C
        if(cartes .eq. 'T') then
          coordflag = RECTILINEAR
          iresult = AVSptr_alloc('field', nxl+nyl+nzl, 4, 0, 
     *         coords, pcoords, ocoords)
          iresult = phoenics_cart_read(coords(ocoords+1),nxl,nyl,nzl,
     *                  xc,yc,zc,nx,ny,nz)
C
C  if the interpolation flag is set, compute the required facial areas
C  and volumes
C
           iresult = phoenics_cart_areas(xc,yc,zc,nx,ny,nz,ax,ay,az)
c    *       ax(oax+1),
c    *       ay(oay+1),az(oaz+1))
c    
           iresult = phoenics_cart_vols(xc,yc,zc,nx,ny,nz,vols)
c    *   (ovols+1))
        else
C 
C  cylindrical case
C
          coordflag = IRREGULAR
          iresult = AVSptr_alloc('field', nxl*nyl*nzl*ncoord, 
     *     4, 0, coords, pcoords, ocoords)
          iresult = phoenics_cyl_coord(coords(ocoords+1),
     *         nxl,nyl,nzl,ncoord,xc,yc,zc,nx,ny,nz,rinner)
C
C  if the interpolation flag is set, compute the required facial areas
C  and volumes
C
         iresult = phoenics_cyl_areas(xc,yc,zc,nx,ny,nz,ax,ay,az,rinner)
c    *     x(oax+1),
c    *     ay(oay+1),az(oaz+1),rinner)
         iresult = phoenics_cyl_vols(xc,yc,zc,nx,ny,nz,vols,rinner)
c    *     vols(ovols+1),rinner)
        endif
      else
         coordflag = IRREGULAR
         iresult = AVSptr_alloc('field', nxl*nyl*nzl*ncoord,
     *    4, 0, coords, pcoords, ocoords)
         iresult = phoenics_bfc_read(coords(ocoords+1),nxl,nyl,nzl,
     *                        ncoord, xyzname)
         if(iresult .eq. 0) then
              phoenics_field_compute = 0
              return
         endif
         iresult = phoenics_bfc_areas(xyzname,nx,ny,nz,ax,ay,az)
c    *        ax(oax+1),
c    *        ay(oay+1),az(oaz+1))
         iresult = phoenics_bfc_vols(xyzname,nx,ny,nz,vols)
c    *    (ovols+1))
      endif
C
C   for interpolation, extrapolate the velocity values to the low
C   faces in each direction and interpolate values to the nodal points
C
       istored = 0
       do i = 1,nphi
        if(store(i) .eq. 'T') then
         istored = istored + 1
         if(i .eq. 3 .or. i .eq. 4) then
            if(bfc .eq. 'F' .and. cartes .eq. 'T') then
              iresult = cart_uvel_extrap(field(ofield+1),istored,nxl,
     *                                   nyl,nzl,nval,xc,nx)
            else
             iresult = uvel_extrap(field(ofield+1),istored,nxl,nyl,nzl,
     *                            nval,coords(ocoords+1),ncoord)
            endif
            iresult = uvel_interp(field(ofield+1),istored,nxl,nyl,nzl,
     *                         nval,ax,nx,ny,nz)
c    *     (oax+1),nx,ny,nz)
         else
          if(i .eq. 4 .or. i .eq. 5) then
           if(bfc .eq. 'F' .and. cartes .eq. 'T') then
            iresult = cart_vvel_extrap(field(ofield+1),istored,nxl,nyl
     *                                 nzl,nval,yc,ny)
           else
            iresult = vvel_extrap(field(ofield+1),istored,nxl,nyl,nzl,
     *                            nval,coords(ocoords+1),ncoord)
           endif
            iresult = vvel_interp(field(ofield+1),istored,nxl,nyl,nzl,
     *                         nval,ay,nx,ny,nz)
c    *       (oay+1),nx,ny,nz)
          else
           if(i .eq. 6 .or. i .eq. 7) then
            if(bfc .eq. 'F' .and. cartes .eq. 'T') then
             iresult = cart_wvel_extrap(field(ofield+1),istored,nxl,nyl,
     *                                  nzl,nval,zc,nz)
            else
             iresult = wvel_extrap(field(ofield+1),istored,nxl,nyl,nzl,
     *                            nval,coords(ocoords+1),ncoord)
            endif
            iresult = wvel_interp(field(ofield+1),istored,nxl,nyl,nzl,
     *                         nval,az,nx,ny,nz)
c    *       (oaz+1),nx,ny,nz)
           else
            print*, nxl,nyl,nzl
            iresult = oth_interp(field(ofield+1),istored,nxl,nyl,nzl,
     *                nval,vols,nx,ny,nz)
C    *        (ovols+1),nx,ny,nz)
           endif
          endif
         endif
        endif
       enddo
c      endif
      print*,' done with interpolation functions'
C
C  if the data set was bfc, resort the velocity components such that
C  the order of the cartesian coordinates is ucrt, vcrt, wcrt or for
C  the two phase case, u2cr, v2cr, w2cr
C
C  first phase
C
      lucrt = 0
      lvcrt = 0
      lwcrt = 0
      if(bfc .eq. 'T') then
        do i = 1,nval
          print*, name_stored(i)
          if(name_stored(i) .eq. 'UCRT') lucrt = i
          if(name_stored(i) .eq. 'VCRT') lvcrt = i
          if(name_stored(i) .eq. 'WCRT') lwcrt = i
        enddo
        iresult = phoenics_bfc_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                   lucrt,lvcrt,lwcrt,name_stored,nphimax)
      endif
C
C  second phase
C
      lucrt = 0
      lvcrt = 0
      lwcrt = 0
      if(bfc .eq. 'T') then
        do i = 1,nval
          if(name_stored(i) .eq. 'U2CR') lucrt = i
          if(name_stored(i) .eq. 'V2CR') lvcrt = i
          if(name_stored(i) .eq. 'W2CR') lwcrt = i
        enddo
        iresult = phoenics_bfc_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                   lucrt,lvcrt,lwcrt,name_stored,nphimax)
      endif
C
C   now reorder the data values for non bfc cases to allow place the velocity
C   components in the proper locations in the field structure
C
      lfst = 0
      if(store(1) .eq. 'T') lfst = lfst + 1
      if(store(2) .eq. 'T') lfst = lfst + 1
 
      if(bfc .eq. 'F') then
        nrep = lfst
        if(nx .eq. 1) then
          nrep = nrep + 1
          store(3) = 'T'
          iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval)
          if(onephs .eq. 'F') then
           nrep = nrep + 1
           store(4) = 'T'
           iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval-1)
          endif
        endif
        if(ny .eq. 1) then
          nrep = nrep + 1
          store(5) = 'T'
          iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval)
          if(onephs .eq. 'F') then
           nrep = nrep + 1
           store(6) = 'T'
           iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval-1)
          endif 
        endif
        if(nz .eq. 1) then
          nrep = nrep + 1
          store(7) = 'T'
          iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval)
          if(onephs .eq. 'F') then
           nrep = nrep + 1
           store(8) = 'T'
           iresult = field_reorder(field(ofield+1),nxl,nyl,nzl,nval,
     *                         name_stored,nrep,nval-1)
          endif 
        endif
      endif

C
C  for the cylindrical case convert the theta, r and z information into
C  x y and z
C
      if (bfc .eq. 'F' .and. cartes .eq. 'F') then
          iresult = phoenics_cyl_conv(field(ofield+1),
     *         nxl,nyl,nzl,nval,xc,yc,zc,nx,ny,nz,rinner,
     *         store,nphi)
C                          ,interp)
      endif
c
c   set labels
c
      write(labels,104) (name_stored(l),';',l=1,nval)
 104  format(50(a4,a1))      
c
c   set the field labels
c
      call AVSfield_set_labels(AVSport_field('field'),labels,';')
C
c  now free up the memory allocated for the volumes and areas
C
c     if(interp .eq. 1) then
c        call AVSptr_free(' ', pax, oax)
c        call AVSptr_free(' ', pay, oay)
c        call AVSptr_free(' ', paz, oaz)
c        call AVSptr_free(' ', pvols,ovols)
c     endif

      phoenics_field_compute = iresult
      return
      end

C **********************************************************************
      integer function phoenics_data_read(field,nxf,nxl,nyf,nyl,
     *                                   nzf,nzl,nv,nval,iphi)
C
C   this function reads in the data from the phi file phiname
C
      include '/usr/avs/include/avs.inc'
      integer nxl,nyl,nzl,nxf,nyf,nzf,nv,iphi,i,j,k,l
      dimension field(nval,nxl,nyl,nzl)
C
C   now read in the data for each z slab
C
      print*, nxf,nxl,nyf,nyl,nzf,nzl,nv
      print*, ' in function data read'
C
      do k = nzf,nzl
C
C   now for each variable
C
        do l = 1,nv
C
C   read in ni*nj variables with j varying fastest
C
           read(iphi,9,end=10) ((field(l,i,j,k),j=nyf,nyl),i=nxf,nxl) 
        enddo
      enddo
 9    format(6(1pe13.6))
C  
C   set value for successful completion
C
c10   print*,field(nv,nxl,nyl,nzl)
c     print*,field(nv,1,1,1)
c     print*,field(nv,2,2,2)
 10   close(iphi)
      print*,i,j,k,l
      phoenics_data_read = 1
      return
      end

      
C **************************************************************
      integer function phoenics_cart_read(coords,nxl,nyl,nzl,xc,yc,zc,
     *           nx,ny,nz)
C                      ,interp)
C
C    this function uses the xc, yc and zc arrays to create the
C    coordinate array for the input AVS field.  
C
      include '/usr/avs/include/avs.inc'
      integer nxl,nyl,nzl,nx,ny,nz
      real coords(nxl+nyl+nzl),xc(0:nx),yc(0:ny),zc(0:nz)
      integer i,j,k,icoord
C,interp
C
C   load the first values of the cell center arrays based depending
C   on the cartes flag
C
      print*, 'in function cart read'
      print*, nxl,nyl,nzl,nx,ny,nz
      xc(0)=0.0
      yc(0)=0.0
      zc(0)=0.0
C
C  now load the coords array with 'cell center' positional data
C
      icoord = 0
C
C  if the interpolation flag is set, use each positional value,
C  otherwise average the coordinate locations
C
c      if(interp .eq. 1) then
       do i = 1,nxl
        icoord = icoord + 1
        coords(icoord) = xc(i-1)
       enddo
       do j = 1,nyl
        icoord = icoord + 1
        coords(icoord) = yc(j-1)
       enddo
       do k = 1,nzl
        icoord = icoord + 1
        coords(icoord) = zc(k-1)
       enddo
c      else
c       do i = 1,nxl
c        icoord = icoord+1
c        coords(icoord) = (xc(i-1) + xc(i))/2.0
c       enddo
c       do j = 1,nyl
c        icoord = icoord+1
c        coords(icoord) = (yc(j-1) + yc(j))/2.0
c       enddo
c       do k = 1,nzl
c        icoord = icoord+1
c        coords(icoord) = (zc(k-1) + zc(k))/2.0
c       enddo
c      endif
C
C  set a return value
C
      phoenics_cart_read = 1
      return
      end

C *****************************************************************
      integer function phoenics_cyl_conv(field,nxl,nyl,nzl,nv,
     *         xc,yc,zc,nx,ny,nz,rinner,store,nphi)
C,interp)
C
C   this function converts the input cell face locations (interpreted
C   as theta,r and z values) into cell center cartesian values and
C   converts the theta and r velocity values (u1,u2 and v1,v2) into
C   cartesian velocity resolutes or loads the coordinate array with
C   unaveraged values, based on the interp flag
C
      include '/usr/avs/include/avs.inc'
      integer nx,ny,nz,nphi,lu1,lu2,lv1,lv2,lfst
      dimension field(nv,nxl,nyl,nzl)
      dimension xc(0:nx),yc(0:ny),zc(0:nz)
      real rinner
      character*1 store(nphi)
C 
C   first locate the velocity indices in the data array
C
      print*, 'in function cyl read'
      print*, nxl,nyl,nzl,nx,ny,nz
      lu1 = 0
      lu2 = 0
      lv1 = 0
      lv2 = 0
      lfst = 0
      if(store(1) .eq. 'T') lfst = lfst + 1
      if(store(2) .eq. 'T') lfst = lfst + 1
      if(store(3) .eq. 'T') lu1 = 1
      if(store(4) .eq. 'T') lu2 = 1
      if(store(5) .eq. 'T') lv1 = 1
      if(store(6) .eq. 'T') lv2 = 1
c     print*, lfst,lu1,lu2,lv1,lv2
c     print*, store(1),store(2),store(3),store(4),store(5),store(6)
C 
C  use these flag values to zero out components that don't exist
C  and convert those that do into cartesian components, since
C  we know that in the phi file u1 is lfst+1 depending on storage
C  of the p1 and p2 fields and v1 is one or two storge postions higher
C  depending on storage of u2 we can use these values to help in the
C  computation and storage of the cartesian components of the velocity
C
      if(store(3) .eq. 'T' .or. store(5) .eq. 'T') then
c     print*,lu1,lu2,lv1,lv2
      print*,1,nyl,nzl
      print*,field(lfst+lu1,1,nyl,nzl)
      print*,field(lfst+lu1+lu2+lv1,1,nyl,nzl)
      do k= 1,nzl
       do j= 1,nyl
        do i = 1,nxl
c         if(interp .eq. 1) then
          theta = xc(i-1)
c         else
c          theta=(xc(i)+xc(i-1))/2.0 
c         endif
          ft1 =-field(lfst+lu1,i,j,k)*sin(theta)*float(lu1) +
     *         field(lfst+lu1+lu2+lv1,i,j,k)*cos(theta)*float(lv1)
          ft3 = field(lfst+lu1,i,j,k)*cos(theta)*float(lu1) +
     *         field(lfst+lu1+lu2+lv1,i,j,k)*sin(theta)*float(lv1) 
          if(lu1 .ne. 0) field(lfst+lu1,i,j,k) = ft1
          if(lv1 .ne. 0) field(lfst+lu1+lu2+lv1,i,j,k) = ft3
         enddo
        enddo
       enddo
      endif
      print*,'converted'
      print*,field(lfst+lu1,1,nyl,nzl)
      print*,field(lfst+lu1+lu2+lv1,1,nyl,nzl)
      if(store(4) .eq. 'T' .or. store(6) .eq. 'T') then
       do k=1,nzl
        do j=1,nyl
         do i = 1,nxl
c           if(interp .eq. 1) then
            theta = xc(i-1)
c           else
c            theta=(xc(i)+xc(i-1))/2.0 
c           endif
           ft1 =-field(lfst+lu1+lu2,i,j,k)*sin(theta)*float(lu2) - 
     *         field(lfst+lu1+lu2+lv1+lv2,i,j,k)*cos(theta)*float(lv2)
           ft3 = field(lfst+lu1+lu2,i,j,k)*cos(theta)*float(lu2) +
     *         field(lfst+lu1+lu2+lv1+lv2,i,j,k)*sin(theta)*float(lv2) 
           if(lu2 .ne. 0) field(lfst+lu1+lu2,i,j,k) = ft1
           if(lv2 .ne. 0) field(lfst+lu1+lu2+lv1+lv2,i,j,k) = ft3
          enddo
         enddo
        enddo
       endif
C
C  now set the return value
C
      phoenics_cyl_conv = 1
      return
      end           
C *****************************************************************
      integer function phoenics_cyl_coord(coords,nxl,nyl,nzl,
     *         ncoord,xc,yc,zc,nx,ny,nz,rinner)
C,interp)
C
C   this function coorderts the input cell face locations (interpreted
C   as theta,r and z values) into cell center cartesian values 
C
	include '/usr/avs/include/avs.inc'
      integer nx,ny,nz
      dimension coords(nxl,nyl,nzl,ncoord)
      dimension xc(0:nx),yc(0:ny),zc(0:nz)
      real rinner
C 
C   first locate the velocity indices in the data array
C
      print*, 'in function cyl coord'
      print*, nxl,nyl,nzl,nx,ny,nz
C 
C  we now convert the coordinate information to cartesian values
C
      do i=1,nxl
c         if(interp .eq. 1) then
          theta = xc(i-1)
c         else
c          theta=(xc(i)+xc(i-1))/2.0 
c         endif
          do j=1,nyl
c            if(interp .eq. 1) then
              radius = rinner + yc(j-1)
c            else
c             radius= rinner + (yc(j-1)+yc(j))/2.0
c            endif           
            do k=1,nzl
               coords(i,j,k,1)=radius*cos(theta)
               coords(i,j,k,2)=radius*sin(theta)
c               if(interp .eq. 1) then
                coords(i,j,k,3)=zc(k-1)
c               else
c                coords(i,j,k,3)=(zc(k-1)+zc(k))/2.0
c               endif
            enddo
         enddo
      enddo   
C
C  now set the return value
C
      phoenics_cyl_coord = 1
      return
      end

C **********************************************************************
      integer function phoenics_bfc_read(coords,nxl,nyl,nzl,ncoord,
     *                                       xyzname)
c,interp)
C
C  this function reads in the coordinate array from the xyzname file
C  and converts it into cell centered postions by averaging adjacent
C  XC, YC and ZC values
C
	include '/usr/avs/include/avs.inc'
      parameter (nx1 = 101, ny1 = 101, nz1 = 101)
      integer nxl,nyl,nzl,ncoord,i,j,k
      dimension coords(nxl,nyl,nzl,ncoord)
      dimension xc(nx1,ny1,nz1), yc(nx1,ny1,nz1), zc(nx1,ny1,nz1)
      character*(*) xyzname
C
C  error exit if file name not specified
C
      print*, ' in function bfc read'
c     print*, nxl,nyl,nzl
      if(xyzname .eq. '') then
          phoenics_bfc_read = 0
          return
      endif
C
C  open the file
C  
      open(unit=28,file=xyzname,status='OLD')
      rewind(28)
      read(28,*) nxc,nyc,nzc
      if(nxc .ne. nxl .or. nyc .ne. nyl .or. nzc .ne. nzl) then
        pheonics_bfc_read = 0
        close(28)
        return
      endif
C
C  now read in the values 
C
      do k = 1,nzc
         read(28,15,end=17) ((xc(i,j,k),j=1,nyc),i=1,nxc)
         read(28,15,end=17) ((yc(i,j,k),j=1,nyc),i=1,nxc)
         read(28,15,end=17) ((zc(i,j,k),j=1,nyc),i=1,nxc)
       enddo
 15   format(5(1pe13.6))
C
C  now average the adjacent values to give 'cell centered' values
C  compress by 1 in each dimension in turn
C
c 17   if(interp .eq. 1) then
 17    do k = 1,nzl
        do j = 1,nyl
         do i = 1,nxl
          coords(i,j,k,1) = xc(i,j,k)
          coords(i,j,k,2) = yc(i,j,k)
          coords(i,j,k,3) = zc(i,j,k)
         enddo
        enddo
       enddo 
c      else
c       do k = 1,nz
c        do j = 1,ny
c         do i = 1,nx
c          coords(i,j,k,1) = .5*(xc(i,j,k)+xc(i+1,j,k))
c          coords(i,j,k,2) = .5*(yc(i,j,k)+yc(i,j+1,k))
c          coords(i,j,k,3) = .5*(zc(i,j,k)+zc(i,j,k+1))
c         enddo
c        enddo
c       enddo
c      endif
C
C   set return flag successful
C
      phoenics_bfc_read = 1
      return
      end
C**********************************************************************
      integer function phoenics_bfc_reorder(field,nxl,nyl,nzl,nv,
     *                            lu,lv,lw,name_stored,nphi)
c
C  this function reorders the data array to reflect the fact that AVS
C  likes the u, v, and w velocity components to be stored in that order
C
      integer nxl,nyl,nzl,nv,lu,lv,lw,low,mid,high,nphi,i,j,k
      dimension field(nv,nxl,nyl,nzl)
      character*4 name_stored(nphi),tname1,tname2,tname3
      real temp1,temp2,temp3
C
C  reorder the data array so that the lu values are first, lv are second
C  and lw are third
C
      print*, ' in function bfc reorder'
c     print*, nxl,nyl,nzl
      low = min(lu,lv,lw)
      high = max(lu,lv,lw)
      mid = lu + lv + lw - low - high
c     if(lu .eq. 0 .and. lv .eq. 0) return
c     if(lv .eq. 0 .and. lw .eq. 0) return
c     if(lu .eq. 0 .and. lw .eq. 0) return
      do k=1,nzl
        do j = 1,nyl
          do i = 1,nxl
c           if(lu .eq. 0) then
c            temp2 = field(lv,i,j,k)
c            temp3 = field(lw,i,j,k)
c            field(mid,i,j,k) = temp2
c            field(high,i,j,k) = temp3
c           else
c            if(lv .eq. 0) then
c             temp1 = field(lu,i,j,k)
c             temp3 = field(lw,i,j,k)
c             field(mid,i,j,k) = temp1
c             field(high,i,j,k) = temp3
c            else
c             if(lw .eq. 0) then
c              temp1 = field(lu,i,j,k)
c              temp2 = field(lv,i,j,k)
c              field(mid,i,j,k) = temp1
c              field(high,i,j,k) = temp2
c             else
               temp1 = field(lu,i,j,k)
               temp2 = field(lv,i,j,k)
               temp3 = field(lw,i,j,k)
               field(low,i,j,k) = temp1
               field(mid,i,j,k) = temp2
               field(high,i,j,k) = temp3
c             endif
c            endif
c           endif
          enddo
        enddo
      enddo
C
C  reorder the name_stored array
C
c     if(lu .eq. 0) then
c      tname2 = name_stored(lv)
c      tname3 = name_stored(lw)
c      name_stored(mid) = tname2
c      name_stored(high) = tname3
c     else
c      if(lv .eq. 0) then
c       tname1 = name_stored(lu)
c       tname3 = name_stored(lw)
c       name_stored(mid) = tname1
c       name_stored(high) = tname3
c      else
c       if(lw .eq. 0) then
c        tname1 = name_stored(lu)
c        tname2 = name_stored(lv)
c        name_stored(mid) = tname1
c        name_stored(high) = tname2
c       else
         tname1 = name_stored(lu)
         tname2 = name_stored(lv)
         tname3 = name_stored(lw)
         name_stored(low) = tname1
         name_stored(mid) = tname2
         name_stored(high)= tname3
c       endif
c      endif
c     endif
c
C  signify successful completion
C
      phoenics_bfc_reorder = 1
      return
      end
C*******************************************************************
C
      integer function phoenics_cart_areas(xc,yc,zc,nx,ny,nz,ax,ay,az)
C
C   this function computes the facial areas ax,ay and az from the 
C   cell corner locations stored in xc, yc and zc
C
      real xc(0:nx), yc(0:ny), zc(0:nz)
      real ax(0:nx,ny,nz), ay(nx,0:ny,nz), az(nx,ny,0:nz)
C
C   compute the x face areas
C
      print*,' in cartes areas function'
      print*, nx,ny,nz
      print*, (xc(i),i=0,nx)
      print*, (yc(j),j=0,ny)
      print*, (zc(k),k=0,nz)
      do i = 0,nx
       do k = 1,nz
        do j = 1,ny
         ax(i,j,k) = (yc(j) - yc(j-1))*(zc(k) - zc(k-1))
         ax(i,j,k) = max(ax(i,j,k),1.e-10)
        enddo
       enddo
       print*,i,((ax(i,jj,kk),jj=1,ny,2),kk=1,nz,2)
      enddo
      do j = 0,ny
       do k = 1,nz
        do i = 1,nx
         ay(i,j,k) = (zc(k) - zc(k-1))*(xc(i) - xc(i-1))
         ay(i,j,k) = max(ay(i,j,k),1.e-10)
        enddo
       enddo
      enddo
      do k = 0,nz
       do j = 1,ny
        do i = 1,nx
         az(i,j,k) = (xc(i) - xc(i-1))*(yc(j) - yc(j-1))
         az(i,j,k) = max(az(i,j,k),1.e-10)
        enddo
       enddo
      enddo
C
      phoenics_cart_areas = 1
      return
      end
C*************************************************************************
C
      integer function phoenics_cart_vols(xc,yc,zc,nx,ny,nz,vols)
C
C   this function computes the cell volumes for cartesian cases, based
C   on the cell corner locations in the xc, yc and zc arrays
C
      real xc(0:nx),yc(0:ny),zc(0:nz)
      real vols(nx,ny,nz)
C
      print*,' in cartes vols function'
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         vols(i,j,k) = (xc(i) - xc(i-1))*
     1                 (yc(j) - yc(j-1))*
     2                 (zc(k) - zc(k-1))
         vols(i,j,k) = max(vols(i,j,k),1.e-10)
        enddo
       enddo
      enddo
C
      phoenics_cart_vols = 1
      return
      end
C***********************************************************************
C
      integer function phoenics_cyl_areas(xc,yc,zc,nx,ny,nz,ax,ay,az,
     1                                    rinner)
C
C   this function computes the facial areas ax, ay and az for the
C   cylindrical case with xc the theta, yc the radial and zc the axial
C   cell face locations
C
      real xc(0:nx), yc(0:ny), zc(0:nz)
      real ax(0:nx,ny,nz), ay(nx,0:ny,nz), az(nx,ny,0:nz)
C
      print*,' in cyl areas function'
      do i = 0,nx
       do k = 1,nz
        do j = 1,ny
         ax(i,j,k) = max(1.e-10,(zc(k) - zc(k-1))*(yc(j) - yc(j-1)))
        enddo
       enddo
      enddo
      do j = 0,ny
       do k = 1,nz
        do i = 1,nx
         ay(i,j,k) = max(1.e-10,(rinner + yc(j))*(xc(i) - xc(i-1))*
     1                                (zc(k) - zc(k-1)))
        enddo
       enddo
      enddo
      do k = 0,nz
       do j = 1,ny
        do i = 1,nx
         az(i,j,k) = max(1.e-10,(rinner + .5*(yc(j) + yc(j-1)))*
     1                                       (xc(i) - xc(i-1))*
     1                                       (yc(j) - yc(j-1)))
        enddo
       enddo
      enddo
C
      phoenics_cyl_areas = 1
      return
      end
C**********************************************************************
C
      integer function phoenics_cyl_vols(xc,yc,zc,nx,ny,nz,vols,rinner)
C
C  this function computes the cell volumes for the cylindrical case using
C  the cell corner information stored in the xc, yc, and zc arrays
C
      real xc(0:nx), yc(0:ny), zc(0:nz)
      real vols(nx,ny,nz)
C
      print*, ' in cyl vols function'
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         vols(i,j,k) = (rinner + .5*(yc(j)+yc(j-1)))*
     1                              (xc(i) - xc(i-1))*
     1                              (yc(j) - yc(j-1))*
     2                              (zc(k) - zc(k-1))
         if(vols(i,j,k) .lt. 1.e-10) print*, i,j,k
        enddo
       enddo
      enddo
C
      phoenics_cyl_vols = 1
      return
      end
C**************************************************************************
C
      integer function phoenics_bfc_areas(xyzname,nx,ny,nz,ax,ay,az)
C
C   this function computes the cell facial areas for the bfc case 
C   using the corner location information in the xc, yc and zc arrays
C   the corners are assumed to be coplanar
C
      parameter (nx1 = 100, ny1 = 100, nz1 = 100)
      character*(*) xyzname
      dimension xc(0:nx1,0:ny1,0:nz1),yc(0:nx1,0:ny1,0:nz1),
     *    zc(0:nx1,0:ny1,0:nz1)
      dimension ax(0:nx,ny,nz), ay(nx,0:ny,nz), az(nx,ny,0:nz)
C
C  error exit if file name not specified
C
      print*,' in bfc areas function'
      if(xyzname .eq. '') then
          phoenics_bfc_areas = 0
          return
      endif
C
C  open the file
C  
      open(unit=28,file=xyzname,status='OLD')
      rewind(28)
      read(28,*) nxc,nyc,nzc
C
C  now read in the values 
C
      do k = 0,nz
         read(28,15,end=17) ((xc(i,j,k),j=0,ny),i=0,nx)
         read(28,15,end=17) ((yc(i,j,k),j=0,ny),i=0,nx)
         read(28,15,end=17) ((zc(i,j,k),j=0,ny),i=0,nx)
       enddo
 15   format(5(1pe13.6))
C
 17   do i = 0,nx
       do k = 1,nz
        do j = 1,ny
         p = sqrt((xc(i,j,k-1) - xc(i,j-1,k))**2
     1           +(yc(i,j,k-1) - yc(i,j-1,k))**2
     2           +(zc(i,j,k-1) - zc(i,j-1,k))**2)
         q = sqrt((xc(i,j,k) - xc(i,j-1,k-1))**2
     1           +(yc(i,j,k) - yc(i,j-1,k-1))**2
     2           +(zc(i,j,k) - zc(i,j-1,k-1))**2)
         a = sqrt((xc(i,j,k-1) - xc(i,j,k))**2
     1           +(yc(i,j,k-1) - yc(i,j,k))**2
     2           +(zc(i,j,k-1) - zc(i,j,k))**2)
         b = sqrt((xc(i,j,k-1) - xc(i,j-1,k-1))**2
     1           +(yc(i,j,k-1) - yc(i,j-1,k-1))**2
     2           +(zc(i,j,k-1) - zc(i,j-1,k-1))**2)
         c = sqrt((xc(i,j-1,k) - xc(i,j-1,k-1))**2
     1           +(yc(i,j-1,k) - yc(i,j-1,k-1))**2
     2           +(zc(i,j-1,k) - zc(i,j-1,k-1))**2)
         d = sqrt((xc(i,j,k) - xc(i,j-1,k))**2
     1           +(yc(i,j,k) - yc(i,j-1,k))**2
     2           +(zc(i,j,k) - zc(i,j-1,k))**2)
         ax(i,j,k) = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
        enddo
       enddo
      enddo
C      
      do j = 0,ny
       do k = 1,nz
        do i = 1,nx
         p = sqrt((xc(i,j,k-1) - xc(i-1,j,k))**2
     1           +(yc(i,j,k-1) - yc(i-1,j,k))**2
     2           +(zc(i,j,k-1) - zc(i-1,j,k))**2)
         q = sqrt((xc(i,j,k) - xc(i-1,j,k-1))**2
     1           +(yc(i,j,k) - yc(i-1,j,k-1))**2
     2           +(zc(i,j,k) - zc(i-1,j,k-1))**2)
         a = sqrt((xc(i,j,k-1) - xc(i,j,k))**2
     1           +(yc(i,j,k-1) - yc(i,j,k))**2
     2           +(zc(i,j,k-1) - zc(i,j,k))**2)
         b = sqrt((xc(i,j,k-1) - xc(i-1,j,k-1))**2
     1           +(yc(i,j,k-1) - yc(i-1,j,k-1))**2
     2           +(zc(i,j,k-1) - zc(i-1,j,k-1))**2)
         c = sqrt((xc(i-1,j,k) - xc(i-1,j,k-1))**2
     1           +(yc(i-1,j,k) - yc(i-1,j,k-1))**2
     2           +(zc(i-1,j,k) - zc(i-1,j,k-1))**2)
         d = sqrt((xc(i,j,k) - xc(i-1,j,k))**2
     1           +(yc(i,j,k) - yc(i-1,j,k))**2
     2           +(zc(i,j,k) - zc(i-1,j,k))**2)
         ay(i,j,k) = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
        enddo
       enddo
      enddo
C      
      do k = 0,nz
       do j = 1,ny
        do i = 1,nx
         p = sqrt((xc(i,j-1,k) - xc(i-1,j,k))**2
     1           +(yc(i,j-1,k) - yc(i-1,j,k))**2
     2           +(zc(i,j-1,k) - zc(i-1,j,k))**2)
         q = sqrt((xc(i,j,k) - xc(i-1,j-1,k))**2
     1           +(yc(i,j,k) - yc(i-1,j-1,k))**2
     2           +(zc(i,j,k) - zc(i-1,j-1,k))**2)
         a = sqrt((xc(i,j-1,k) - xc(i,j,k))**2
     1           +(yc(i,j-1,k) - yc(i,j,k))**2
     2           +(zc(i,j-1,k) - zc(i,j,k))**2)
         b = sqrt((xc(i,j-1,k) - xc(i-1,j-1,k))**2
     1           +(yc(i,j-1,k) - yc(i-1,j-1,k))**2
     2           +(zc(i,j-1,k) - zc(i-1,j-1,k))**2)
         c = sqrt((xc(i-1,j,k) - xc(i-1,j-1,k))**2
     1           +(yc(i-1,j,k) - yc(i-1,j-1,k))**2
     2           +(zc(i-1,j,k) - zc(i-1,j-1,k))**2)
         d = sqrt((xc(i,j,k) - xc(i-1,j,k))**2
     1           +(yc(i,j,k) - yc(i-1,j,k))**2
     2           +(zc(i,j,k) - zc(i-1,j,k))**2)
         az(i,j,k) = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
        enddo
       enddo
      enddo
C
      phoenics_bfc_areas = 1
      return
      end
C**********************************************************************
C
      integer function phoenics_bfc_vols(xyzname,nx,ny,nz,vols)
C 
C   this function computes the cell volumes for the bfc case based
C   on the cell corner location information stored in the xc, yc and
C   zc arrays.  the cell volume is computed approximately from the
C   average of the volumes calculated using the facial areas in each
C   direction multiplied by the average height in the direction 
C   under consideration.  this should give results close to the 
C   actual volumes
C
      parameter (nx1 = 100, ny1 = 100, nz1 = 100)
      character*(*) xyzname
      dimension xc(0:nx1,0:ny1,0:nz1),yc(0:nx1,0:ny1,0:nz1),
     *    zc(0:nx1,0:ny1,0:nz1)
      dimension vols(nx,ny,nz)
C
C  error exit if file name not specified
C
      print*,' in bfc volumes function'
      if(xyzname .eq. '') then
          phoenics_vols = 0
          return
      endif
C
C  open the file
C  
      open(unit=28,file=xyzname,status='OLD')
      rewind(28)
      read(28,*) nxc,nyc,nzc
C
C  now read in the values 
C
      do k = 0,nz
         read(28,15,end=17) ((xc(i,j,k),j=0,ny),i=0,nx)
         read(28,15,end=17) ((yc(i,j,k),j=0,ny),i=0,nx)
         read(28,15,end=17) ((zc(i,j,k),j=0,ny),i=0,nx)
       enddo
 15   format(5(1pe13.6))
C
C
 17   do i = 1,nx
       do k = 1,nz
        do j = 1,ny
         p = sqrt((xc(i,j,k-1) - xc(i,j-1,k))**2
     1           +(yc(i,j,k-1) - yc(i,j-1,k))**2
     2           +(zc(i,j,k-1) - zc(i,j-1,k))**2)
         q = sqrt((xc(i,j,k) - xc(i,j-1,k-1))**2
     1           +(yc(i,j,k) - yc(i,j-1,k-1))**2
     2           +(zc(i,j,k) - zc(i,j-1,k-1))**2)
         a = sqrt((xc(i,j,k-1) - xc(i,j,k))**2
     1           +(yc(i,j,k-1) - yc(i,j,k))**2
     2           +(zc(i,j,k-1) - zc(i,j,k))**2)
         b = sqrt((xc(i,j,k-1) - xc(i,j-1,k-1))**2
     1           +(yc(i,j,k-1) - yc(i,j-1,k-1))**2
     2           +(zc(i,j,k-1) - zc(i,j-1,k-1))**2)
         c = sqrt((xc(i,j-1,k) - xc(i,j-1,k-1))**2
     1           +(yc(i,j-1,k) - yc(i,j-1,k-1))**2
     2           +(zc(i,j-1,k) - zc(i,j-1,k-1))**2)
         d = sqrt((xc(i,j,k) - xc(i,j-1,k))**2
     1           +(yc(i,j,k) - yc(i,j-1,k))**2
     2           +(zc(i,j,k) - zc(i,j-1,k))**2)
         ax = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
         pdis = .25*((xc(i,j,k) - xc(i-1,j,k))+
     1               (xc(i,j-1,k) - xc(i-1,j-1,k))+
     2               (xc(i,j,k-1) - xc(i-1,j,k-1))+
     3               (xc(i,j-1,k-1) - xc(i-1,j-1,k-1)))
         vols(i,j,k) = ax * pdis
        enddo
       enddo
      enddo
C      
      do j = 1,ny
       do k = 1,nz
        do i = 1,nx
         p = sqrt((xc(i,j,k-1) - xc(i-1,j,k))**2
     1           +(yc(i,j,k-1) - yc(i-1,j,k))**2
     2           +(zc(i,j,k-1) - zc(i-1,j,k))**2)
         q = sqrt((xc(i,j,k) - xc(i-1,j,k-1))**2
     1           +(yc(i,j,k) - yc(i-1,j,k-1))**2
     2           +(zc(i,j,k) - zc(i-1,j,k-1))**2)
         a = sqrt((xc(i,j,k-1) - xc(i,j,k))**2
     1           +(yc(i,j,k-1) - yc(i,j,k))**2
     2           +(zc(i,j,k-1) - zc(i,j,k))**2)
         b = sqrt((xc(i,j,k-1) - xc(i-1,j,k-1))**2
     1           +(yc(i,j,k-1) - yc(i-1,j,k-1))**2
     2           +(zc(i,j,k-1) - zc(i-1,j,k-1))**2)
         c = sqrt((xc(i-1,j,k) - xc(i-1,j,k-1))**2
     1           +(yc(i-1,j,k) - yc(i-1,j,k-1))**2
     2           +(zc(i-1,j,k) - zc(i-1,j,k-1))**2)
         d = sqrt((xc(i,j,k) - xc(i-1,j,k))**2
     1           +(yc(i,j,k) - yc(i-1,j,k))**2
     2           +(zc(i,j,k) - zc(i-1,j,k))**2)
         ay = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
         pdis = .25*((yc(i,j,k) - yc(i,j-1,k))+
     1               (yc(i-1,j,k) - yc(i-1,j-1,k))+
     2               (yc(i,j,k-1) - yc(i,j-1,k-1))+
     3               (yc(i-1,j,k-1) - yc(i-1,j-1,k-1)))
         vols(i,j,k) = vols(i,j,k) + ay * pdis
        enddo
       enddo
      enddo
C      
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         p = sqrt((xc(i,j-1,k) - xc(i-1,j,k))**2
     1           +(yc(i,j-1,k) - yc(i-1,j,k))**2
     2           +(zc(i,j-1,k) - zc(i-1,j,k))**2)
         q = sqrt((xc(i,j,k) - xc(i-1,j-1,k))**2
     1           +(yc(i,j,k) - yc(i-1,j-1,k))**2
     2           +(zc(i,j,k) - zc(i-1,j-1,k))**2)
         a = sqrt((xc(i,j-1,k) - xc(i,j,k))**2
     1           +(yc(i,j-1,k) - yc(i,j,k))**2
     2           +(zc(i,j-1,k) - zc(i,j,k))**2)
         b = sqrt((xc(i,j-1,k) - xc(i-1,j-1,k))**2
     1           +(yc(i,j-1,k) - yc(i-1,j-1,k))**2
     2           +(zc(i,j-1,k) - zc(i-1,j-1,k))**2)
         c = sqrt((xc(i-1,j,k) - xc(i-1,j-1,k))**2
     1           +(yc(i-1,j,k) - yc(i-1,j-1,k))**2
     2           +(zc(i-1,j,k) - zc(i-1,j-1,k))**2)
         d = sqrt((xc(i,j,k) - xc(i-1,j,k))**2
     1           +(yc(i,j,k) - yc(i-1,j,k))**2
     2           +(zc(i,j,k) - zc(i-1,j,k))**2)
         az = .25*sqrt(4.0*p*p*q*q - (b*b+d*d - a*a- c*c)**2)
         pdis = .25*((yc(i,j,k) - yc(i,j,k-1))+
     1               (yc(i-1,j,k) - yc(i-1,j,k-1))+
     2               (yc(i,j-1,k) - yc(i,j-1,k-1))+
     3               (yc(i-1,j-1,k) - yc(i-1,j-1,k-1)))
         vols(i,j,k) = (vols(i,j,k) + az * pdis)/3.0
        enddo
       enddo
      enddo
C
      phoenics_bfc_vols = 1
      return
      end
C***************************************************************************
C
      integer function uvel_extrap(field,iu,nxl,nyl,nzl,nv,
     *                                        coords,ncoord)
C
C   this function extrapolates the u velocty components to east-most
C   i = 1 face of the domain and to the west most i=nx face of the domain
C   by simple linear extrapolation
C
      dimension field(nv,nxl,nyl,nzl),coords(nxl,nyl,nzl,ncoord)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*, ' in uvel_extrap'
      if(nxl .gt. 2) then
       do k = 2,nzl
        do j = 2,nyl
          dx21 = .25 * (coords(2,j,k,1) - coords(1,j,k,1) +
     1                coords(2,j-1,k,1) - coords(1,j-1,k,1) +
     2                coords(2,j,k-1,1) - coords(1,j,k-1,1) +
     3                coords(2,j-1,k-1,1) - coords(1,j-1,k-1,1))
          dx32 = .25 * (coords(3,j,k,1) - coords(2,j,k,1) +
     1                coords(3,j-1,k,1) - coords(2,j-1,k,1) +
     2                coords(3,j,k-1,1) - coords(2,j,k-1,1) +
     3                coords(3,j-1,k-1,1) - coords(2,j-1,k-1,1))
          dxr = dx21 / dx32
          field(iu,1,j,k) = field(iu,2,j,k) + dxr*field(iu,2,j,k) 
     1                                    - dxr*field(iu,3,j,k)
         enddo
        enddo
        do k = 2,nzl
         do j = 2,nyl
          dxnx1 = .25 *(coords(nxl,j,k,1) - coords(nxl-1,j,k,1) +
     1               coords(nxl,j-1,k,1) - coords(nxl-1,j-1,k,1) +
     2               coords(nxl,j,k-1,1) - coords(nxl-1,j,k-1,1) +
     3               coords(nxl,j-1,k-1,1) - coords(nxl-1,j-1,k-1,1))
          dxnx2 = .25 *(coords(nxl-1,j,k,1) - coords(nxl-2,j,k,1) +
     1               coords(nxl-1,j-1,k,1) - coords(nxl-2,j-1,k,1) +
     2               coords(nxl-1,j,k-1,1) - coords(nxl-2,j,k-1,1) +
     3               coords(nxl-1,j-1,k-1,1) - coords(nxl-2,j-1,k-1,1))
          dxr = dxnx1 / dxnx2
          field(iu,nxl,j,k) = field(iu,nxl-1,j,k) 
     1                 +dxr*field(iu,nxl-1,j,k)-dxr*field(iu,nxl-2,j,k)
         enddo
       enddo
      else
       do k = 2,nzl
         do j = 2,nyl
           field(iu,1,j,k) = field(iu,2,j,k)
           field(iu,nxl,j,k) = field(iu,nxl-1,j,k)
         enddo
       enddo
      endif
C
      uvel_extrap = 1
      return
      end 
C***************************************************************************
C
      integer function vvel_extrap(field,iv,nxl,nyl,nzl,nv,
     *                                        coords,ncoord)
C
C   this function extrapolates the v velocity components to south-most
C   j = 1 face of the domain and to the north most j = nyl face of the domain
C   by simple linear extrapolation
C
      dimension field(nv,nxl,nyl,nzl),coords(nxl,nyl,nzl,ncoord)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*,' in vvel extrap'
      if(nyl .gt. 2) then
       do k = 2,nzl
        do i = 2,nxl
          dy21 = .25 * (coords(i,2,k,2) - coords(i,1,k,2) +
     1                  coords(i-1,2,k,2) - coords(i-1,1,k,2) +
     2                  coords(i,2,k-1,2) - coords(i,1,k-1,2) +
     3                  coords(i-1,2,k-1,2) - coords(i-1,1,k-1,2))
          dy32 = .25 * (coords(i,3,k,2) - coords(i,2,k,2) +
     1                  coords(i-1,3,k,2) - coords(i-1,2,k,2) +
     2                  coords(i,3,k-1,2) - coords(i,2,k-1,2) +
     3                  coords(i-1,3,k-1,2) - coords(i-1,2,k-1,2))
          dyr = dy21 / dy32
          field(iv,i,1,k) = field(iv,i,2,k) + dyr*field(iv,i,2,k) 
     1                                    - dyr*field(iv,i,3,k)
         enddo
        enddo
        do k = 2,nzl
         do j = 2,nyl
          dynyl1 = .25 *(coords(i,nyl,k,2) - coords(i,nyl-1,k,2) +
     1              coords(i-1,nyl,k,2) - coords(i-1,nyl-1,k,2) +
     2              coords(i,nyl,k-1,2) - coords(i,nyl-1,k-1,2) +
     3              coords(i-1,nyl,k-1,2) - coords(i-1,nyl-1,k-1,2))
          dynyl2 = .25 *(coords(i,nyl-1,k,2) - coords(i,nyl-2,k,2) +
     1              coords(i-1,nyl-1,k,2) - coords(i-1,nyl-2,k,2) +
     2              coords(i,nyl-1,k-1,2) - coords(i,nyl-2,k-1,2) +
     3              coords(i-1,nyl-1,k-1,2) - coords(i-1,nyl-1,k-1,2))
          dyr = dynyl1 / dynyl2
          field(iv,i,nyl,k) = field(iv,i,nyl-1,k)
     1                 +dyr*field(iv,i,nyl-1,k)-dyr*field(iv,i,nyl-2,k)
         enddo
       enddo
      else
       do k = 2,nzl
         do i = 2,nxl
           field(iv,i,1,k) = field(iv,i,2,k)
           field(iv,i,nyl,k) = field(iv,i,nyl-1,k)
         enddo
       enddo
      endif
C
      vvel_extrap = 1
      return
      end
C***************************************************************************
C
      integer function wvel_extrap(field,iw,nxl,nyl,nzl,nv,
     *                                        coords,ncoord)
C
C   this function extrapolates the w velocty components to low-most
C   k = 1 face of the domain and to the high most k=nzl face of the domain
C   by simple linear extrapolation
C
      dimension field(nv,nxl,nyl,nzl),coords(nxl,nyl,nzl,ncoord)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*,' in wvel extrap'
      if(nzl .gt. 2) then
       do j = 2,nyl
        do i = 2,nxl
          dz21 = .25 * (coords(i,j,2,3) - coords(i,j,1,3) +
     1                  coords(i-1,j,2,3) - coords(i-1,j,1,3) +
     2                  coords(i,j-1,2,3) - coords(i,j-1,1,3) +
     3                  coords(i-1,j-1,2,3) - coords(i-1,j-1,1,3))
          dz32 = .25 * (coords(i,j,3,3) - coords(i,j,2,3) +
     1                  coords(i-1,j,3,3) - coords(i-1,j,2,3) +
     2                  coords(i,j-1,3,3) - coords(i,j-1,2,3) +
     3                  coords(i-1,j-1,3,3) - coords(i-1,j-1,2,3))
          dzr = dz21 / dz32
          field(iw,i,j,1) = field(iw,i,j,2) + dzr*field(iw,i,j,2) 
     1                                    - dzr*field(iw,i,j,3)
         enddo
        enddo
        do j = 2,nyl
         do i = 2,nxl
          dznzl1 = .25 *(coords(i,j,nzl,3) - coords(i-i,j,nzl-1,3) +
     1              coords(i-1,j,nzl,3) - coords(i-1,j,nzl-1,3) +
     2              coords(i,j,k-1,3) - coords(i-i,j,nzl-1,3) +
     3              coords(i,j-1,k-1,3) - coords(i-1,j,nzl-1,3))
          dznzl2 = .25 *(coords(i,j,nzl-1,3) - coords(i,j,nzl-2,3) +
     1              coords(i-1,j,nzl-1,3) - coords(i-1,j,nzl-2,3) +
     2              coords(i,j-1,nzl-1,3) - coords(i,j-1,nzl-2,3) +
     3              coords(i-1,j-1,nzl-1,3) - coords(i-1,j,nzl-2,3))
          dzr = dznzl1 / dznzl2
          field(iw,i,j,nzl) = field(iw,i,j,nzl-1)
     1              + dzr*field(iw,i,j,nzl-1)-dzr*field(iw,i,j,nzl-2)
         enddo
       enddo
      else
       do j = 2,nyl
         do i = 2,nxl
           field(iw,i,j,1) = field(iw,i,j,2)
           field(iw,i,j,nzl) = field(iw,i,j,nzl-1)
         enddo
       enddo
      endif
C
      wvel_extrap = 1
      return
      end
C***********************************************************************
C
      integer function uvel_interp(field,iv,nxl,nyl,nzl,nv,ax,nx,ny,nz)
C
C   this function interpolates the u component velocity data in the
C   data array into nodal point located information, based on the
C   x face cell areas stored in the ax array.  This is an adaptation of
C   the interpolation routine found in the mpgs conversion program
C
      parameter (nxt = 101, nyt = 101, nzt = 101)
      integer nxl,nyl,nzl,nx,ny,nz,nv,iv
      dimension field(nv,nxl,nyl,nzl), ax(0:nx,ny,nz)
      dimension vfunc(nxt*nyt*nzt),sfunc(nxt*nyt*nzt)
C
C  first zero out the accumulation arrays and initialize the vfunc array
C
      print*, ' in uvel interp'
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         vfunc(nn) = field(iv,i,j,k)
         field(iv,i,j,k) = 0.0
         sfunc(nn) = 0.0
        enddo
       enddo
      enddo
C
C  now step through the domain, using surface area X velocity products
C  to average the velocity to each nodal point
C
      nn = 0
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) + vfunc(nn+nxl+nxl*nyl)*
     *                                               ax(i-1,j,k)
         sfunc(nn) = sfunc(nn) + ax(i-1,j,k)
C
         field(iv,i+1,j,k) = field(iv,i+1,j,k)+vfunc(nn+1+nxl+nxl*nyl)*
     *                                               ax(i,j,k)
         sfunc(nn+1) = sfunc(nn+1) + ax(i,j,k)
C
         field(iv,i,j+1,k) = field(iv,i,j+1,k) + vfunc(nn+nxl+nxl*nyl)*
     *                                               ax(i-1,j,k)
         sfunc(nn+nxl) = sfunc(nn+nxl) + ax(i-1,j,k)
C
         field(iv,i,j,k+1) = field(iv,i,j,k+1) + vfunc(nn+nxl+nxl*nyl)*
     *                                               ax(i-1,j,k)
         sfunc(nn+nxl*nyl) = sfunc(nn+nxl*nyl) + ax(i-1,j,k)
C
         field(iv,i+1,j+1,k) = field(iv,i+1,j+1,k)
     *                           +vfunc(nn+1+nxl+nxl*nyl)*ax(i,j,k)
         sfunc(nn+1+nxl) = sfunc(nn+1+nxl) + ax(i,j,k)
C
         field(iv,i,j+1,k+1) = field(iv,i,j+1,k+1)
     *                           +vfunc(nn+nxl+nxl*nyl)*ax(i-1,j,k)
         sfunc(nn+nxl+nxl*nyl) = sfunc(nn+nxl+nxl*nyl) + ax(i-1,j,k)
C
         field(iv,i+1,j,k+1) = field(iv,i+1,j,k+1)
     *                           +vfunc(nn+1+nxl+nxl*nyl)*ax(i,j,k)
         sfunc(nn+1+nxl*nyl) = sfunc(nn+1+nxl*nyl) + ax(i,j,k)
C
         field(iv,i+1,j+1,k+1) = field(iv,i+1,j+1,k+1)
     *                           +vfunc(nn+1+nxl+nxl*nyl)*ax(i,j,k)
         sfunc(nn+1+nxl+nxl*nyl) = sfunc(nn+1+nxl+nxl*nyl) + ax(i,j,k)
        enddo
        nn = nn + 1
       enddo
       nn = nn + nxl
      enddo
C
C  now divide the sums by the cell face areas
C
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) / sfunc(nn)
        enddo
       enddo
      enddo
C
C  return on successful completion
C
      uvel_interp = 1
      return
      end
 
C***********************************************************************
C
      integer function vvel_interp(field,iv,nxl,nyl,nzl,nv,ay,nx,ny,nz)
C
C   this function interpolates the v component velocity data in the
C   data array into nodal point located information, based on the
C   x face cell areas stored in the ay array.  This is an adaptation of
C   the interpolation routine found in the mpgs conversion program
C
      parameter (nxt = 101, nyt = 101, nzt = 101)
      integer nxl,nyl,nzl,nx,ny,nz,nv,iv
      dimension field(nv,nxl,nyl,nzl), ay(nx,0:ny,nz)
      dimension vfunc(nxt*nyt*nzt),sfunc(nxt*nyt*nzt)
C
C  first zero out the accumulation arrays and initialize the vfunc array
C
      print*,' in vvel interp'
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         vfunc(nn) = field(iv,i,j,k)
         field(iv,i,j,k) = 0.0
         sfunc(nn) = 0.0
        enddo
       enddo
      enddo
C
C  now step through the domain, using surface area X velocity products
C  to average the velocity to each nodal point
C
      nn = 0
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) + vfunc(nn+1+nxl*nyl)*
     *                                               ay(i,j-1,k)
         sfunc(nn) = sfunc(nn) + ay(i,j-1,k)
C
         field(iv,i+1,j,k) = field(iv,i+1,j,k) + 
     *                    vfunc(nn+1+nxl*nyl)*ay(i,j-1,k)
         sfunc(nn+1) = sfunc(nn+1) + ay(i,j-1,k)
C
         field(iv,i,j+1,k) = field(iv,i,j+1,k) + 
     *                    vfunc(nn+1+nxl+nxl*nyl)*ay(i,j,k)
         sfunc(nn+nxl) = sfunc(nn+nxl) + ay(i,j,k)
C
         field(iv,i,j,k+1) = field(iv,i,j,k+1) + 
     *                    vfunc(nn+1+nxl*nyl)*ay(i,j-1,k)
         sfunc(nn+nxl*nyl) = sfunc(nn+nxl*nyl) + ay(i,j-1,k)
C
         field(iv,i+1,j+1,k) = field(iv,i+1,j+1,k) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*ay(i,j,k)
         sfunc(nn+1+nxl) = sfunc(nn+1+nxl) + ay(i,j,k)
C
         field(iv,i+1,j,k+1) = field(iv,i+1,j,k+1) + 
     *                          vfunc(nn+1+nxl*nyl)*ay(i,j-1,k)
         sfunc(nn+1+nxl*nyl) = sfunc(nn+1+nxl*nyl) + ay(i,j-1,k)
C
         field(iv,i,j+1,k+1) = field(iv,i,j+1,k+1) + 
     *                          vfunc(nn+nxl+nxl*nyl)*ay(i,j,k)
         sfunc(nn+nxl+nxl*nyl) = sfunc(nn+nxl+nxl*nyl) + ay(i,j,k)
C
         field(iv,i+1,j+1,k+1) = field(iv,i+1,j+1,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*ay(i,j,k)
         sfunc(nn+1+nxl+nxl*nyl) = sfunc(nn+1+nxl+nxl*nyl)+ay(i,j,k)
        enddo
        nn = nn + 1
       enddo
       nn = nn + nxl
      enddo
C
C  now divide the sums by the cell face areas
C
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) / sfunc(nn)
        enddo
       enddo
      enddo
C
C  return on successful completion
C
      print*,'locating ',field(3,1,nyl,nzl)
      vvel_interp = 1
      return
      end
C***********************************************************************
C
      integer function wvel_interp(field,iv,nxl,nyl,nzl,nv,az,nx,ny,nz)
C
C   this function interpolates the u component velocity data in the
C   data array into nodal point located information, based on the
C   x face cell areas stored in the ax array.  This is an adaptation of
C   the interpolation routine found in the mpgs conversion program
C
      parameter (nxt = 101, nyt = 101, nzt = 101)
      integer nxl,nyl,nzl,nx,ny,nz,nv,iv
      dimension field(nv,nxl,nyl,nzl), az(nx,ny,0:nz)
      dimension vfunc(nxt*nyt*nzt),sfunc(nxt*nyt*nzt)
C
C  first zero out the accumulation arrays and initialize the vfunc array
C
      print*,' in wvel interp'
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         vfunc(nn) = field(iv,i,j,k)
         field(iv,i,j,k) = 0.0
         sfunc(nn) = 0.0
        enddo
       enddo
      enddo
C
C  now step through the domain, using surface area X velocity products
C  to average the velocity to each nodal point
C
      nn = 0
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) + vfunc(nn+1+nxl)*
     *                                               az(i,j,k-1)
         sfunc(nn) = sfunc(nn) + az(i,j,k-1)
C
         field(iv,i+1,j,k) = field(iv,i+1,j,k) + vfunc(nn+1+nxl)*
     *                                               az(i,j,k-1)
         sfunc(nn+1) = sfunc(nn+1) + az(i,j,k-1)
C
         field(iv,i,j+1,k) = field(iv,i,j+1,k) + vfunc(nn+1+nxl)*
     *                                               az(i,j,k-1)
         sfunc(nn+nxl) = sfunc(nn+nxl) + az(i,j,k-1)
C
         field(iv,i,j,k+1) = field(iv,i,j,k+1)+vfunc(nn+1+nxl+nxl*nyl)*
     *                                               az(i,j,k)
         sfunc(nn+nxl*nyl) = sfunc(nn+nxl*nyl) + az(i,j,k)
C
         field(iv,i+1,j+1,k) = field(iv,i+1,j+1,k) + 
     *                          vfunc(nn+1+nxl)*az(i,j,k-1)
         sfunc(nn+1+nxl) = sfunc(nn+1+nxl) + az(i,j,k-1)
C
         field(iv,i+1,j,k+1) = field(iv,i+1,j,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*az(i,j,k)
         sfunc(nn+1+nxl*nyl) = sfunc(nn+1+nxl*nyl) + az(i,j,k)
C
         field(iv,i,j+1,k+1) = field(iv,i,j+1,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*az(i,j,k)
         sfunc(nn+nxl+nxl*nyl) = sfunc(nn+nxl+nxl*nyl) + az(i,j,k)
C
         field(iv,i+1,j+1,k+1) = field(iv,i+1,j+1,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*az(i,j,k)
         sfunc(nn+1+nxl+nxl*nyl) = sfunc(nn+1+nxl+nxl*nyl) + az(i,j,k)
        enddo
        nn = nn + 1
       enddo
       nn = nn + nxl
      enddo
C
C  now divide the sums by the cell face areas
C
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) / sfunc(nn)
        enddo
       enddo
      enddo
C
C  return on successful completion
C
      wvel_interp = 1
      return
      end
C***********************************************************************
C
      integer function oth_interp(field,iv,nxl,nyl,nzl,nv,vols,nx,ny,nz)
C
C   this function interpolates the scalar field values to the nodal 
C   points based on a volume average with the cell volumes stored in the
c   vols routine.  This is an adaptation of the interpolation routine 
C   found in the mpgs conversion program
c
      parameter (nxt = 101, nyt = 101, nzt = 101)
      integer nxl,nyl,nzl,nx,ny,nz,nv,iv,nn
      real field(nv,nxl,nyl,nzl), vols(nx,ny,nz)
      real sfunc(nxt*nyt*nzt),vfunc(nxt*nyt*nzt)
C
C  first zero out the accumulation arrays and initialize the vfunc array
C
      print*,' in oth interp'
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
         vfunc(nn) = field(iv,i,j,k)
         field(iv,i,j,k) = 0.0
         sfunc(nn) = 0.0
        enddo
       enddo
      enddo
C
C  now step through the domain, using cell volume X velocity products
C  to average the values to each nodal point
C
      nn = 0
      do k = 1,nz
       do j = 1,ny
        do i = 1,nx
C        print*,i,j,k
         nn = nn + 1
         field(iv,i,j,k) = field(iv,i,j,k) + vfunc(nn+1+nxl+nxl*nyl)*
     *                                               vols(i,j,k)
         sfunc(nn) = sfunc(nn) + vols(i,j,k)
C
         field(iv,i+1,j,k) = field(iv,i+1,j,k)+vfunc(nn+1+nxl+nxl*nyl)*
     *                                               vols(i,j,k)
         sfunc(nn+1) = sfunc(nn+1) + vols(i,j,k)
C
         field(iv,i,j+1,k) = field(iv,i,j+1,k)+vfunc(nn+1+nxl+nxl*nyl)*
     *                                               vols(i,j,k)
         sfunc(nn+nxl) = sfunc(nn+nxl) + vols(i,j,k)
C
         field(iv,i,j,k+1) = field(iv,i,j,k+1)+vfunc(nn+1+nxl+nxl*nyl)*
     *                                               vols(i,j,k)
         sfunc(nn+nxl*nyl) = sfunc(nn+nxl*nyl) + vols(i,j,k)
C
         field(iv,i+1,j+1,k) = field(iv,i+1,j+1,k) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*vols(i,j,k)
         sfunc(nn+1+nxl) = sfunc(nn+1+nxl) + vols(i,j,k)
C
         field(iv,i+1,j,k+1) = field(iv,i+1,j,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*vols(i,j,k)
         sfunc(nn+1+nxl*nyl) = sfunc(nn+1+nxl*nyl) + vols(i,j,k)
C
         field(iv,i,j+1,k+1) = field(iv,i,j+1,k+1) + 
     *                          vfunc(nn+1+nxl+nxl*nyl)*vols(i,j,k)
         sfunc(nn+nxl+nxl*nyl) = sfunc(nn+nxl+nxl*nyl) + vols(i,j,k)
C
         field(iv,i+1,j+1,k+1) = field(iv,i+1,j+1,k+1) + 
     *                           vfunc(nn+1+nxl+nxl*nyl)*vols(i,j,k)
         sfunc(nn+1+nxl+nxl*nyl) = sfunc(nn+1+nxl+nxl*nyl) + vols(i,j,k)
c        if( k.eq. 10) then
c           print*,i,j,k,nn
c           print*, vfunc(nn+1+nxl+nxl*nyl),vols(i,j,k)
c           print*, nn,sfunc(nn),nn+1,sfunc(nn+1)
c           print*, nn+nxl,sfunc(nn+nxl),nn+nxl*nyl,sfunc(nn+nxl*nyl)
c           print*, nn+1+nxl,sfunc(nn+1+nxl)
c           print*, nn+1+nxl*nyl,sfunc(nn+1+nxl*nyl)
c           print*, nn+nxl+nxl*nyl,sfunc(nn+nxl+nxl*nyl)
c           print*, nn+1+nxl+nxl*nyl,sfunc(nn+1+nxl+nxl*nyl)
c        endif
        enddo
        nn = nn + 1
       enddo
       nn = nn + nxl
      enddo
C
C  now divide the sums by the cell face areas
C
      nn = 0
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         nn = nn + 1
c        if(sfunc(nn) .lt. 1.e-7) print*,i,j,k
         field(iv,i,j,k) = field(iv,i,j,k) / sfunc(nn)
        enddo
       enddo
      enddo
C
C  return on successful completion
C
      oth_interp = 1
      return
      end
 
C***************************************************************************
C
      integer function cart_uvel_extrap(field,iu,nxl,nyl,nzl,nv,xc,nx)
C
C   this function extrapolates the u velocty components to east-most
C   i = 1 face of the domain and to the west most i=nx face of the domain
C   by simple linear extrapolation for the cartesian case
C
      dimension field(nv,nxl,nyl,nzl),xc(0:nx)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*, ' in uvel_extrap'
      if(nxl .gt. 2) then
       do k = 2,nzl
        do j = 2,nyl
          dx21 = xc(1)-xc(0)
          dx32 = xc(2)-xc(1)
          dxr = dx21 / dx32
          field(iu,1,j,k) = field(iu,2,j,k) + dxr*field(iu,2,j,k) 
     1                                    - dxr*field(iu,3,j,k)
         enddo
        enddo
        do k = 2,nzl
         do j = 2,nyl
          dxnx1 = xc(nx)-xc(nx-1)
          dxnx2 = xc(nx-1)-xc(nx-2)
          dxr = dxnx1 / dxnx2
          field(iu,nxl,j,k) = field(iu,nxl-1,j,k) +
     1             dxr*field(iu,nxl-1,j,k) -dxr*field(iu,nxl-2,j,k)
         enddo
       enddo
      else
       do k = 2,nzl
         do j = 2,nyl
           field(iu,1,j,k) = field(iu,2,j,k)
           field(iu,nxl,j,k) = field(iu,nxl-1,j,k)
         enddo
       enddo
      endif
C
      cart_uvel_extrap = 1
      return
      end 
C***************************************************************************
C
      integer function cart_vvel_extrap(field,iv,nxl,nyl,nzl,nv,yc,ny)

C
C   this function extrapolates the v velocity components to south-most
C   j = 1 face of the domain and to the north most j = nyl face of the domain
C   by simple linear extrapolation
C
      dimension field(nv,nxl,nyl,nzl),yc(0:ny)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*,' in vvel extrap'
      if(nyl .gt. 2) then
       do k = 2,nzl
        do i = 2,nxl
          dy21 = yc(1)-yc(0) 
          dy32 = yc(2)-yc(1)
          dyr = dy21 / dy32
          field(iv,i,1,k) = field(iv,i,2,k) + dyr*field(iv,i,2,k) 
     1                                    - dyr*field(iv,i,3,k)
         enddo
        enddo
        do k = 2,nzl
         do j = 2,nyl
          dynyl1 = yc(ny)-yc(ny-1)
          dynyl2 = yc(ny-1)-yc(ny-2)
          dyr = dynyl1 / dynyl2
          field(iv,i,nyl,k) = field(iv,i,nyl-1,k)+
     1              dyr*field(iv,i,nyl-1,k) -dyr*field(iv,i,nyl-2,k)
         enddo
       enddo
      else
       do k = 2,nzl
         do i = 2,nxl
           field(iv,i,1,k) = field(iv,i,2,k)
           field(iv,i,nyl,k) = field(iv,i,nyl-1,k)
         enddo
       enddo
      endif
C
      cart_vvel_extrap = 1
      return
      end
C***************************************************************************
C
      integer function cart_wvel_extrap(field,iw,nxl,nyl,nzl,nv,zc,nz)
C
C   this function extrapolates the w velocty components to low-most
C   k = 1 face of the domain and to the high most k=nzl face of the domain
C   by simple linear extrapolation
C
      dimension field(nv,nxl,nyl,nzl),zc(0:nz)
C
C  start the loops at 2 since we stored the data shifted one up in all
C  directions
C
      print*,' in wvel extrap'
      if(nzl .gt. 2) then
       do j = 2,nyl
        do i = 2,nxl
          dz21 = zc(1)-zc(0)
          dz32 = zc(2)-zc(1)
          dzr = dz21 / dz32
          field(iw,i,j,1) = field(iw,i,j,2) + dzr*field(iw,i,j,2) 
     1                                    - dzr*field(iw,i,j,3)
         enddo
        enddo
        do j = 2,nyl
         do i = 2,nxl
          dznzl1 = zc(nz)-zc(nz-1)
          dznzl2 = zc(nz-1)-zc(nz-2)
          dzr = dznzl1 / dznzl2
          field(iw,i,j,nzl) = field(iw,i,j,nzl-1) + 
     1               dzr*field(iw,i,j,nzl-1) - dzr*field(iw,i,j,nzl-2)
         enddo
       enddo
      else
       do j = 2,nyl
         do i = 2,nxl
           field(iw,i,j,1) = field(iw,i,j,2)
           field(iw,i,j,nzl) = field(iw,i,j,nzl-1)
         enddo
       enddo
      endif
C
      cart_wvel_extrap = 1
      return
      end
C************************************************************************
C
      integer function field_reorder(field,nxl,nyl,nzl,nval,name_stored,
     *                                     nrep,nlast)
C
C   this function moves the nlast scalar field to the nrep postion and
C   reorders everything appropriately
C
      dimension field(nval,nxl,nyl,nzl)
      character*4 tname,name_stored(nlast)
      print*, 'in field reorder'
C
      do k = 1,nzl
       do j = 1,nyl
        do i = 1,nxl
         temp = field(nrep,i,j,k)
         field(nrep,i,j,k) = field(nlast,i,j,k)
         field(nlast,i,j,k) = temp
        enddo
       enddo
      enddo
C
C  now change the names array
c   
      tname = name_stored(nrep)
      name_stored(nrep) = name_stored(nlast)
      name_stored(nlast) = tname
c
      field_reorder = 1
      return
      end