define_surface.f90

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

subroutine define_surface (params,surface,threadinfo,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
use threads

implicit none

include 'mpif.h'

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

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
     
integer ierr,iproc,nproc,is,noctree,nnode,nleaves,nface
integer nlsf,ncl
integer vermajor,verminor,verstat,verrev
logical isostasyout,isobcout,nestout,compactionout
character*7 ic
double precision 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
      call int_to_char(ic,2,is)
      if (iproc.eq.0) call show_time(total,step,inc,1,'surface '//ic(1:2)//'$')
         if (surface(is)%activation_time.ge.0.d0) then
            ! if activation_time ge 0 than the surface is simply made as a copy of the top (1st) surface
            if (is.eq.0) call stop_run ('First surface must be active at all times (activation_time<0)$')
            surface(is)%nsurface=surface(1)%nsurface
            surface(is)%levelt=surface(1)%levelt
            surface(is)%nt=surface(1)%nt
            surface(is)%itype=surface(1)%itype
            surface(is)%surface_type=surface(1)%surface_type
            surface(is)%rand=surface(1)%rand
            allocate (surface(is)%x(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
            allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
            allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
            allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
            allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
            allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
            allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%r),'dp',+1)
            allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','define_surface',size(surface(is)%s),'dp',+1)
            allocate (surface(is)%icon(3,surface(is)%nt),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','define_surface',size(surface(is)%icon),'int',+1)
            surface(is)%x=surface(1)%x
            surface(is)%y=surface(1)%y
            surface(is)%z=surface(1)%z
            surface(is)%xn=surface(1)%xn
            surface(is)%yn=surface(1)%yn
            surface(is)%zn=surface(1)%zn
            surface(is)%r=surface(1)%r
            surface(is)%s=surface(1)%s
            surface(is)%icon=surface(1)%icon(:,1:surface(is)%nt)
         else
            ! otherwise create_surface is invoked
!            surface(is)%nsurface=(2**surface(is)%levelt+1)*(2**surface(is)%levelt+1)
            surface(is)%nsurface=(2**surface(is)%levelt+2)*(2**surface(is)%levelt+2)   !opla
!            surface(is)%nt=(2**surface(is)%levelt)*(2**surface(is)%levelt)*2
            surface(is)%nt=(2**surface(is)%levelt +1)*(2**surface(is)%levelt +1)*2      ! opla
            allocate (surface(is)%x(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
            allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
            allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
            allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
            allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
            allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
            allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%r),'dp',+1)
            allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','define_surface',size(surface(is)%s),'dp',+1)
            allocate (surface(is)%icon(3,surface(is)%nt),stat=threadinfo%err)
            if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','define_surface',size(surface(is)%icon),'int',+1)
            call create_surf(surface(is),is,params%debug)
            surface(is)%z=surface(is)%z/params%vex
            surface(is)%zn=surface(is)%zn/params%vex
         endif
      current_time=0.d0
   end do
else
   open (9,file=trim(params%restartfile),status='old',form='unformatted')
   ! Read version number
   read (9) vermajor,verminor,verstat,verrev
   if (iproc==0) then
     if (vermajor/=params%vermajor .or. verminor/=params%verminor) then
       write(*,'(a,i1,a,i1,a,i1,a,i1)') 'DOUAR-WSMP v',params%vermajor,'.',params%verminor,' cannot be restarted using output from v',vermajor,'.',verminor
       write(*,'(a)') 'because the output formats are not compatible. To restart with this output'
       write(*,'(a,i1,a,i1,a)') 'please use DOUAR-WSMP v',params%vermajor,'.',params%verminor,'.'
       stop
     endif
   endif

   ! Read optional output flags
   read (9) isostasyout,isobcout,nestout,compactionout

   read (9) noctree,nnode,nleaves,nface,nlsf,ncl,current_time

   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)

   if (isostasyout) read (9)
   if (compactionout) read (9)

   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)
   read (9)

   if (compactionout) read (9)

   read (9)
   read (9)

   read (9)
   read (9)
   read (9)
   read (9)
   read (9)

   read (9)

   ! 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=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
      allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
      allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
      allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
      allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
      allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
      allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%x),'dp',+1)
      allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','define_surface',size(surface(is)%x),'dp',+1)
      allocate (surface(is)%icon(3,surface(is)%nt),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','define_surface',size(surface(is)%icon),'int',+1)

      read (9) surface(is)%r
      read (9) surface(is)%s
      read (9) surface(is)%x
      read (9) surface(is)%y
      read (9) surface(is)%z
      read (9) surface(is)%xn
      read (9) surface(is)%yn
      read (9) surface(is)%zn
      read (9)
      read (9)
      read (9)
      read (9) surface(is)%icon

! correct surface for vertical exageration

      surface(is)%z=surface(is)%z/params%vex
      surface(is)%zn=surface(is)%zn/params%vex

   end do
   close (9)
end if

return
end

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