c23456789012345678901234567890123456789012345678901234567890123456789012


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

c     ******************************************************************
c
c     This routine quantities in the tracer transport 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,
     >		f_Pr

      real*8	S  ,
     >		vol,
     >		dp




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

c---- Advective fluxes

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

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

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

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

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

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


c---- Diffusive fluxes

      if (e_pos.lt.0.d0) then

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

      endif


c---- Production flux 
c     
c     Notes
c     -----
c     S   is in  [m2]
c     vol is in  [m3]
c     sv  is in  [Sv s/m3]
c     prod is in [mmol P/(m3 s)]

      S   = (r_0**2) * (sin(lat(j+1))-sin(lat(j))) * dx(i)
      vol = S * dz(k) 

      dp  = p0 - (p(i  ,j  ,k  )+
     >		  p(i+1,j  ,k  )+
     >		  p(i+1,j  ,k+1)+
     >		  p(i  ,j  ,k+1)+
     >		  p(i  ,j+1,k  )+
     >		  p(i+1,j+1,k  )+
     >		  p(i+1,j+1,k+1)+
     >		  p(i  ,j+1,k+1)) / 8.d0 

      f_Pr = - prod * vol * sv * dp 

      if (tag_ppp(1:3).eq.'o18') f_Pr = 0.d0


c      f_Pr = - prod * vol * sv  




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

      if (opt.eq.'R') then

        conse_p = f_A +
     >            f_B +
     >            f_C +
     >            f_D +
     >            f_E +
     >            f_F + f_Pr 

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

        conse_p = f_A

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

        conse_p = f_B

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

        conse_p = f_C

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

        conse_p = f_D

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

        conse_p = f_E

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

        conse_p = f_F

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

     


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

      return

      end

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