subroutine hyperbolic

**********************************************************************
*                                                                    *
*       This subroutine is to solve the                              *
*       hyperbolic equation of the density n(i,j).                   *
*                                                                    *
**********************************************************************

        parameter (nz=126,nx=106)
        parameter (pi=3.141592654)
        real n,nt
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        COMMON /WX1/DZ,DX,Z0,X0,MOZ,MOX
        common /times/dt,t,tt
        common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx
        common /wx4/u1x0,u1z0,xk,zk,w,zki,uki
        common /wavelengths/lx,lz
        common /wx5/vin(-5:nz),ven(-5:nz)
        common /wx61/rvi(-5:nz),rix1(-5:nz),riz1(-5:nz),ri2(-5:nz)
        common /wx62/rve(-5:nz),rex1(-5:nz),rez1(-5:nz),re2(-5:nz)
        common /wx7/ti(-5:nz),te(-5:nz),sini,cosi
        common /wx8/epc,dperp
        common /wx10/produ(-5:nz),recom(-5:nz)
        common /wx32/dntd(-5:nz,-5:nx)
        common /wx33/dxz
        common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx)
        common /wx35/u1x0h(-5:nz),u1z0h(-5:nz),zkr(-5:nz)
        dimension nt(-5:nz,-5:nx)

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

*********************************************************************  
*                                                                   *
*       Calculate the amplitude of the gravity wave.                *
*       The vertical wavelength and amplitude of the gravity wave   *
*	vary with height.					    *
*                                                                   *
*********************************************************************  

        do 300, i=-5,nz
          do 300, j=-5,nx    

c        u1x0h=-u1z0h*zkr/xk           

            z=dz*real(i-1)-z0/2.0
            x=dx*real(j-1)-x0/2.0
            u1z(i,j)=u1z0h(i)*cos(xk*x+zkr(i)*z-w*t)
 300        u1x(i,j)=u1x0h(i)*cos(xk*x+zkr(i)*z-w*t)


**********************************************************************
*                                                                    *
*       Determine the time step 'dt'.                                *
*                                                                    *
**********************************************************************

        dt=0.0
 310    eps1=0.0
        do 312, i=2,moz-1
          do 312, j=1,mox
            vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*sini/((1.0+rvi(i)**2)*b)
     $  +sini*cosi*gg/((1.0+rvi(i)**2)*vin(i))
     $  -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
            vix1=abs(vix1/dx)
            viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*cosi/((1.0+rvi(i)**2)*b)
     $  -(rvi(i)*rvi(i)+sini**2)*gg/((1.0+rvi(i)**2)*vin(i))
     $  +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
            viz1=abs(viz1/dz)
 312        dt=max(dt,vix1,viz1)
        dt=0.35/dt
        if (dt.gt.10.0) then
          dt=10.0
          end if
        if (dt.lt.1.0) then
          dt=1.0
          end if

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

        do 322, i=-3,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=-3,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.                                            *
*                                                                    *
**********************************************************************

        do 350, j=1,mox
          n(2,j)=n(4,j)-(cosi/sini)*(n(3,j+1)-n(3,j-1))*dz/dx
          n(1,j)=n(3,j)-(cosi/sini)*(n(2,j+1)-n(2,j-1))*dz/dx
          n(moz-1,j)=n(moz-3,j)
     $  +(cosi/sini)*(n(moz-2,j+1)-n(moz-2,j-1))*dz/dx
          n(moz,j)=n(moz-2,j)
     $  +(cosi/sini)*(n(moz-1,j+1)-n(moz-1,j-1))*dz/dx
 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)

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

**********************************************************************
*                                                                    *
*       The functions used above.                                    *
*                                                                    *
**********************************************************************

        function c1(i,j)

        if (af(i,j).ge.0.0) then
          c1=min(r1(i,j+1),r0(i,j))
         else
          c1=min(r1(i,j),r0(i,j+1))
          end if
        return
        end

**********************************************************************

        function c2(i,j)

        if (ag(i,j).ge.0.0) then
          c2=min(r1(i+1,j),r0(i,j))
         else
          c2=min(r1(i,j),r0(i+1,j))
          end if
        return
        end

**********************************************************************

        function r1(i,j)

        parameter (nz=126,nx=106)
        common /wx32/dntd(-5:nz,-5:nx)
        common /wx33/dxz
        p(i,j)=max(0.0,ag(i-1,j))-min(0.0,ag(i,j))+
     $  max(0.0,af(i,j-1))-min(0.0,af(i,j))
        q(i,j)=(wmax(i,j)-dntd(i,j))*dxz
        if (p(i,j).gt.0.0) then
          r1=min(1.0,q(i,j)/p(i,j))
         else
          r1=0.0
          end if
        return
        end

**********************************************************************

        function r0(i,j)

        parameter (nz=126,nx=106)
        common /wx32/dntd(-5:nz,-5:nx)
        common /wx33/dxz
        p(i,j)=max(0.0,ag(i,j))-min(0.0,ag(i-1,j))+
     $  max(0.0,af(i,j))-min(0.0,af(i,j-1))
        q(i,j)=(dntd(i,j)-wmin(i,j))*dxz
        if (p(i,j).gt.0.0) then
          r0=min(1.0,q(i,j)/p(i,j))
         else
          r0=0.0
          end if
        return
        end

**********************************************************************

        function f(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx
        common /wx4/u1x0,u1z0,xk,zk,w,zki,uki
        common /wx5/vin(-5:nz),ven(-5:nz)
        common /wx61/rvi(-5:nz),rix1(-5:nz),riz1(-5:nz),ri2(-5:nz)
        common /wx62/rve(-5:nz),rex1(-5:nz),rez1(-5:nz),re2(-5:nz)
        common /wx7/ti(-5:nz),te(-5:nz),sini,cosi
        common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx)
        vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*sini/((1.0+rvi(i)*rvi(i))*b)
     $  +sini*cosi*gg/((1.0+rvi(i)*rvi(i))*vin(i))
     $  -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
        f=vix1*n(i,j)
        return
        end

**********************************************************************

        function g(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx
        common /wx4/u1x0,u1z0,xk,zk,w,zki,uki
        common /wx5/vin(-5:nz),ven(-5:nz)
        common /wx61/rvi(-5:nz),rix1(-5:nz),riz1(-5:nz),ri2(-5:nz)
        common /wx62/rve(-5:nz),rex1(-5:nz),rez1(-5:nz),re2(-5:nz)
        common /wx7/ti(-5:nz),te(-5:nz),sini,cosi
        common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx)
        viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*cosi/((1.0+rvi(i)*rvi(i))*b)
     $  -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i))
     $  +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
        g=viz1*n(i,j)
        return
        end

**********************************************************************

        function fl(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /times/dt,t,tt
        fl=(0.5*dt*(f(i,j+1)+f(i,j))-0.25*dx*(n(i,j+1)-n(i,j)))*dz
        return
        end

**********************************************************************

        function gl(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /times/dt,t,tt
        gl=(0.5*dt*(g(i+1,j)+g(i,j))-0.25*dz*(n(i+1,j)-n(i,j)))*dx
        return
        end

**********************************************************************

        function fh(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /times/dt,t,tt
        fh=(7.0*(f(i,j+1)+f(i,j))-(f(i,j+2)+f(i,j-1))
     $  )*dt*dz/12.0
        return
        end

**********************************************************************

        function gh(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx1/dz,dx,z0,x0,moz,mox
        common /times/dt,t,tt
        gh=(7.0*(g(i+1,j)+g(i,j))-(g(i+2,j)+g(i-1,j))
     $  )*dt*dx/12.0
        return
        end

**********************************************************************

        function ag(i,j)

        parameter (nz=126,nx=106)
        common /wx32/dntd(-5:nz,-5:nx)
        ag=gh(i,j)-gl(i,j)
        t1=ag*(dntd(i+1,j)-dntd(i,j))
        t2=ag*(dntd(i+2,j)-dntd(i+1,j))
        t3=ag*(dntd(i,j)-dntd(i-1,j))
        if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) ag=0.0
        return
        end

**********************************************************************

        function af(i,j)

        parameter (nz=126,nx=106)
        common /wx32/dntd(-5:nz,-5:nx)
        af=fh(i,j)-fl(i,j)
        t1=af*(dntd(i,j+1)-dntd(i,j))
        t2=af*(dntd(i,j+2)-dntd(i,j+1))
        t3=af*(dntd(i,j)-dntd(i,j-1))
        if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) af=0.0
        return
        end

**********************************************************************

        function wmax(i,j)

        parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx32/dntd(-5:nz,-5:nx)
        wa(i,j)=max(n(i,j),dntd(i,j))
        wmax=max(wa(i-1,j),wa(i,j),wa(i,j-1),wa(i,j+1),
     $  wa(i+1,j))
        return
        end

**********************************************************************

        function wmin(i,j)
         parameter (nz=126,nx=106)
        real n
        common n(-5:nz,-5:nx),phi(-5:nz,-5:nx)
        common /wx32/dntd(-5:nz,-5:nx)
        wb(i,j)=min(n(i,j),dntd(i,j))
        wmin=min(wb(i-1,j),wb(i,j),wb(i+1,j),wb(i,j-1),
     $  wb(i,j+1))
        return
        end