stitch_nest.f90

! Stitch nest
!
! Combines velocity fields from a coarse model and n nests at one level finer to
! update the surface positions and temperatures in the coarse model
!
! dwhipp - 06/11

program stitch_nest

use definitions
use threads

implicit none

include 'mpif.h'

! Variable type declaration
character (len=4) :: cistep,cnnest

double precision :: current_time,inc,step,total,xss,yss,zss
double precision :: xmin,xmax,ymin,ymax,zmin,zmax

integer :: err,i,ierr,iproc,istep,iter,j,material0,nnest,nproc

type (box),dimension(:),allocatable :: boxes
type (cloud),dimension(:),allocatable :: ss_cl
type (cloud) :: ls0_cl,ls_cl
type (cross_section),dimension(:),allocatable :: sections
type (face),dimension(6) :: cube_faces
type (material),dimension(:),allocatable :: mat
type (nest_info),dimension(:),allocatable :: nest
type (octreelsf) :: olsf
type (octreesolve) :: ls_osolve
type (octreev),dimension(:),allocatable :: ss_ov
type (octreev) ls0_ov,ls_ov
type (parameters) :: params
type (sheet),dimension(:,:),allocatable :: ss_surf,ss_surf0
type (sheet),dimension(:),allocatable :: ls0_surf,ls0_surf0,ls_surf,ls_surf0
type (thread) :: threadinfo
type (void) :: ls_vo
type (ziso) :: zi

!-------------------------------------------------------------------------------

!write (*,*) 'Program started'
call show_time (total,step,inc,0,'Program started$')

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

! Init some variables
current_time=0.d0

! Read/echo current time step and # of nests from command line
call getarg (1,cistep)
call getarg (2,cnnest)
write (*,*) 'Received command line arguments:'
write (*,*) 'istep = ',trim(cistep)
write (*,*) 'nnest = ',trim(cnnest)

! Convert character values to integers
read (cistep,'(i4)') istep
read (cnnest,'(i4)') nnest

! Read the coarse model input file
params%infile=trim('input_ls.txt')
call read_controlling_parameters (params,threadinfo,'main')
allocate (ls0_surf  (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf in main$')
allocate (ls0_surf0 (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf0 in main$')
allocate (mat (0:params%nmat),stat=err); if (err/=0) call stop_run ('Error alloc mat in main$')
allocate (nest (nnest),stat=err); if (err/=0) call stop_run ('Error alloc nest in main$')
if (params%nboxes.gt.0) then
  allocate (boxes(params%nboxes),stat=err) ; if (err/=0) call stop_run ('Error alloc boxes arrays$')
else
  allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
endif
if (params%nsections.gt.0) then
  allocate (sections(params%nsections),stat=err) ; if (err/=0) call stop_run ('Error alloc cross_section arrays$')
else
  allocate (sections(1),stat=err) ! necessary to avoid nil size argument
end if

call read_input_file (params,threadinfo,material0,mat,ls0_surf,boxes,sections,  &
                      cube_faces,nest(1))

params%irestart=0
! Modify params%restartfile to read previous restart file
if (istep > 1) then
  call int_to_char (cistep,4,istep-1)
  ! This may not always be the OUT directory...need to keep it general
  params%irestart=istep
  params%restartfile='LSOUT/time_'//cistep//'.bin'
endif

! Read the previous coarse model restart file if istep > 1, define surface,
! cloud and ov
call define_surface (params,ls0_surf,threadinfo,total,step,inc,current_time)
call define_cloud (ls0_cl,params,zi)
call define_ov (ls0_ov,params,threadinfo)

! Read the coarse model restart file from the current time step
call int_to_char (cistep,4,istep)
! As above, this may now always be in the OUT directory
params%restartfile='LSOUT/time_'//cistep//'.bin'

! Read the current coarse model restart file, define surface, cloud and ov
allocate (ls_surf  (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf in main$')
allocate (ls_surf0 (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf0 in main$')
call read_restart (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1))

deallocate(mat,boxes,sections)

allocate (ss_surf(params%ns,nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_surf in main$')
allocate (ss_surf0(params%ns,nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_surf in main$')
allocate (ss_cl(nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_cl in main$')
allocate (ss_ov(nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_ov in main$')

! Loop through fine model(s) and read their current restart files
do i=1,nnest
  ! Read the coarse model input file
  params%infile=trim('input_ss.txt')
  call read_controlling_parameters (params,threadinfo,'main')
  allocate (mat (0:params%nmat),stat=err); if (err/=0) call stop_run ('Error alloc mat in main$')
  if (params%nboxes.gt.0) then
    allocate (boxes(params%nboxes),stat=err) ; if (err/=0) call stop_run ('Error alloc boxes arrays$')
  else
    allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
  endif
  if (params%nsections.gt.0) then
    allocate (sections(params%nsections),stat=err) ; if (err/=0) call stop_run ('Error alloc cross_section arrays$')
  else
    allocate (sections(1),stat=err) ! necessary to avoid nil size argument
  endif

  call read_input_file (params,threadinfo,material0,mat,ss_surf(:,i),boxes,sections,  &
                        cube_faces,nest(i))

  ! Modify params%restartfile to read previous restart file
  params%restartfile='SSOUT/time_'//cistep//'.bin'

  call define_surface (params,ss_surf(:,i),threadinfo,total,step,inc,current_time)
  call define_cloud (ss_cl(i),params,zi)
  call define_ov (ss_ov(i),params,threadinfo)
  !call read_restart (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1))
enddo

do i=1,nnest
  xmin=nest(i)%xminls
  xmax=nest(i)%xminls+nest(i)%sselemx
  ymin=nest(i)%yminls
  ymax=nest(i)%yminls+nest(i)%sselemy
  zmin=nest(i)%zminls
  zmax=nest(i)%zminls+nest(i)%sselemz
  do j=1,ls_osolve%nnode
    xss=(ls_osolve%x(i)-nest(i)%xminls)/nest(i)%sselemx
    yss=(ls_osolve%y(i)-nest(i)%yminls)/nest(i)%sselemy
    zss=(ls_osolve%z(i)-nest(i)%zminls)/nest(i)%sselemz
    if(ls_osolve%x(i)<=xmax .and. ls_osolve%x(i)>=xmin .and. ls_osolve%y(i)<=ymax .and. ls_osolve%y(i)>=ymin .and. ls_osolve%z(i)<=zmax .and. ls_osolve%z(i)>=zmin) then
      call octree_interpolate_three (3,ss_ov(i)%octree,ss_ov(i)%noctree,ss_ov(i)%icon,       &
                                   ss_ov(i)%nleaves,ss_ov(i)%nnode,xss,yss,zss,         &
                                   ss_ov(i)%unode,ls_osolve%u(i),ss_ov(i)%vnode,ls_osolve%v(i),ss_ov(i)%wnode, &
                                   ls_osolve%w(i))
    endif
  enddo
enddo

! Create ov arrays for move_surface
ls_ov%noctree=ls_osolve%noctree
ls_ov%nnode=ls_osolve%nnode
ls_ov%nleaves=ls_osolve%nleaves

allocate (ls_ov%octree(ls_ov%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%octree in stitch_nest$')
allocate (ls_ov%x(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%x in stitch_nest$')
allocate (ls_ov%y(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%y in stitch_nest$')
allocate (ls_ov%z(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%z in stitch_nest$')
allocate (ls_ov%unode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%unode in stitch_nest$')
allocate (ls_ov%vnode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%vnode in stitch_nest$')
allocate (ls_ov%wnode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%wnode in stitch_nest$')
allocate (ls_ov%wnodeiso(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%wnodeiso in stitch_nest$')
allocate (ls_ov%icon(8,ls_ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%icon in stitch_nest$')

ls_ov%octree=ls_osolve%octree
ls_ov%x=ls_osolve%x
ls_ov%y=ls_osolve%y
ls_ov%z=ls_osolve%z
ls_ov%unode=ls_osolve%u
ls_ov%vnode=ls_osolve%v
ls_ov%wnode=ls_osolve%w
ls_ov%wnodeiso=ls_osolve%wiso
ls_ov%icon=ls_osolve%icon

write (*,*) 'Before calls to move surface'

do i=1,params%ns
  if (current_time+tiny(current_time).ge. ls_surf(i)%activation_time) then
    write (*,*) 'Before call to move surface ',i
    ! Need to make sure to use params%dt from the ls model here!!!
    call move_surface (ls0_surf(i),ls0_surf(i),0,0,ls_ov,params%dt,params,istep,i)
    write (*,*) 'After call to move surface ',i
  else
    ls_surf(i)%x=ls_surf(1)%x
    ls_surf(i)%y=ls_surf(1)%y
    ls_surf(i)%z=ls_surf(1)%z
    ls_surf(i)%xn=ls_surf(1)%xn
    ls_surf(i)%yn=ls_surf(1)%yn
    ls_surf(i)%zn=ls_surf(1)%zn
    ls_surf(i)%r=ls_surf(1)%r
    ls_surf(i)%s=ls_surf(1)%s
  endif

  ! We apply erosion
  if (material0.eq.0.and.params%erosion) then
    if (current_time+tiny(current_time).ge. ls_surf(i)%activation_time) then    
      call erosion (ls_surf(i),olsf,i,params)
    else
      if (iproc.eq.0) write(8,*) 'Surface ',i,' is attached to Surface0'
      ls_surf(i)%z=ls_surf(1)%z   
    endif
  end if
  if (iproc.eq.0) write (8,*) 'Min-Max z surf ',i,':',minval(ls_surf(i)%z),maxval(ls_surf(i)%z)

  allocate(ls_surf(i)%u(ls_surf(i)%nsurface))
  allocate(ls_surf(i)%v(ls_surf(i)%nsurface))
  allocate(ls_surf(i)%w(ls_surf(i)%nsurface))

enddo

! Need to do isostasy???

! Write new output file
call write_global_output (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1),'final')

deallocate(ls_surf,ls_surf0,ss_surf,ss_surf0)

call mpi_finalize (ierr)

! Sign off
write (*,*) 'Program execution complete'
end