subroutine ob_freeze(g_,gp1_,k_,kp1_)

use tensors
type(GridFunction), dimension(3,3)     :: g_,k_,gp1_,kp1_
integer i,j
integer nx,ny,nz
integer, dimension(3) :: shp

! why do you need it here?
  shp = shape(g_(1,1)%d)
  nx = shp(1) 
  ny = shp(2) 
  nz = shp(3) 

    do i=1,3
      do j=i,3
	gp1_(i,j)%d(1 ,2:Ny-1,2:Nz-1) = g_(i,j)%d(1 ,2:Ny-1,2:Nz-1)
	gp1_(i,j)%d(Nx,2:Ny-1,2:Nz-1) = g_(i,j)%d(Nx,2:Ny-1,2:Nz-1)
	gp1_(i,j)%d(: ,1     ,2:Nz-1) = g_(i,j)%d(: ,1     ,2:Nz-1)
	gp1_(i,j)%d(: ,Ny    ,2:Nz-1) = g_(i,j)%d(: ,Ny    ,2:Nz-1)
	gp1_(i,j)%d(: ,:     ,1     ) = g_(i,j)%d(: ,:     ,1     )
	gp1_(i,j)%d(: ,:     ,Nz    ) = g_(i,j)%d(: ,:     ,Nz    )
			 !
	kp1_(i,j)%d(1 ,2:Ny-1,2:Nz-1) = k_(i,j)%d(1 ,2:Ny-1,2:Nz-1)
	kp1_(i,j)%d(Nx,2:Ny-1,2:Nz-1) = k_(i,j)%d(Nx,2:Ny-1,2:Nz-1)
	kp1_(i,j)%d(: ,1     ,2:Nz-1) = k_(i,j)%d(: ,1     ,2:Nz-1)
	kp1_(i,j)%d(: ,Ny    ,2:Nz-1) = k_(i,j)%d(: ,Ny    ,2:Nz-1)
	kp1_(i,j)%d(: ,:     ,1     ) = k_(i,j)%d(: ,:     ,1     )
	kp1_(i,j)%d(: ,:     ,Nz    ) = k_(i,j)%d(: ,:     ,Nz    )
      end do
    end do
end
