subroutine hyperbolic

**********************************************************************
*                                                                    *
*       This subroutine is to solve the                              *
*       hyperbolic equation of the density n(i,j).                   *
*                                                                    *
*       Version 0.5 - Last Modified 7/10/04                          *
*                                                                    *
**********************************************************************

        INCLUDE "fpvar.inc"

        dimension nt(ST:NZ,ST:NX)

        re=0.05
        eps0=1.0
        epmax=1.0e18
        dxz=dx*dz

**********************************************************************
*                                                                    *
*       Determine the time step 'dt' by first finding the            *
*       maximum velocity in any given cell in either the x-          *
*       or z-direction.                                              *
*                                                                    *
**********************************************************************

        DT=0.0
 310    EPS1=0.0
        DO 312, I=3,MOZ-2
           DO 312, J=1,MOX
              VIX1=ABS(F(I,J)/(N(I,J)*DX))
              VIZ1=ABS(G(I,J)/(N(I,J)*DZ))
 312          DT=MAX(DT,VIX1,VIZ1)
        DT=0.35/DT
        IF (DT.GT.10.) DT=10.
        IF (DT.LT.1.)  DT=1.

**********************************************************************
*                                                                    *
*       Calculate the time advanced solution for lower order flux.   *
*                                                                    *
**********************************************************************

        DO 322, I=ST+2,NZ-2
           AVEN=0.0

	   DO 320, J=1,MOX-1   
 320          AVEN=AVEN+N(I,J)
 	   AVEN=AVEN/REAL(MOX-1)    

           do 322, j=ST+2,nx-2
              dntd(i,j)=n(i,j)+dt*(produ(i)-recom(i)*n(i,j))
     $  -(FL(I,J)-FL(I,J-1)+GL(I,J)-GL(I-1,J))/DXZ
     $  -dt*((cosi*xk-sini*zkr(i))**2)*(n(i,j)-aven)
     $  *ti(i)/(1.0+rvi(i)**2)
 322         continue

**********************************************************************
*                                                                    *
*       Calculate the time advanced solution for higher order flux   *
*       in the inner region (i=3:moz-2, j=1:mox).                    *
*                                                                    *
**********************************************************************

        do 330, i=3,moz-2
           do 330, j=1,mox
              nt(i,j)=dntd(i,j)-(c1(i,j)*af(i,j)-c1(i,j-1)*af(i,j-1)
     $  +c2(i,j)*ag(i,j)-c2(i-1,j)*ag(i-1,j))/dxz
 330          eps1=max(eps1,abs(nt(i,j)-n(i,j)))

**********************************************************************
*                                                                    *
*       Put the values of nt(i,j) to n(i,j).                         *
*                                                                    *
**********************************************************************

      do 340, i=3,moz-2
         do 340, j=1,mox
 340        n(i,j)=nt(i,j)

**********************************************************************
*                                                                    *
*       Calculate the values of n(2,j), n(1,j), n(moz-1,j), and      *
*       n(moz,j) for j=1,mox under suitable boundary conditions.     *
*       These values are not calculated in the calculations for      *
*       the inner region.                                            *
*                                                                    *
**********************************************************************

c        do 350, j=1,mox
c           n(2,j)=n(4,j)-(cosi/sini)*(n(3,j+1)-n(3,j-1))*dz/dx
c           n(1,j)=n(3,j)-(cosi/sini)*(n(2,j+1)-n(2,j-1))*dz/dx
c           n(moz-1,j)=n(moz-3,j)
c     $  +(cosi/sini)*(n(moz-2,j+1)-n(moz-2,j-1))*dz/dx
c           n(moz,j)=n(moz-2,j)
c     $  +(cosi/sini)*(n(moz-1,j+1)-n(moz-1,j-1))*dz/dx
c 350       continue

**********************************************************************
*                                                                    *
*       Restrict sudden change of n(i,j) caused by numerical error.  *
*                                                                    *
**********************************************************************

        do 360, i=1,moz
           do 360, j=1,mox
              if (n(i,j).lt.0.0) then
                 n(i,j)=0.25*(n(i+1,j)+n(i-1,j)+n(i,j+1)+n(i,j-1))
                 end if
 360       continue

**********************************************************************
*                                                                    *
*       Periodic boundary condition in the horizontal direction.     *
*                                                                    *
**********************************************************************

 370    do 372, j=-5,0
           do 372, i=1,moz
 372          n(i,j)=n(i,mox+j-1)
        do 374, j=mox+1,nx
           do 374, i=1,moz
 374          n(i,j)=n(i,j-mox+1)

        open (status='old',access='append',unit=8,file='trouble.dat')
        write(8,*)'end of hyperbolic'
        close(8)
        return
        end