!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 threads
use definitions
implicit none
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 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,is_plas
character*7 ic
double precision s,p,t,x,y,z,u,v,w,wiso,xlsf,e2d,don,con,activation_time,dilatr
double precision epr,espr,evisc
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
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)
call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%r),'dp',+1)
allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
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)
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)
call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%r),'dp',+1)
allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
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)
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 (9) noctree, &
nnode, &
nleaves, &
nface, &
nlsf, &
ncl, &
current_time
read (9) (x, &
y, &
z, &
u, &
v, &
w, &
wiso, &
(xlsf,j=1,nlsf),&
t, &
p, &
s, &
kx, &
ky, &
kz, &
kt, &
i=1,nnode)
read (9) ((icon,k=1,8),epr,espr,con,e2d,evisc,is_plas,dilatr,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=threadinfo%err)
call heap (threadinfo,'surface(is)%x','define_surface',size(surface(is)%x),'dp',+1)
allocate (surface(is)%y(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%y','define_surface',size(surface(is)%y),'dp',+1)
allocate (surface(is)%z(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%z','define_surface',size(surface(is)%z),'dp',+1)
allocate (surface(is)%xn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%xn','define_surface',size(surface(is)%xn),'dp',+1)
allocate (surface(is)%yn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%yn','define_surface',size(surface(is)%yn),'dp',+1)
allocate (surface(is)%zn(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%zn','define_surface',size(surface(is)%zn),'dp',+1)
allocate (surface(is)%r(surface(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface(is)%r','define_surface',size(surface(is)%x),'dp',+1)
allocate (surface(is)%s(surface(is)%nsurface),stat=threadinfo%err)
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)
call heap (threadinfo,'surface(is)%icon','define_surface',size(surface(is)%icon),'int',+1)
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)
! 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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|