fortran 90 trilinear interpolation

[sithif]

Привет,

написал 3D линейную интерполяцию для функции вида $f=f(x,y,z)$, которая задана таблично $x, y, z, f$.
Про сам метод есть, например, здесь paulbourke.net/miscellaneous/interpolation/index.html.
Вроде все работает для разного числа точек и для разного шага по направлениям.
Но мне нужно сделать быстрее код, так как эта функция используется очень часто.
Мне интересно на сколько "долгие" стейтменты типа,

xmin = minval(array(1,:))

Может мне лучше отказаться от них вообще и задавать их руками там же где и nx, ny, nz?

Используется синтаксис Fortran

 
! common data
module array_mod
        real, dimension(:,:), allocatable :: array
        ! number of greed points
        integer :: nx = 5
        integer :: ny = 10
        integer :: nz = 15
end module array_mod

! interface for data manipulation       
module field_in_mod
        interface field_in_inter
                subroutine field_in(i,field)
                        use dynamic
                        implicit none
                        integer :: i, k
                        character(7) :: field
                end subroutine field_in
        end interface field_in_inter
end module field_in_mod

! interpolation function interface
module trilinear_mod
        interface trilinear_int
                function trilinear(x,y,z)
                        use array_mod                   
                        implicit none
                        real :: x,y,z,dx,dy,dz,xmin,xmax,ymin,ymax,zmin,zmax
                        real :: trilinear
                        integer :: a,b,c,d
                end function trilinear
        end interface trilinear_int
end module trilinear_mod

! test program 
program test
        ! use modules
        use array_mod
        use field_in_mod
        use trilinear_mod
        implicit none
        integer :: i = 4  ! number of rows
        integer :: j      ! number of lines
        character(7) :: in = "in01.in" ! file name
        integer :: a,b,c 
        real, dimension(3) :: coord=(/3.1*0.15, 8.2*0.13, 7.4*0.17/)
        
        j = nx*ny*nz

        ! create input file
        open(unit=12, file=in, form='formatted')        
        do a=1,5,1
                do b=1,10,1
                        do c=1,15,1
                                write(12,fmt=*) a*0.15, b*0.13, c*0.17, 5.0
                        end do
                end do
        end do
        close(12)
        
        allocate(array(i,j))
        ! create i x j array
        call field_in(j,in)
        ! 3d interpolation
        write(*,*) trilinear(coord(1), coord(2), coord(3))
        deallocate(array)
        
end program test

function trilinear(x,y,z)
        use array_mod   
        implicit none
        ! paulbourke.net/miscellaneous/interpolation/index.html
        real :: x,y,z,dx,dy,dz,xmin,xmax,ymin,ymax,zmin,zmax
        real :: trilinear
        integer :: a,b,c,d
        
        xmin = minval(array(1,:))       
        xmax = maxval(array(1,:))       
        ymin = minval(array(2,:))       
        ymax = maxval(array(2,:))       
        zmin = minval(array(3,:))       
        zmax = maxval(array(3,:))       

        dx = (xmax-xmin)/(nx-1)
        dy = (ymax-ymin)/(ny-1)
        dz = (zmax-zmin)/(nz-1) 
        
        a = int((x-xmin)/dx)  ! (a)-time x value has changed
        b = int((y-ymin)/dy)  ! (b)-time y value has changed
        c = int((z-zmin)/dz)  ! (c)-time z value has changed

        d = ny*nz*a + nz*b + c + 1 ! current file position      

        x = x-(a+1)*dx
        y = y-(b+1)*dy 
        z = z-(c+1)*dz 
        
        !V000|V100|V010|V001|V101|V011|V110|V111

        trilinear   =   1.0/(dx*dy*dz)*(array(4,d)*(dx-x)*(dy-y)*(dz-z) + &
                        array(4,d+ny*nz)*x*(dy-y)*(dz-z) + &
                        array(4,d+nz)*(dx-x)*y*(dz-z) + &
                        array(4,d+1)*(dx-x)*(dy-y)*z + &
                        array(4,d+nz*ny+1)*x*(dy-y)*z + &
                        array(4,d+nz+1)*(dx-x)*y*z + &
                        array(4,d+ny*nz+nz)*x*y*(dz-z) + &
                        array(4,d+ny*nz+nz+1)*x*y*z)    

end function trilinear

subroutine field_in(i, field)
        use array_mod
        implicit none
        integer :: i, k
        character(7) :: field
        open(unit=10, file=field, form='formatted')
        ! suggestion => do do, field_in(i,j,field) in new version
        do k=1,i,1
                read(10,*)      array(1,k), & ! x
                                array(2,k), & ! y
                                array(3,k), & ! z
                                array(4,k)
                                
        end do
        close(10)
end subroutine field_in