!/*@@
!
! $Id: initcn.f90,v 1.17 1997/07/21 13:52:05 haupt Exp $
!
!@@*/

subroutine initcn(params,m,d,nl,t, &
     & g_,k_,alpha_,beta_,mask_,mask_old_)
  use tensors
  use tempstorage
  use input_params
  implicit none
  real*8, dimension(100)                 :: params
  real*8, dimension(:)                   :: d
  integer, dimension(:)                  :: m
  integer                                :: nl
  integer                                :: t
  type(GridFunction)                     :: alpha_,mask_,mask_old_
  type(GridFunction), dimension(3)       :: beta_
  type(GridFunction), dimension(3,3)     :: g_,k_
  intent(IN)    :: t
  intent(INOUT) :: alpha_,mask_,mask_old_,beta_,g_,k_, &
       & m,d,nl,params

  type(GridFunction), pointer                     :: temp1_
  type(GridFunction), pointer                     :: temp2_,temp3_,temp4_
  type(GridFunction), dimension(:), pointer       :: tempv_,coords_
  type(GridFunction), dimension(:,:), pointer     :: tempt_
  type(GridFunction), dimension(:,:,:), pointer   :: c_
  
  integer compute_xyz
!DAGH$ user(INTEGER compute_xyz;)


  interface
     subroutine periodic_obc_symm2(params,gf)
       use tensors
       use input_params
       implicit none
       type(GridFunction),dimension(3,3),intent(INOUT)  :: gf
       real*8, dimension(100),intent(IN)                :: params
     end subroutine periodic_obc_symm2

     subroutine remask( g_, k_, mask_, mask_old_)
       use input_params
       use tensors
       implicit none
       type(GridFunction),dimension(3,3) :: g_, k_
       type(GridFunction)                :: mask_, mask_old_
     end subroutine remask

     subroutine horizon(x,y,z,r,params,mask)
       use input_params
       real*8,dimension(:,:,:)              :: mask
       real*8,dimension(100)                :: params
       real*8,dimension(:,:,:)              :: x,y,z,r
       real*8 s
     end subroutine horizon

     subroutine flags(params,mask,m,d,nl)
       use input_params
       implicit none
       integer                       :: nl
       real*8,dimension(100)         :: params
       real*8,dimension(:,:,:)       :: mask
       real*8,dimension(:)           :: d
       integer,dimension(:)          :: m
     end subroutine flags
     
     subroutine setshift(params,g_,k_,coords_,tempv_,mask,beta_,compute_xyz)
       use input_params
       use tensors
       real*8, dimension(100)               :: params
       type(GridFunction),dimension(3)      :: beta_,coords_,tempv_
       type(GridFunction),dimension(3,3)    :: g_, k_
       real*8,dimension(:,:,:)              :: mask
       integer                              :: compute_xyz
       intent(IN)    :: compute_xyz, params, g_, k_, mask
       intent(INOUT) :: beta_, coords_, tempv_
     end subroutine setshift

     subroutine setlapse(params,g_,beta_,coords_,tempv_,mask,alpha,compute_xyz)
       use input_params
       use tensors
       implicit none
       real*8, dimension(100)               :: params
       type(GridFunction), dimension(3)     :: beta_,tempv_,coords_
       type(GridFunction),dimension(3,3)    :: g_
       real*8,dimension(:,:,:)              :: mask, alpha
       integer                              :: compute_xyz
       intent(IN)    :: g_,beta_,params,mask,compute_xyz
       intent(OUT)   :: alpha
       intent(INOUT) :: coords_,tempv_
     end subroutine setlapse

     subroutine initial(params,tempt_,g_,k_,c_,tempv_)
       use input_params
       use tensors
       implicit none
       real*8, dimension(100) :: params
       type(GridFunction), dimension(3)   :: tempv_ 
       type(GridFunction), dimension(3,3) :: g_,k_,tempt_
       type(GridFunction), dimension(3,3,3):: c_
     end subroutine initial
     
  end interface
  
  params(P_CURRENTTIME) = 0.d0
!DAGH$ user(params[P_CURRENTTIME]=0;)

  !----------------------------
  ! Compute initial g and k.
  !----------------------------
  c_     => tempconnection()
  tempt_ => tempsymmtensor()
  tempv_ => tempvector()
  call initial(params,tempt_,g_,k_,c_,tempv_)
  call freevector(tempv_)
  call freesymmtensor(tempt_)
  call freeconnection(c_)
  
  !----------------------------------------
  ! Compute mask
  !----------------------------------------
  temp1_   => tempscalar()
  temp2_   => tempscalar()
  temp3_   => tempscalar()
  temp4_   => tempscalar()
  call horizon(temp1_%d,temp2_%d,temp3_%d,temp4_%d,params,mask_%d)
  call freescalar(temp1_)
  call freescalar(temp2_)
  call freescalar(temp3_)
  call freescalar(temp4_)
  mask_old_%d = mask_%d
!DAGH$ equals(mask_old_,mask_)

  !-------------------------------------------------\
  ! Use periodic boundary conditions                 >
  !-------------------------------------------------/
  if(params(P_OUTER_BOUNDARY).eq. 3.d0) then
     call periodic_obc_symm2(params,g_)
     call periodic_obc_symm2(params,k_)
  end if
  
  !----------------------------------------
  ! Set variables in interior of hole
  !----------------------------------------
  call remask(g_,k_,mask_,mask_old_) 
  
  !----------------------------------------------------------------
  ! Set up 1d arrays for taking derivatives near the hole boundary
  !----------------------------------------------------------------
  call flags(params,mask_%d,m,d,nl)
  write(33,*) 'total # of boundary points =',nl

  !--------------------------------------
  ! Set initial shift
  !--------------------------------------
  compute_xyz=1
!DAGH$ user(compute_xyz=1;)

  coords_ => tempvector()
  tempv_  => tempvector()
  call setshift(params,g_,k_,coords_,tempv_,mask_%d,beta_,compute_xyz)
  call freevector(tempv_)
  call freevector(coords_)
  !--------------------------------------
  ! Set initial lapse
  !--------------------------------------
  tempv_  => tempvector()
  coords_ => tempvector()
  call setlapse(params,g_,beta_,coords_,tempv_,mask_%d,alpha_%d,compute_xyz)
  call freevector(tempv_)
  call freevector(coords_)

end subroutine initcn

