!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DEFINE_SURFACE    Nov. 2006                                     |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_surface (params,surface,total,step,inc,current_time)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! if irestart=0, this routine allocates and creates the ns surfaces present 
! in the system. Otherwise, it reads form a user supplied file name the 
! surfaces as they were at the end of a previous run. In this case, since the 
! run output files contain all the octree+lsf+cloud+surface informations, 
! the routine first reads dummy parameters until it gets to the real 
! interesting surface information .

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
!use mpi

implicit none

include 'mpif.h'

type(parameters) params
type (sheet),dimension(params%ns)::surface
double precision total, step, inc, current_time

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
     
integer err,ierr,iproc,nproc,iface,ioc,i,j,k,is,noctree,nnode,nleaves,nface,nt
integer icon,inn,inl,nlsf,ncl,nf,inr,kx,ky,kz,kt
logical inf,wlif
character*7 ic
double precision s,p,t,x,y,z,u,v,w,xlsf,e2d,don,con,activation_time

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

if (params%irestart.eq.0) then
   do is=1,params%ns
      write (ic(1:2),'(i2)') is
      if (iproc.eq.0) call show_time(total,step,inc,1,'surface '//ic(1:2)//'$')
      surface(is)%nsurface=(2**surface(is)%levelt+1)*(2**surface(is)%levelt+1)
      surface(is)%nt=(2**surface(is)%levelt)*(2**surface(is)%levelt)*2
      allocate (surface(is)%x(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%x in define_surface$')
      allocate (surface(is)%y(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%y in define_surface$')
      allocate (surface(is)%z(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%z in define_surface$')
      allocate (surface(is)%xn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%xn in define_surface$')
      allocate (surface(is)%yn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%yn in define_surface$')
      allocate (surface(is)%zn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%zn in define_surface$')
      allocate (surface(is)%r(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%r in define_surface$')
      allocate (surface(is)%s(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface%s in define_surface$')
      allocate (surface(is)%icon(3,surface(is)%nt),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface%icons in define_surface$')
      call create_surf(surface(is),is,params%debug)
      current_time=0.d0
   end do
else
   open (9,file=trim(params%restartfile),status='old')   
   read (9,*) noctree, &
              nnode,   &
              nleaves, &
              nface,   &
              nlsf,    &
              ncl,     &
              current_time 
   read (9,*) (x,              &
               y,              &
               z,              & 
               u,              &
               v,              &
               w,              &
               (xlsf,j=1,nlsf),&
               t,              &
               p,              &
               s,              &
               kx,             &
               ky,             &
               kz,             &
               kt,             &
               i=1,nnode)
   read (9,*) ((icon,k=1,8),con,e2d,wlif,i=1,nleaves)
   read (9,*) (ioc,k=1,noctree)                    
   read (9,*) ((iface,k=1,9),i=1,nface)
   read (9,*) (inn, &
               inl, &
               nf,  &
               inr, &
               inf, &
               i=1,nnode)
   read (9,*) (nf,i=1,nface)
   ! read surface information
   do is=1,nlsf                                    
      read (9,*) surface(is)%nsurface, &
                 activation_time,      &
                 surface(is)%nt
      allocate (surface(is)%x(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%x in define_surface$')
      allocate (surface(is)%y(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%y in define_surface$')
      allocate (surface(is)%z(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%z in define_surface$')
      allocate (surface(is)%xn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%xn in define_surface$')
      allocate (surface(is)%yn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%yn in define_surface$')
      allocate (surface(is)%zn(surface(is)%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc surface%zn in define_surface$')
      allocate (surface(is)%r(surface(is)%nsurface), stat=err); if (err.ne.0) call stop_run ('Error alloc surface%r in define_surface$')
      allocate (surface(is)%s(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface%s in define_surface$')
      allocate (surface(is)%icon(3,surface(is)%nt),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface%icons in define_surface$')

      read (9,*) (surface(is)%r(i),  &
                  surface(is)%s(i),  &
                  surface(is)%x(i),  &
                  surface(is)%y(i),  &
                  surface(is)%z(i),  &
                  surface(is)%xn(i), &
                  surface(is)%yn(i), &
                  surface(is)%zn(i), &
                  u,v,w,             & 
                  i=1,surface(is)%nsurface)
      read (9,*) ((surface(is)%icon(k,i),k=1,3),i=1,surface(is)%nt)

   end do
   close (9)
end if

return
end

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|