subroutine calc_flux(q, x, FF, Nx, Ng)

c------------------------------------------------------------
c This routine calculates the numerical fluxes
c at every cell boundary---FF(i) is the flux at 
c the cell boundary x_(i+1/2).
c------------------------------------------------------------

         implicit      none

         character*9   cdnm
         parameter   ( cdnm = 'calc_flux' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       neq
         parameter   ( neq = 1 )

         integer       Nx
         integer       Ng

         real*8        q      ( Nx ) 

         real*8        rq_iph 

         real*8        rqL_iph 
         real*8        rqR_iph 

         real*8        x  ( Nx )
         
         real*8        FF ( Nx ) 

         real*8        FF_iph 

         real*8        A

         real*8        lambda

         real*8        eta

         real*8        omega

         real*8        fL_iph 
         real*8        fR_iph 

         integer       i

c============================================================
c============================================================

         if ( ltrace ) then
            write(0,*) cdnm,':     Nx=', Nx
            write(0,*) cdnm,':     Ng=', Ng
         endif
c------------------------------------------------------------
c Note that we don't need to calculate the fluxes
c for boundaries between the ghost cells since
c the ghost cells are not updated using the equation 
c of motion.
c------------------------------------------------------------
         
         do i = Ng, Nx - Ng 
              if ( ltrace ) then
                 write(0,*) cdnm,': x(',i,') =',x(i) 
              endif
              
c------------------------------------------------------------
c Reconstrucction.
c Note that here x(i) is the position of the cell centres
c the output is the recostructed variables at x_iph.
c------------------------------------------------------------

              call recos_qL (q(i-1),  q(i),  q(i+1), 
     &                       x(i-1),  x(i),  x(i+1), 
     &                       rqL_iph )

              call recos_qR (q(i),  q(i+1),  q(i+2), 
     &                       x(i),  x(i+1),  x(i+2), 
     &                       rqR_iph )

c------------------------------------------------------------
c Calculation of the Characteristic Structure.
c------------------------------------------------------------

              rq_iph  = 0.5d0 * (rqR_iph  + rqL_iph) 

              call calc_A      (A,       rq_iph) 
              call calc_lambda (lambda,  A)
              call calc_eta    (eta,     A) 
              call calc_omega  (omega,   eta, 
     &                          rqL_iph, rqR_iph)

c------------------------------------------------------------
c Calculation of the Numerical Flux (FF)
c Note that FF is centred at x_iph
c------------------------------------------------------------

              call calc_f  (fR_iph, rqR_iph)
              call calc_f  (fL_iph, rqL_iph)
              call calc_FF (FF_iph, fL_iph,  fR_iph,
     &                      eta,    lambda,  omega)

              FF(i)  = FF_iph
         enddo 
         return
         end