c23456789012345678901234567890123456789012345678901234567890123456789012


      REAL*8 FUNCTION conse_r(opt,i,j,k) 

c     ******************************************************************
c
c     This routine calculates quantities in the density equation
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variables

      character	opt*1

      integer	i,j,k


c---- Local variables

      real*8    f_A,
     >          f_B,
     >          f_C,
     >          f_D,
     >          f_E,f_E_d,
     >          f_F,f_F_d




c     ----------------------------
c     Fluxes in governing equation
c     ----------------------------

c---- Advective fluxes

      f_A = -V(i  ,j  ,k  ) * (r(i  ,j  ,k  )+
     >                         r(i+1,j  ,k  )+
     >                         r(i+1,j  ,k+1)+
     >                         r(i  ,j  ,k+1))/4.d0

      f_B = +U(i+1,j  ,k  ) * (r(i+1,j  ,k  )+
     >                         r(i+1,j+1,k  )+
     >                         r(i+1,j+1,k+1)+
     >                         r(i+1,j  ,k+1))/4.d0

      f_C = +V(i  ,j+1,k  ) * (r(i  ,j+1,k  )+
     >                         r(i+1,j+1,k  )+
     >                         r(i+1,j+1,k+1)+
     >                         r(i  ,j+1,k+1))/4.d0

      f_D = -U(i  ,j  ,k  ) * (r(i  ,j  ,k  )+
     >                         r(i  ,j+1,k  )+
     >                         r(i  ,j+1,k+1)+
     >                         r(i  ,j  ,k+1))/4.d0

      f_E = -W(i  ,j  ,k  ) * (r(i  ,j  ,k  )+
     >                         r(i+1,j  ,k  )+
     >                         r(i+1,j+1,k  )+
     >                         r(i  ,j+1,k  ))/4.d0

      f_F = +W(i  ,j  ,k+1) * (r(i  ,j  ,k+1)+
     >                         r(i+1,j  ,k+1)+
     >                         r(i+1,j+1,k+1)+
     >                         r(i  ,j+1,k+1))/4.d0


c---- Diffusive fluxes

      if (e_den.lt.0.d0) then

        F_E = F_E + F_r(i,j,k  ) 
        F_F = F_F - F_r(i,j,k+1) 

      endif




c     --------------------
c     Function assignation
c     --------------------

      if (opt.eq.'R') then

        conse_r = f_A +
     >            f_B +
     >            f_C +
     >            f_D +
     >            f_E +
     >            f_F

      else if (opt.eq.'A') then

        conse_r = f_A

      else if (opt.eq.'B') then

        conse_r = f_B

      else if (opt.eq.'C') then

        conse_r = f_C

      else if (opt.eq.'D') then

        conse_r = f_D

      else if (opt.eq.'E') then

        conse_r = f_E

      else if (opt.eq.'F') then

        conse_r = f_F

      else
 
        write(*,*)
        write(*,*) 'BAD ARGUMENT IN conse_r'
        write(8,*)
        stop       'STOP'
  
      endif

     


c     --------------------
c     End of the function 
c     --------------------

      return

      end

c     ******************************************************************
