c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE cov_sta 

c     ******************************************************************
c
c     This routine computes the error covariance matrix for posterior x
c
c     ******************************************************************

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

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

c---- Local variables

      integer	i,j,k,kcomp,n,n1,n2,n3,n_tot,kco

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

      real*8 	pp(idim_n_o)	

      real*8	MAT_min,
     >		MAT_max,sum

      real*8	amiss

      parameter (amiss=-9.d0)
 
      integer	m,m_n

      character	string0*12

      logical	flag

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

      integer   values(8),nrhs,info




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

c---- Write

      write(*,*)
      write(*,*)
      write(*,*)
      write(*,*) '------ In cov_sta ------'
      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     Constraint matrix A (=E(2)) 
c
c     Note
c     ----
c     Formula (3.7.13, Wunsch, 1996) is used
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---- First n_li rows are statements for linear model constraints 

      if (iflag.ne.2) then

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

      endif


c---- Subsequent n_zon rows are statements for linear observational constraints

      if (iflag.ne.2) then

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

      endif 


c---- Last n_nl rows are statements of linearized non-linear equations
c
c      if (iflag.ne.2) then
c
c        n     = 0
c        do n1 = n_li+n_zon+1,n_li+n_zon+n_nl
c           n    = n + 1
c           do n2=1,n_o
c              AA(n1,n2) = FF(n,n2)
c           enddo
c        enddo
c
c      endif
c
c
c---- Case if iflag = 2
c
c      if (iflag.eq.2) then
c
c        do n1=1,n_nl
c        do n2=1,n_o
c           AA(n1,n2) = FF(n1,n2)
c        enddo
c        enddo
c
c      endif
c
c
c
c
c     ------------------------------------------
c     Covariance matrix for model errors (=R(2))
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^{-1} ...'


c---- First n_li rows are relative to errors in linear model constraints 

      if (iflag.ne.2) then

        do n=1,n_li
           C_AA_1(n) = Ca_mod_1(n)
        enddo

      endif


c---- Subsequent n_zon rows are statements for linear observational constraints

      if (iflag.ne.2) then

         n    = 0
         do n1=n_li+1,n_li+n_zon
            n          = n + 1
            C_AA_1(n1) = Ca_zon_1(n)
         enddo

      endif 


c---- Last n_nl rows are relative to statements of linearized non-linear equations
c
c      if (iflag.ne.2) then
c
c          n    = 0
c          do n1=n_li+n_zon+1,n_li+n_zon+n_nl
c             n          = n + 1
c             C_AA_1(n1) = Cf_1(n)
c          enddo
c 
c      endif
c 
c
c---- Case if iflag = 2
c
c      if (iflag.eq.2) then
c
c        do n=1,n_nl
c           C_AA_1(n) = Cf_1(n)
c        enddo
c  
c      endif
c
c
c---- Total number of rows of constraint matrix

c      n_tot = n_li + n_zon + n_nl

      n_tot = n_li + n_zon




c     --------------------------------------
c     Matrix to invert: Co^{-1} + A'Ca^{-1}A 
c     -------------------------------------- 

c---- Matrix A' x Ca^{-1}

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

      do n1=1,n_o
      do n2=1,n_tot
         MAT2(n1,n2) = AA(n2,n1)*C_AA_1(n2)
      enddo
      enddo


c---- Matrix product A' CA^{-1} A

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

      m   = 10
      m_n = n_o / m

      do i=1,m
         n1 = 1  + (i  -1)*m_n
         n2 = n1 + (m_n-1)
         if (i.eq.m) n2 = n_o
c        -------------------------------------------------
         MAT1(1:n_o,n1:n2) = matmul(MAT2(1:n_o  ,1:n_tot),
     >                              AA  (1:n_tot,n1:n2  ))
c        -------------------------------------------------
         write(*,*) '  Sub-multiplication # ',i
      enddo
      write(*,*)


c---- Sum Co^{-1} + A'CA^{-1}A

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

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




c     ----------------------------------------------
c     Inversion of Co^{-1} + A' C_A^{-1} A 
c
c     Note
c     ----
c     Inversion is done using Cholesky decomposition
c     ----------------------------------------------

c---- Testing if Co^{-1} + A' C_A^{-1} A is symmetric

      open (80,file='sym1')

      kco = 0

      do n1=1,n_o
      do n2=1,n_o
         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(*,*) '  Co^{-1} + AT C_A^{-1} A is strictly symmetric' 
      else
        write(*,*) '  Non-symmetric elements in sym1' 
      endif 
      write(*,*)


c---- Cholesky factorize Co^{-1} + A' CA^{-1}

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     ---------------------------------------
      call DPOTRF('L',n_o,MAT1,idim_n_o,info)
c     --------------------------------------- 

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


c---- Compute selected columns of error covariance matrix

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

c--   Initialization

      do n1=1,n_o
      do n2=1,n_o
         Coa(n1,n2) = amiss 
      enddo
      enddo 


c--   Solving a sequence of linear equations Ax = b

c++
      do n2=1,n_o     	! Loop on columns 
c++

c        Test for computing (co)variances

         flag    = .false.
         string0 = string(n2)

         if (string0(1:2).eq.'VV') then

           if (string0(7:8).eq.'08' .or.
     >         string0(7:8).eq.'36' .or.
     >         string0(7:8).eq.'42') flag = .true.

         endif


c        Compute (co)variances

         if (flag) then

           do n1=1,n_o
              MAT2(n1,1) = 0.d0
           enddo
           MAT2(n2,1) = 1.d0
           nrhs       = 1
c          ----------------------------------------------------------
           call DPOTRS('L',n_o,nrhs,MAT1,idim_n_o,MAT2,idim_n_o,info)
c          ----------------------------------------------------------

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

           do n1=1,n_o
              Coa(n1,n2) = MAT2(n1,1)
           enddo
           write(*,*) '  ',string0

         endif

c++
      enddo
c++




c     -------------------------------
c     Assignation of posterior errors
c     -------------------------------

c---- Interactive write

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


c---- Initialization

      kcomp = 0


c---- Errors in density

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (flag_r(i,j,k)) then
           kcomp = kcomp + 1
           if (Coa(kcomp,kcomp).ne.amiss) then 
             e_r(i,j,k) = sqrt(Coa(kcomp,kcomp)) 
           else
             e_r(i,j,k) =      Coa(kcomp,kcomp)
           endif
         endif

      enddo
      enddo
      enddo


c---- Errors in transport component U

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (flag_U(i,j,k)) then
           kcomp = kcomp + 1
           if (Coa(kcomp,kcomp).ne.amiss) then
             e_U(i,j,k) = sqrt(Coa(kcomp,kcomp))
           else
             e_U(i,j,k) =      Coa(kcomp,kcomp)
           endif
         endif

      enddo
      enddo
      enddo


c---- Errors in transport component V

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (flag_V(i,j,k)) then
           kcomp = kcomp + 1
           if (Coa(kcomp,kcomp).ne.amiss) then
             e_V(i,j,k) = sqrt(Coa(kcomp,kcomp)) 
           else
             e_V(i,j,k) =      Coa(kcomp,kcomp)
           endif 
         endif

      enddo
      enddo
      enddo


c---- Errors in transport component W

      do i=1  ,itot
      do j=1  ,jtot
      do k=1+1,ktot+1           ! W(i,j,1) is NOT in state vector

         if (flag_W(i,j,k)) then
           kcomp = kcomp + 1
           if (Coa(kcomp,kcomp).ne.amiss) then
             e_W(i,j,k) = sqrt(Coa(kcomp,kcomp))
           else
             e_W(i,j,k) =      Coa(kcomp,kcomp)
           endif
         endif

      enddo
      enddo
      enddo


c---- Errors in tracer 

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (flag_p(i,j,k)) then
           kcomp = kcomp + 1
           if (Coa(kcomp,kcomp).ne.amiss) then
             e_p(i,j,k) = sqrt(Coa(kcomp,kcomp)) 
           else
             e_p(i,j,k) =      Coa(kcomp,kcomp)
           endif
         endif

      enddo
      enddo
      enddo




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 cov_sta ------'
      write(*,*)
      write(*,*)




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

      return

      end

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

