define_cloud.f90

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

subroutine define_cloud (cl,params,bcdef,nest,density_str,threadinfo)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! if irestart=0, this routine allocates and creates the cloud of points present 
! in the system. Otherwise it reads from 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+icloud+surface informations, the routine 
! first reads dummy parameters until it gets to the real interesting 
! cloud information 

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

use definitions
!use mpi
use threads

implicit none

include 'mpif.h'

type (cloud) cl
type (parameters) params
type (bc_definition) bcdef
type (nest_info) :: nest
type (string) :: density_str(2**params%levelmax_oct,2**params%levelmax_oct)
type (thread) threadinfo

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision x,y,z,u,v,w,p,r,s,t,e2d,xlsf,xn,yn,zn,con,wiso,ev,epr,espr,yr
double precision current_time,activation_time,dilatr,fa,wcompact,cd,zi
double precision ssx,ssy,ssz,xls,yls,zls,eps,str_den,str_denp
logical inf,wlif,ip,isostasyout,isobcout,nestout,compactionout
integer err,iproc,nproc,ierr,i,j,k,nl,nf,nn,nr,kx,ky,kz,kt,icon,str_length
integer ioc,iface,ilsf,noctree,nnode,nface,nlsf,nleaves,ns,nt,restart_res,matnum
integer vermajor,verminor,verstat,verrev

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

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

eps=1.d-10

if (params%irestart.eq.0) then
   cl%np=0
   allocate (cl%x(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x in define_cloud$')
   allocate (cl%y(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y in define_cloud$')
   allocate (cl%z(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z in define_cloud$')
   allocate (cl%x0(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in define_cloud$')
   allocate (cl%y0(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y0 in define_cloud$')
   allocate (cl%z0(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z0 in define_cloud$')
   allocate (cl%x0mp(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0mp in define_cloud$')
   allocate (cl%y0mp(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y0mp in define_cloud$')
   allocate (cl%z0mp(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z0mp in define_cloud$')
   allocate (cl%strain(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%strain in define_cloud$')
   allocate (cl%lsf0(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%lsf0 in define_cloud$')
   allocate (cl%temp(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%temp in define_cloud$')
   allocate (cl%press(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%press in define_cloud$')
   allocate (cl%e2dp(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%e2dp in define_cloud$')
   allocate (cl%tag(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in define_cloud$')
   allocate (cl%matnum(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%matnum in define_cloud$')
   allocate (cl%ematnump(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%ematnump in define_cloud$')
   if (params%isobc) then
     allocate (bcdef%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc bcdef%zisodisp in define_cloud$')
     allocate (bcdef%zisoinc(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc bcdef%zisoinc in define_cloud$')
     bcdef%zisodisp=0.d0
     bcdef%zisoinc=0.d0
   endif

else
   open (19,file=trim(params%restartfile),status='old',form='unformatted')
   ! Read version number
   read (19) vermajor,verminor,verstat,verrev

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

   read (19) noctree, &
               nnode,   &
               nleaves, &
               nface,   &
               nlsf,    &
               cl%np,   &
               current_time     

   allocate (cl%x(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x in define_cloud$')
   allocate (cl%y(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y in define_cloud$')
   allocate (cl%z(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z in define_cloud$')
   allocate (cl%x0(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in define_cloud$')
   allocate (cl%y0(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y0 in define_cloud$')
   allocate (cl%z0(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z0 in define_cloud$')
   allocate (cl%x0mp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0mp in define_cloud$')
   allocate (cl%y0mp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y0mp in define_cloud$')
   allocate (cl%z0mp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z0mp in define_cloud$')
   allocate (cl%strain(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%strain in define_cloud$')
   allocate (cl%lsf0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%lsf0 in define_cloud$')
   allocate (cl%temp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%temp in define_cloud$')
   allocate (cl%press(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%press in define_cloud$')
   allocate (cl%e2dp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%e2dp in define_cloud$')
   allocate (cl%tag(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in define_cloud$')
   allocate (cl%matnum(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%matnum in define_cloud$')
   allocate (cl%ematnump(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%ematnump in define_cloud$')
   if (params%isobc) then
     allocate (bcdef%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc bcdef%zisodisp in define_cloud$')
     allocate (bcdef%zisoinc(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc bcdef%zisoinc in define_cloud$')
     bcdef%zisodisp=0.d0
     bcdef%zisoinc=0.d0
   endif

   read (19) (x,               &
                y,               &
                z,               &
                u,               &
                v,               &
                w,               &
                (xlsf,j=1,nlsf), &
                t,               &
                p,               &
                s,               &
                kx,              &
                ky,              &
                kz,              &
                kt,              &
                i=1,nnode)

   if (isostasyout) read (19) (wiso,i=1,nnode)
   if (compactionout) read (19) (wcompact,i=1,nnode)

   read (19) ((icon,k=1,8),epr,espr,con,e2d,ev,ip,dilatr,matnum,yr,fa,wlif,i=1,nleaves)

   if (compactionout) read (19) (cd,i=1,nleaves)

   read (19) (ioc,k=1,noctree)                    
   read (19) ((iface,k=1,9),i=1,nface)
   read (19) (nn,nl,nf,nr,inf,i=1,nnode)
   read (19) (nf,i=1,nface)
   do ilsf=1,nlsf                                  
      read (19) ns, &
                  activation_time, &
                  nt
      read (19) (r,  &
                   s,  &
                   x,  &
                   y,  &
                   z,  &
                   xn, &
                   yn, &
                   zn, &
                   u,v,w, &
                   i=1,ns)
      read (19) ((icon,k=1,3),i=1,nt)
   enddo
   
   read (19) (cl%x(i),      & 
                cl%y(i),      &
                cl%z(i),      &
                cl%x0(i),     &
                cl%y0(i),     & 
                cl%z0(i),     &
                cl%strain(i), &
                cl%lsf0(i),   &
                cl%temp(i),   &
                cl%press(i),  &
                cl%e2dp(i),   &
                cl%tag(i),    &
                cl%matnum(i), &
                cl%ematnump(i),&
                i=1,cl%np)

   ! correct for vertical exaggeration
   cl%z=cl%z/params%vex
   cl%z0=cl%z0/params%vex

   ! read isostasy basal displacement array (optional) - dwhipp 07.13
   if (isobcout) then
     read (19) restart_res
     if (params%isobc) then
       if (restart_res /= 2**params%levelmax_oct) then
         if (iproc==0) write (*,'(a)') 'Error: Restart file and current model have different resolutions for the'
         if (iproc==0) write (*,'(a)') '       isostasy basal displacement arrays!'
         if (iproc==0) write (*,'(a)') 'This is not supported and DOUAR cannot continue.'
         call stop_run ('Error: Different octree resolutions on restart in define_cloud$')
       endif
       read (19) ((bcdef%zisodisp(i,j),j=1,2**params%levelmax_oct+1),          &
                  i=1,2**params%levelmax_oct+1)
     else
       read (19) ((zi,j=1,restart_res+1),i=1,restart_res+1)
     endif
   elseif (params%isobc) then
     if (iproc==0) write (*,'(a)') 'Warning: Restarting without previous isostasy basal displacement data!'
     if (iproc==0) write (*,'(a)') '         Assuming zero basal displacement.'
   endif

   ! read nest info (optional) - dwhipp 07.13
   if (nestout) read(19) ssx,ssy,ssz,xls,yls,zls
   if (params%nest .and. nestout) then
     if (abs(ssx-nest%sselemx) > eps .or. abs(ssy-nest%sselemy) > eps .or.     &
         abs(ssz-nest%sselemz) > eps .or. abs(xls-nest%xminls) > eps .or.      &
         abs(yls-nest%yminls) > eps .or. abs(zls-nest%zminls) > eps) then         
       if (iproc==0) write (*,'(a)') 'Error: Nest dimensions differ on restart'
       if (iproc==0) write (*,'(a)') 'This is not supported and DOUAR cannot continue.'
       call stop_run ('Error: Different nest resolutions on restart in define_cloud$')
     endif
   endif

   ! read compaction density strings (optional) - dwhipp 07.13
   if (compactionout) then
     read (19) restart_res
     if (params%compaction) then
       if (restart_res /= 2**params%levelmax_oct) then
         if (iproc==0) write (*,'(a)') 'Error: Restart file and current model have different resolutions for the'
         if (iproc==0) write (*,'(a)') '       compaction string arrays!'
         if (iproc==0) write (*,'(a)') 'This is not supported and DOUAR cannot continue.'
         call stop_run ('Error: Different octree resolutions on restart in define_cloud$')
       endif
       do j=1,2**params%levelmax_oct
         do i=1,2**params%levelmax_oct
           read (19) str_length
           allocate (density_str(i,j)%density(str_length),stat=threadinfo%err)
           allocate (density_str(i,j)%densityp(str_length),stat=threadinfo%err)
           read (19) (density_str(i,j)%density(k),density_str(i,j)%densityp(k),  &
                      k=1,str_length)
         enddo
       enddo
     else
       do j=1,restart_res
         do i=1,restart_res
           read (19) str_length
           read (19) (str_den,str_denp,k=1,str_length)
         enddo
       enddo
     endif
   elseif (params%compaction) then
     if (iproc==0) write (*,'(a)') 'Warning: Restart file does not contain density string output!'
     if (iproc==0) write (*,'(a)') 'Assuming no previous sediment compaction.'
   endif

   close (19)
endif

if (iproc.eq.0) write (8,*) cl%np,' particles in cloud'

return
end

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