program eps_modo
use sep
use sep3d_struct_mod
use basic_util_mod
type(sep3d) :: dat
real, allocatable :: eps(:,:,:)
integer :: n(5),i1,i2,i3
real :: o(5),d(5),eps0,eps_grad,z,eps_2

call sep_init(SOURCE)
call init_sep3d("in",dat,"INPUT")
call from_param("zero",eps0,1.)
call from_param("pow1",eps_grad,0.)
call from_param("pow2",eps_2,0.)
n=1
n(1:sep3d_ndims(dat))=dat%n(1:sep3d_ndims(dat))
o(1)=dat%o(1)
d(1)=dat%d(1)

allocate(eps(n(1),n(2),n(3)))
do i3=1,n(3)
  do i2=1,n(2)
    do i1=1,n(1)
      z=o(1)+d(1)*(i1-1)
      eps(i1,i2,i3)=eps0+eps_grad*z+z**2*eps_2
    end do
  end do
end do
do i5=1,n(5)
do i4=1,n(5)
call srite("out",eps,size(eps)*4)
end do
end do

end program






