c23456789012345678901234567890123456789012345678901234567890123456789012


      PROGRAM main 

c     ******************************************************************
c
c     This program combines (paleo)observations with model constraints 
c
c     Method: Method of total inversion of Tarantola & Valette (1982)
c     ------                               Mercier             (1986)
c
c     ******************************************************************


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

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

c---- Local variables

      integer	i,j,k,kmax,kcomp,kco1,kco2

      integer	n,n1,n2

      real*8	J_tot 




c     -----------------------------
c     Reading inputs in START file 
c     -----------------------------

c---- Extension of START file

      write(*,*)
      write(*,*) '-------------------------------------'
      write(*,*) 'Enter the extension of the START file'
      write(*,*) '-------------------------------------'
      write(*,*)

      read (*,*)  tag

      write(*,*)
      write(*,*) 'Name of START file: ',tag
      write(*,*)


c---- Opening START file

      open (08,file='start.'//tag,status='old')


c---- Readings in START file 

      read(8,*) e_r0            ! Error in a priori r(,,) 
      read(8,*) e_U0 		! Error in a priori U(,,) 
      read(8,*) e_V0		! Error in a priori V(,,) 
      read(8,*) e_W0            ! Error in a priori W(,,) 
      read(8,*) e_p0		! Error in a priori p(,,)
c--
      read(8,*) e_zon0		! Errors in a prior volume transports
      read(8,*) e_wbc0	
      read(8,*) e_mow0
      read(8,*) e_vem0
      read(8,*) e_hun0
      read(8,*) e_rac0
      read(8,*) e_eqc0
      read(8,*) e_vfz0
c--
      read(8,*) e_mas		! Error in a priori mass conservation 
      read(8,*) e_twu		!                   thermal wind equation (u)
      read(8,*) e_twv           !                   thermal wind equation (v)
      read(8,*) e_vor		!                   linear vorticity balance
      read(8,*) e_den           !                   density equation
      read(8,*) e_pos           !                   tracer  equation
c--
      read(8,*) b_zon(3)	! Error in a priori integrated meridional transports
      read(8,*) b_zon(2)
      read(8,*) b_zon(1)
c--
      read(8,*) b_wbc0_up,b_wbc0_do	! A prior volume transports
      read(8,*) b_mow0
      read(8,*) b_vem0
      read(8,*) b_hun0
      read(8,*) b_rac0
      read(8,*) b_eqc0
      read(8,*) b_vfz0
c--
      read(8,*) lnm		! A priori level of no motion
c--
c
c     Notes
c     -----
c     iflag = 0 means only r(WOCE)             are used
c     iflag = 1 means      r(WOCE) and p(WOCE) are used
c     iflag = 2 means (U,V,W) are read from previous inversion 
c                     linear constraints (t. wind) are     used 
c                     non-linear constraints for r are     used 
c                     non-linear constraints for p are not used 
c     iflag = 3 means (U,V,W) are read from previous inversion
c                     linear constraints           are not used
c                     non-linear constraints for r are not used
c                     non-linear constraints for p are     used

      read(8,*) iflag

      if (iflag.eq.2) then
        read(8,*) tag_UVW
        read(8,*) tag_rrr
      endif

      if (iflag.eq.3) then
        read(8,*) tag_UVW
        read(8,*) tag_ppp
      endif


c---- Closing START file

      close (8)




c     -------------------------------
c     Open files, model grid, & flags
c     -------------------------------

c---- Opening files

c     ------------
      call opening
c     ------------


C---- Model grid

c     -----------------------
      call modgrid(1,'AT002')
c     -----------------------


c---- Flags

c     ---------------------
      call flag_vs('AT002')
c     ---------------------




c     -------------------------------------------
c     Computing vectors/matrices before iteration 
c     -------------------------------------------

c---- String identifier for state variables

c     ------------
      call x_ident
c     ------------


c---- Path of WBC

c     Resolution of 3 degrees

c     ---------------------------------------
c      call pathwbc(20,34    ,2,-2,1     ,kco2)
c      call pathwbc(21,jtot-1,3,-1,kco2+1,n_wbc)
c     ---------------------------------------

c     Resolution of 2 degrees

c     -----------------------------------------
      call pathwbc(30,51    ,2,-2,1     ,kco1 )
      call pathwbc(33,jtot-3,3,-2,kco1  ,kco2 )
      call pathwbc(31,jtot-1,3,-2,kco2  ,n_wbc)
c     -----------------------------------------


c---- Coefficient matrices

c     ---------------------
      call A_modxx
      call A_zonxx('AT002')
c     ---------------------
   

c---- Covariance matrices

c     ---------------------
      call Co_1xxx('AT002')
      call Ca_modx
      call Ca_zonx
      call Cfxxxxx
c     ---------------------


c---- Interactive writing of dimensions

      write(*,*) 'n_o  =',n_o
      write(*,*) 'n_li =',n_li
      write(*,*) 'n_zon=',n_zon
      write(*,*) 'n_nl =',n_nl
      write(*,*)
      write(*,*) 'n_r=',n_r
      write(*,*) 'n_U=',n_U
      write(*,*) 'n_V=',n_V
      write(*,*) 'n_W=',n_W
      write(*,*) 'n_p=',n_p
      write(*,*)
      write(*,*) 'n_mas=',n_mas
      write(*,*) 'n_twu=',n_twu
      write(*,*) 'n_twv=',n_twv
      write(*,*) 'n_vor=',n_vor
      write(*,*)
      write(*,*) 'n_wbc=',n_wbc
      write(*,*)


c---- Covariance matrix C_oA^{-1}

      if (iflag.eq.3) then

        do n1=1,n_o
        do n2=1,n_o
           Coa(n1,n2) = 0.d0
           if (n1.eq.n2) Coa(n1,n2) = 1.d0 / Co_1(n1)
        enddo
        enddo

      else 

c       ------------
        call Coa_new
c       ------------
 
      endif


c---- Assignation of values to state variables

      if (iflag.eq.2) then

c       -----------------------
        call assignx('var_den')
c       -----------------------

      else if (iflag.eq.3) then

c       -----------------------
        call assignx('var_pos')
c       -----------------------

      else 

c       ---------------------
        call x_prior('AT002')
c       ---------------------

      endif


c---- Write model set-up, prior xo, and error in prior xo 

c     ---------------
      call writing(0)
c     ---------------

      stop 'Look at prior files'


c---- Write posterior x, error in posterior x, model residuals, and data misfits

c     --------------------
      call residua('dum')
      call writing(1)
c     --------------------

c      stop 'Look at residual files'


c---- Prior state 

c     ------------
      call xo_fill
c     ------------




c     -------------------
c     Iterative procedure
c     -------------------

c---- State used to initialize inversion 

      do n=1,n_o
         x(n) = xo(n)
      enddo
c     ------------
c      call noising
c     ------------


c---- General loop

      kcomp = 0 

c     #####################
      do while (kcomp.le.2) 
c     #####################


c---- Iteration number

      kcomp = kcomp + 1


c---- Computing vector of non-linear functions

c     -------------
      call fffffff
c     -------------


c---- Computing Jacobian matrix of non-linear functions

c     -------------
      call FFxxxxx
c     ------------- 


c---- Computing and writing objective function

c     ------------
      call cost_fu
c     ------------ 

c--

      J_tot = J_o + J_zon + J_li + J_nl

      write(70,'(i3,4e10.3)') kcomp,J_o           
     >                             ,J_zon + J_li
     >                             ,J_nl          
     >                             ,J_tot
     
      write(*,*)
      write(*,*)
      write(*,200) kcomp,J_tot
200   format('Iteration # ',i2,' - J=',e20.10)
      write(*,*)


c---- Computing vector beta

c     ------------
      call betaxxx
c     ------------


c---- Computing vector alpha

c     -------------------
      call alphaxx(kcomp)
c     -------------------


c---- Computing new state vector x and assigning new value to each state variable 

c     ------------
      call new_xxx
      call new_sta
c     ------------


c---- End of general loop

c     ######
      enddo
c     ######




c     -----------------------------------------
c     - Writing results
c     - Error covariance matrix for posterior x
c     - Compute residuals
c     - Write results
c     -----------------------------------------

c---- Write results

c     ---------------
      call writing(1)
c     ---------------


c---- Covariance matrix for posterior

c     ----------------------------
      if (iflag.ne.2 .and.
     >	  iflag.ne.3) call cov_sta
c     ----------------------------


c---- Residuals

c     -------------------
      call residua('nor')
c     -------------------


c---- Write results
      
c     ---------------
      call writing(1)
c     ---------------




c     ---------------------------
c     Advertising of normal exit
c     ---------------------------

      kmax=50
      do k=50,1,-1
         if (tag(k:k).eq.' ') kmax=k-1
         if (kmax    .eq.k  ) goto 2 
      enddo
2     continue

      write(*,*)
      write(*,*) 'Exit OK - Extension of result files = ',tag(1:kmax) 
      write(*,*)




c     --------------------
c     End of main program
c     --------------------
 
      end

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

