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

program analysis

  use prm , only : IMAX, JMAX, KMAX, dx
  use com_var, only : get_frame, x, zk, zk_half
  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)
  real(8)   :: Uvec(1:IMAX, JMAX, 0:KMAX)
  real(8)   :: Usfc(1:IMAX, JMAX)
  real(8)   :: Uscle(1:IMAX, JMAX)


  integer :: i, j , k

  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


  Uvec(0, :, :)    = 0.0D0
  Uvec(IMAX, :, :) = 0.0D0

  do j = 1, JMAX
     do i = 1,IMAX-1
        Uscle(i, j) = - 2.0D0 * (rho *grav)** Npw     &
        &             * ((var(i+1, j) - var(i, j))/dx)**Npw    &
        &             * ((var(i+1, j) +var(i, j))/2.0D0)**(Npw+1.0D0)

     ! $B:G2<AX(B
        Uvec(i, j, 1) =  Uscle(i, j) * A_con &
          &       * (zk_half(1) - zk_half(0))/2.0D0 &
          &       * ((1 - zk(1))** Npw  + 1.0D0) /2.0D0

     ! $B:G2<AX$+$i:G>eAX(B
        do k = 2, KMAX
           Uvec(i, j, k) = Uvec(i, j, k-1)                    &
                &  + Uscle(i, j)  * A_con * (zk(k) - zk(k-1))  &
                &  * ((1 - zk(k-1))**Npw + (1 -zk(k))**Npw)/2.0D0
        end do

     ! $BI=LLB.EY(B( 2003/11/26 $BDI5-(B)
        Usfc(i, j) = Uvec(i, j, KMAX)                        &
             & + Uscle(i, j) * A_con                           &
             & * (zk_half(KMAX) - zk_half(KMAX-1))/2.0D0       &
             & * ((1 - zk(KMAX))**Npw)/2.0D0
        
        
     end do
  end do



!  do j = 1, JMAX
!     do i =1, IMAX-1


!        Usfc(i, j) = - 2.0D0 * (rho *grav)** Npw                           &
!             &     * ((var(i+1, j) - var(i, j))/dx)**Npw             &
!             &     * ((var(i+1, j) + var(i, j))/2.0D0)**(Npw+1.0D0)  &
!             &     *  A_con /(Npw + 1.0D0)                
                ! $B3FAX$N7W;;(B
!        do k = 1, KMAX
!           Uvec(i, j, k) = Usfc(i, j)                         &
!        &             * (1 - (1 - zk(k))**(Npw+1.0D0)) 
!        end do

        
!     end do
!  end do





  write(6,*) IMAX/2 +1
  write(10, rec=1) REAL(VAR)
  write(10, rec=2) REAL(FLUX(1:IMAX, 1:JMAX))
  do i =1, KMAX
     write(10, rec=i+2) REAL(uvec(1:IMAX, 1:JMAX, i)*s_a)
  end do

  close(10)

  
  
  
  
  

end program analysis
