!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! $BI=Bj(B: 2 $B<!85I9>2$N2r@O2r(B
! $BMzNr(B: 2003/11/06
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program analysis

  use prm , only : IMAX, JMAX, dx
  use com_var, only : get_frame, x
  use prmphys, only : s_a, rho, grav, A_con, Npw
  use Tbound, only: get_Ms_con, Ms

  implicit none

  
  real(8)  :: VAR(1:IMAX, JMAX)
  real(8)   :: FLUX(0:IMAX, JMAX)
  real(8)   :: alpha(1:IMAX, JMAX)
  
  integer :: i, j 

  call get_frame
  call get_Ms_con
  
  write(6,*)x
  alpha  = (20 * Ms*s_a/((rho * grav)**Npw * A_con*s_a))**(1D0/8D0)
  write(6, *)alpha(1,2)
  
  open(10, file='state.dat', form='unformatted', access='direct'&
       & ,recl=4*IMAX*JMAX)
  do j = 1, JMAX-1
     ! x=1 $B$+$i(B $BNN0h$NCf?4(B
     do i = 1, IMAX/2+1
        VAR(i, j)  =  alpha(i, j) * &
             &  (x(IMAX/2+1)**(4.0D0/3.0D0) - (x(IMAX/2+1)- x(i))**(4.0D0/3.0D0))**(3.0D0/8.0D0)!+ alpha(i, j) *x(IMAX/2+1)**(1.0D0/2.0D0)


        FlUX(i, j) = Ms(i, j)*(x(i)-x(IMAX/2+1)+dx/2.0D0)
        
     end do
     do i = IMAX/2+1, IMAX
        VAR(i, j)  = alpha(i, j) * &
&         (((x(IMAX)-x(IMAX/2+1)))**(4.0D0/3.0D0) -(x(i)-x(IMAX/2+1))**(4.0D0/3.0D0))**(3.0D0/8.0D0)

        FlUX(i, j) = Ms(i, j)*(x(i)-x(IMAX/2+1)+dx/2.0D0) 
     end do
  end do

  Flux(0, 3) = 0.0D0
  Flux(IMAX, 3) = 0.0D0 
  Flux(IMAX/2, 3) = -dx/2.0D0*Ms(IMAX/2+1, 3)
  do i = IMAX/2, 1, -1
     Flux(i-1, 3) = Flux(i, 3) -dx*Ms(i, 3)
     if(Flux(i-1, 3) >0)then
        Flux(i-1, 3)=0.0D0
     end if
  end do
  Flux(IMAX/2+1, 3) = +dx/2.0D0*Ms(IMAX/2+1, 3)
  do i = IMAX/2+1, IMAX-1
     Flux(i+1, 3) = Flux(i,3) + dx * Ms(i+1, 3)
     if(Flux(i+1, 3) < 0)then
        Flux(i+1, 3) = 0.0D0
     end if
  end do

  write(6,*) IMAX/2 +1
  write(10, rec=1) REAL(VAR)
  write(10, rec=2) REAL(FLUX(1:IMAX, 1:JMAX))

  close(10)

  
  
  
  
  

end program analysis
