c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE cost_fu 

c     ******************************************************************
c
c     This routine computes the different components of the cost function
c
c     ******************************************************************

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

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

c---- Local variables

      integer	n,n1,n2

      real*8	vec(idim_n_o),
     >		Ax (idim_n_o)




c     ------------------------
c     Component relative to xo 
c     ------------------------

c---- Product (x-xo)' Co^{-1}
 
      do n=1,n_o
         vec(n) = (x(n)-xo(n))*Co_1(n)
      enddo 


c---- Product (x-xo)' Co^{-1} (x-xo)

      J_o = 0.d0
      do n=1,n_o
         J_o = J_o + vec(n)*(x(n)-xo(n))
      enddo




c     ----------------------------------------------
c     Component relative to linear model constraints 
c     ----------------------------------------------

c---- Product of A_mod and x

      do n1=1,n_li

         Ax(n1) = 0.d0
         do n2=1,n_o
            Ax(n1) = 
     >	    Ax(n1) + A_mod(n1,n2)*x(n2)
         enddo

      enddo


c---- Product (A_mod x)' Ca_mod^{-1}

      do n=1,n_li
         vec(n) = Ax(n)*Ca_mod_1(n) 
      enddo


c---- Product (A_mod x)' Ca_mod^{-1} (A_mod x)

      J_li = 0.d0
      do n = 1,n_li
         J_li = J_li + vec(n)*Ax(n)
      enddo




c     ------------------------------------------------------
c     Component relative to linear observational constraints
c     ------------------------------------------------------

c---- Product of A and x

      do n1=1,n_zon
         Ax(n1) = 0.d0
         do n2=1,n_o
            Ax(n1) = Ax(n1) + A_zon(n1,n2)*x(n2)
         enddo
      enddo 


c---- Difference Ax - b

      do n=1,n_zon
         Ax(n) = Ax(n) - b_zon(n)
      enddo 


c---- Product (Ax-b)' CA^{-1} 

      do n=1,n_zon
         vec(n) = Ax(n)*Ca_zon_1(n)
      enddo


c---- Product (Ax-b)' CA^{-1} (Ax-b)

      J_zon = 0.d0
      do n  = 1,n_zon
         J_zon = J_zon + vec(n)*Ax(n)
      enddo




c     -------------------------------------------------
c     Component relative to nonlinear model constraints 
c     -------------------------------------------------

c---- Product f' Cf^{-1}

      do n=1,n_nl
         vec(n) = f(n)*Cf_1(n) 
      enddo


c---- Product f' Cf^{-1} f

      J_nl = 0.d0
      do n = 1,n_nl
         J_nl = J_nl + vec(n)*f(n)
      enddo
 
 


c     -------------------
c     End of the routine
c     -------------------

      return

      end

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

