
      MODULE types
      
!    symbolic names for 4, 2, and 1 byte integers
	
	INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
	INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
	INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)

! and for single and double precision reals
!	INTEGER, PARAMETER :: SP = KIND(1.0)
!	INTEGER, PARAMETER :: DP = KIND(1.0D0)
	INTEGER, PARAMETER :: SP = KIND(1.0D0)

      END MODULE types


	
module hpfprocessors
!HPF$ PROCESSORS PROC(4)
end module hpfprocessors


module diff_interface

INTERFACE
subroutine diff_x(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_x

subroutine diff_y(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_y

subroutine diff_z(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_z

END INTERFACE
end module diff_interface







subroutine inverse(g11,g12,g13,g22,g23,g33, &
                   gup11,gup12,gup13,gup22,gup23,gup33)

!----------------------------------------------------------------------------
!  Routine to calculate the inverse 3-metric. Returns the determinant of
!  the metric also. 
!----------------------------------------------------------------------------

use hpfprocessors
use types

real(SP), dimension(:,:,:),intent(IN)  :: g11,g12,g13,g22,g23,g33
real(SP), dimension(:,:,:),intent(OUT) :: gup11,gup12,gup13,gup22,gup23,gup33

!HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
!HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
!HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33

real(SP), dimension(size(g11,1),size(g11,2),size(g11,3)) :: detg
!HPF$ ALIGN WITH g11  :: detg

integer(I4B) ::  i, j
real(SP),parameter     :: zero=0.0d0, value=1.0d0


!   first calculate the determinant

      detg = g11*g22*g33 - g11*g23*g23  -  &
             g12*g12*g33 + g12*g13*g23  +  &
             g13*g12*g23 - g13*g13*g22


!if(ANY(detg==0)) then 
!  write(6,*) 'determinant = 0. Sorry, exit!'
!  stop
!endif 

      gup11 = (g22 * g33 - g23*g23) / detg
      gup12 = (g13 * g23 - g12*g33) / detg
      gup13 = (g12 * g23 - g13*g22) / detg
      gup22 = (g11 * g33 - g13*g13) / detg
      gup23 = (g13 * g12 - g11*g23) / detg
      gup33 = (g11 * g22 - g12*g12) / detg

return
end subroutine inverse 
