c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE Coa_new

c     ******************************************************************
c
c     This routine sets the matrix CoA 
c
c     ******************************************************************

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

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

c---- Local variables

      integer	n,n1,n2,n3,n_tot,kco

      real*8	MAT_min,MAT_max,sum
     
      real*8	pp(idim_n_o)

      real*8	AA     (idim_n_li+(jtot+1),idim_n_o         ),
     >		C_AA   (idim_n_li+(jtot+1)                  )

      real*8	MATT(idim_n_li+(jtot+1),idim_n_li+(jtot+1))

      integer	i,j,k,m,m_n	

      character date*8,time*10,zone*5

      integer   values(8),info




c     --------------------
c     Interactive writing
c     --------------------

c---- Write

      write(*,*)
      write(*,*)
      write(*,*) '------ In CoA_new ------'
      write(*,*)


c---- Time

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) 'Enter time=',time(1:2)//':'//
     >                         time(3:4)//':'//
     >                         time(5:6)
      write(*,*)




c     ----------------------------
c     Matrix of linear constraints
c     ----------------------------

c---- Interactive writing

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get A ...'


c---- Total number of linear constraints

      n_tot = n_li + n_zon


c---- Fill A(,) with coefficients of model constraints

      do n1=1,n_li
      do n2=1,n_o
         AA(n1,n2) = A_mod(n1,n2)
      enddo
      enddo 


c---- Fill A(,) with coefficients of observational constraints

      if (iflag.ne.2 .and. iflag.ne.3) then

        n = 0
        do n1=n_li+1,n_tot
           n = n + 1
           do n2=1,n_o
              AA(n1,n2) = A_zon(n,n2)
           enddo
        enddo

      endif




c     ---------------------------------------------- 
c     Error covariance matrix for linear constraints
c     ----------------------------------------------

c---- Interactive writing

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get CA ...'


c---- Fill CA(,) with errors in model constraints

      do n=1,n_li
         C_AA(n) = 1. / Ca_mod_1(n)
      enddo


c---- Fill CA(,) with errors in observational constraints 

      if (iflag.ne.2 .and. iflag.ne.3) then

         n = 0
         do n1=n_li+1,n_tot
            n = n +1
            C_AA(n1) = 1. / Ca_zon_1(n)
         enddo

      endif
 



c     --------------------------------
c     Matrix to invert: A Co A'  +  CA 
c     -------------------------------- 

c---- Matrix Co A'

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get Co x AT ...' 

      do n1=1,n_o
      do n2=1,n_tot
         Coa(n1,n2) = AA(n2,n1) / Co_1(n1)
      enddo
      enddo


c---- Matrix product A x (Co A')

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get A x Co x AT ...'

      MAT1(1:n_tot,1:n_tot) = matmul(AA (1:n_tot,1:n_o  ),
     >				     Coa(1:n_o  ,1:n_tot))


c---- Sum A Co A' + CA

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get A x Co x AT  +  CA ...'

      do n=1,n_tot
         MAT1(n,n) = MAT1(n,n) + C_AA(n)
      enddo




c     ---------------------------------
c     Cholesky factorize A Co A'  +  CA 
c     ---------------------------------

c---- Number & fraction of non-zero elements in A Co A'  +  CA 

      kco = 0

      do n1=1,n_tot
      do n2=1,n_tot
         if (MAT1(n1,n2).ne.0.d0) kco = kco + 1
      enddo
      enddo

      write(*,*  )
      write(*,116) kco 
      write(*,117) dfloat(kco) / dfloat(n_tot*n_tot)

116   format('   Number   of non-zero elements = ',i7   )
117   format('   Fraction of non-zero elements = ',e10.3)


c---- Minimum and maximum values in matrix 

      MAT_min = 1.d20
      MAT_max =-1.d20

      do n1=1,n_tot
      do n2=1,n_tot
         MAT_min = min(MAT_min,MAT1(n1,n2))
         MAT_max = max(MAT_max,MAT1(n1,n2))
      enddo
      enddo

      write(*,*  )
      write(*,118) MAT_min
      write(*,119) MAT_max

118   format('   Minimum = ',e10.3)
119   format('   Maximum = ',e10.3)


c---- Test symmetry of A Co A'  +  CA 

      open (80,file='sym0')

      kco = 0

      do n1=1,n_tot
      do n2=1,n_tot
         if (MAT1(n1,n2).ne.MAT1(n2,n1)) then
            kco = kco + 1
            write(80,*) kco,MAT1(n1,n2)-MAT1(n2,n1)
         endif
      enddo
      enddo 

      close (80) 

      write(*,*)
      if (kco.eq.0) then
        write(*,*) '  A Co AT  +  CA is strictly symmetric' 
      else
        write(*,*) '  Non-symmetric elements in sym0'
      endif 
      write(*,*)


c---- Cholesky factorize A Co A' + CA 

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//
     >           '  Get Cholesky factor L ...'


c     Use ESSL routine to Cholesky factorize

c     -----------------------------------------
      call DPOTRF('L',n_tot,MAT1,idim_n_o,info)
c     -----------------------------------------

      if (info.ne.0) then
         write(*,*)
         write(*,*) 'Factorization did not succeed'
         write(*,*)
         stop 'OK'
      endif


c---- Invert Cholesky factor L

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get L^{-1} ...'


c     Put zero above diagonal of L 
c         p()  along diagonal of L

      do n1=1,n_tot
      do n2=1,n_tot
         if (n1.lt.n2) then
            MAT1(n1,n2) = 0.d0
         endif
      enddo
      enddo


c     Use ESSL routine to invert L

c     ---------------------------------------------
      call DTRTRI('L','N',n_tot,MAT1,idim_n_o,info)
c     ---------------------------------------------

      if (info.ne.0) then
         write(*,*)
         write(*,*) 'Matrix is singular and cannot be inverted'
         write(*,*) 
         stop 'OK'
      endif




c     --------------------------------------------
c     Matrix Coa = Co - Co A' (ACoA'+CA)^{-1} A Co
c                = Co - Co (L^{-1} A)' L^{-1} A Co
c     --------------------------------------------

c---- Matrix L^{-1} x A

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get L^{-1} x A ...'

c--

      m   = 10
      m_n = n_o / m

c++
      do i=1,m
c++
                  
         n1 = 1  + (i  -1)*m_n
         n2 = n1 + (m_n-1)

         if (i.eq.m) n2 = n_o

         MAT2(1:n_tot,n1:n2) = matmul(MAT1(1:n_tot,1:n_tot),
     >                                AA  (1:n_tot,n1:n2  ))
     
         if (mod(i,100).eq.0) then
c          -----------------------------------------
           call date_and_time(date,time,zone,values)
c          -----------------------------------------
           write(*,*) '  ',time(1:2)//':'//
     >                     time(3:4)//':'//
     >                     time(5:6)//'  Sub-multiplication # ',i 
         endif

c++
      enddo
c++


c---- Transpose of L^{-1} x A

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get transpose of L^{-1} x A ...'

c--

      do n1=1,n_o
      do n2=1,n_tot
         MAT1(n1,n2) = MAT2(n2,n1)
      enddo
      enddo


c---- Product (L^{-1} A)'  x  L^{-1} x A 

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//'  Get (L^{-1} A)T  x  L^{-1} A ...'

c--

      m   = 10
      m_n = n_o / m

c++
      do i=1,m
c++

         n1 = 1  + (i  -1)*m_n
         n2 = n1 + (m_n-1)

         if (i.eq.m) n2 = n_o

         Coa(1:n_o,n1:n2) = matmul(MAT1(1:n_o  ,1 :n_tot),
     >				   MAT2(1:n_tot,n1:n2  ))

         if (mod(i,100).eq.0) then
c          -----------------------------------------
           call date_and_time(date,time,zone,values)
c          -----------------------------------------
           write(*,*) '  ',time(1:2)//':'//
     >                     time(3:4)//':'//
     >                     time(5:6)//'  Sub-multiplication # ',i
         endif

c++
      enddo
c++


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

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//
     >		 '  Get Co   x   (L^{-1} A)T x L^{-1} A   x   Co ...'

c--

      do n1=1,n_o
      do n2=1,n_o
         Coa(n1,n2) = Coa(n1,n2) / (Co_1(n1)*Co_1(n2))
      enddo
      enddo


c---- Matrix Coa = Co - CoA' (ACoA'+CA)^{-1} ACo

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) time(1:2)//':'//
     >           time(3:4)//':'//
     >           time(5:6)//
     >	         '  Get Coa = Co  -  Co x (L^{-1}A)T (L^{-1}A) x ACo ...'

      do n1=1,n_o
      do n2=1,n_o
         if (n1.eq.n2) then
           Coa(n1,n2) = 1./Co_1(n1) - Coa(n1,n2)
         else
           Coa(n1,n2) =             - Coa(n1,n2)
         endif
      enddo
      enddo




c     ----------------------------
c     Diagnostics of the inversion
c     ----------------------------

c---- GOTO statement

      goto 111 


c---- Original matrix to invert: A' x CA^{-1}

      write(*,*)
      write(*,*) 'Check: Get AT  x  CA^{-1} ...'

      do n1=1,n_o
      do n2=1,n_tot
         MAT1(n1,n2) = AA(n2,n1) / C_AA(n2)
      enddo
      enddo

     
c---- Original matrix to invert: A' CA^{-1}  x  A

      write(*,*) 'Check: Get AT CA^{-1}  x  A ...'

      MAT2(1:n_o,1:n_o) = matmul(MAT1(1:n_o  ,1:n_tot),
     >				 AA  (1:n_tot,1:n_o  ))


c---- Original matrix to invert: Co^{-1}  +  A' CA^{-1} A 
    
      write(*,*) 'Check: Get Co^{-1}  +  AT CA^{-1} A ...'

      do n=1,n_o
         MAT2(n,n) = MAT2(n,n) + Co_1(n)
      enddo


c---- Computing matrix product

      write(*,*) 'Check: Get I_est ...' 

      MAT1(1:n_o,1:n_o) = matmul(Coa (1:n_o,1:n_o),
     >			         MAT2(1:n_o,1:n_o))


c---- Computing and writing the sum of diagonal elements of 
c     matrix product

      sum = 0.d0
      do n1=1,n_o
         sum = sum + MAT1(n1,n1)
      enddo

      write(*,*) 'Check:'
      write(*,*) '   Sum  of diagonal elements of I_est=',sum
      write(*,*) '   Size of problem                   =',n_o


c---- Writing diagonal and off-diagonal elements of matrix product

      write(*,*) 'Check: Write I_est in on_diagonal & of_diagonal ...' 

      open (50,file='on_diagonal.'//tag)

      n3 = 0 

      do n1=1,n_o
      do n2=1,n_o
         if (n1.eq.n2) then
           n3 = n3 + 1
           write(50,*) n3,MAT1(n1,n2)
         endif
      enddo
      enddo

      close (50)

c--
c     
c      open (50,file='of_diagonal.'//tag)
c 
c      n3 = 0 
c
c      do n1=1,n_o
c      do n2=1,n_o
c         if (n1.ne.n2) then
c           n3 = n3 + 1
c           write(50,*) n3,MAT1(n1,n2)
c         endif
c      enddo
c      enddo

      stop 'OK OLIVIER'




c     --------------------
c     Interactive writing
c     --------------------

c---- Time

c     -----------------------------------------
111   call date_and_time(date,time,zone,values)
c     -----------------------------------------


c---- Write

      write(*,*)
      write(*,*) 'Exit time=',time(1:2)//':'//
     >                        time(3:4)//':'//
     >                        time(5:6)

      write(*,*)
      write(*,*) '------ Out CoA_new ------'
      write(*,*)




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

      return

      end

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

