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

subroutine write_global_output (params,istep,iter,current_time,osolve,ov, &
                                vo,surface,cl,bcdef,nest,outputtype)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! subroutine to generate large global output files in subdirectory OUT
! the files contain:
! - array lengths
! - nodal values
! - icon array
! - octree array
! - bad faces array
! - void arrays
! - surface information
! istep is time step
! osolve is solve octree
! surface are surfaces
! ns is number of surfaces
! see code for details and order of output
! Note that if you modify this routine, you may also need to modify
! define_ov.f90, define_cloud.f90, and define_surface.f90

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

use definitions
!use mpi

implicit none

include 'mpif.h'

type (parameters) params
integer istep,iter
double precision current_time
type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (sheet) surface(params%ns)
type (cloud) cl
type (bc_definition) bcdef
type (nest_info) :: nest
character*5 outputtype

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision,dimension(:),allocatable :: compaction_density,wcompact
integer i,j,is,k,err
integer iproc,nproc,ierr
character*4 cistep,citer
character(len=5) :: outdir

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

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

if (iproc.eq.0) then
   if (params%nest) then
     outdir=nest%ssoutdir
   else
     outdir='OUT'
   endif
   call int_to_char (cistep,4,istep)
   select case (outputtype)
      case ('debug')
         call int_to_char (citer,4,iter)
         open (9,file=trim(outdir)//'/time_'//cistep//'_'//citer//'.bin',status='unknown',form='unformatted')
      case ('final')
         open (9,file=trim(outdir)//'/time_'//cistep//'.bin',status='unknown',form='unformatted')
!         open (1234,file='OUT/time_'//cistep//'.dat',status='replace')
      case default
         call stop_run ('pb with argument in write_global_output$')
   end select

   ! Create temporary arrays for (optional) compaction output
   allocate (compaction_density(osolve%nleaves),stat=err)
   allocate (wcompact(osolve%nnode),stat=err)             
   if (params%compaction) then
     do i=1,osolve%nleaves
       compaction_density(i)=osolve%compaction_density(i)
     enddo
     do i=1,osolve%nnode
       wcompact(i)=osolve%wcompact(i)
     enddo
     !do i=1,10
     !  write (*,*) 'compaction_density(i): ',compaction_density(i)
     !  write (*,*) 'osolve%compaction_density(i): ',osolve%compaction_density(i)
     !  write (*,*) ''
     !enddo
     !compaction_density=osolve%compaction_density
     !wcompact=osolve%wcompact
   else
     compaction_density=0.d0
     wcompact=0.d0
   endif

   ! write version number
   write (9) params%vermajor,params%verminor,params%verstat,params%verrev

   ! write array lengths
   write (9) osolve%octree(3), &
               osolve%nnode,     &
               osolve%nleaves,   &
               osolve%nface,     &
               osolve%nlsf,      &
               cl%np,            &
               current_time

   ! write info on octree solve nodes (x,y,z,u,v,w,lsf,temp)
   write (9) (osolve%x(i),                           &
                osolve%y(i),                         &  
                osolve%z(i)*params%vex,              &
                osolve%u(i),                         &
                osolve%v(i),                         & 
                osolve%w(i),                         &
                osolve%wiso(i),                      &
                wcompact(i),                         &
                (osolve%lsf(i,j),j=1,osolve%nlsf),   &
                osolve%temp(i),                      &
                ov%temporary_nodal_pressure(i),      &
                osolve%strain(i),                    &
                osolve%kfix((i-1)*3+1),              &
                osolve%kfix((i-1)*3+2),              &
                osolve%kfix((i-1)*3+3),              &
                osolve%kfixt(i),                     &
                i=1,osolve%nnode)

   ! write icon array
   write (9) ((osolve%icon(k,i),k=1,8),     &
                 osolve%pressure(i),        &
                 osolve%spressure(i),       &
                 osolve%crit(i),            &
                 osolve%e2d(i),             &
                 osolve%eviscosity(i),      &
                 osolve%is_plastic(i),      &
                 osolve%dilatr(i),          &
                 osolve%matnum(i),          &
                 osolve%yield_ratio(i),     &
                 osolve%frict_angle(i),     &
                 compaction_density(i),     &
                 ov%whole_leaf_in_fluid(i), &
                 i=1,osolve%nleaves)

   ! write octree information
   write (9) (osolve%octree(i),i=1,osolve%octree(3))

   ! write bad face information
   write (9) ((osolve%iface(k,i),k=1,9),i=1,osolve%nface)

   ! write void information
   write (9) (vo%node(i),    &
                vo%leaf(i),    &
                vo%ftr(i),     &
                vo%rtf(i),     &
                vo%influid(i), &
                i=1,osolve%nnode)

   ! write bad faces information
   write (9) (vo%face(i),i=1,osolve%nface)

   ! write surface information (r,s,x,y,z,xn,yn,zn)
   do is=1,osolve%nlsf
      write (9) surface(is)%nsurface,          &
                  surface(is)%activation_time, &
                  surface(is)%nt
      select case (outputtype)
      case ('debug')
      write (9) (surface(is)%r(i),  &
                   surface(is)%s(i),  &
                   surface(is)%x(i),  &
                   surface(is)%y(i),  &
                   surface(is)%z(i)*params%vex,  &
                   surface(is)%xn(i) ,&
                   surface(is)%yn(i), &
                   surface(is)%zn(i)*params%vex, &
                   i=1,surface(is)%nsurface)
      case ('final')
      write (9) (surface(is)%r(i),  &
                   surface(is)%s(i),  &
                   surface(is)%x(i),  &
                   surface(is)%y(i),  &
                   surface(is)%z(i)*params%vex,  &
                   surface(is)%xn(i) ,&
                   surface(is)%yn(i), &
                   surface(is)%zn(i)*params%vex, &
                   surface(is)%u(i),  &
                   surface(is)%v(i),  &
                   surface(is)%w(i),  &
                   i=1,surface(is)%nsurface)
      end select
      write (9) ((surface(is)%icon(j,i),j=1,3),i=1,surface(is)%nt)
   enddo

   ! write cloud information
   write (9) (cl%x(i),      &
                cl%y(i),      &
                cl%z(i)*params%vex,      &
                cl%x0(i),     &
                cl%y0(i),     &
                cl%z0(i)*params%vex,     &
                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)

   ! write isostasy basal displacement array - dwhipp 11/09
   write (9) 2**params%levelmax_oct
   if (params%isobc) then
     write (9) ((bcdef%zisodisp(i,j)+surface(osolve%nlsf)%sp01,                   &
               j=1,2**params%levelmax_oct+1),i=1,2**params%levelmax_oct+1)
   else
     write (9) ((0.d0+surface(osolve%nlsf)%sp01,j=1,2**params%levelmax_oct+1), &
               i=1,2**params%levelmax_oct+1)
   end if

   ! write nested model info - dwhipp 05/11
   if (params%nest) then
     write(9) nest%sselemx,nest%sselemy,nest%sselemz,nest%xminls,nest%yminls,&
              nest%zminls
   else
     write(9) 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0
   endif

   close (9)

   deallocate (compaction_density)
   deallocate (wcompact)

end if

return
end

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