C
C**CASCO 
C** the skeleton of the CASCO linear inverse model
C
C casco4a: Nov 2002.  Add elevation observations to adjoint
C casco3b: Aug 1999.  w2 and both gradient and conjugate gradient descent 
C
CGH CASCO assumes a cold start 23 Oct 1998
C
C Nov 27, 1998 casco2 replaces casco1 
C         changed iteration loop
C         gutted moody_start and moody_stop 
C
C
      PROGRAM MAIN
      include 'CASCO.DIM'
      include 'wind.h'
      integer itermax
      integer iter 
      integer method_descent
      logical converged

C  functions
      logical convergence 
      integer max_iterations, descent_code

      iter = 0
      call casco_info(6)
      call casco_input
      itermax = max_iterations() 
      write(*,*)'casco2: maximum number of iterations ', itermax
      method_descent = descent_code()
      write(*,*)'Method of descent ', method_descent
      call startup
      call get_friction(iter)
      call station_calc
      call set_times
C
C        call getdata_moody to read the file(s) containing the
C             unexplained observations and fill the appropriate arrays.
C        call measure_func to compute the measurements functionals
C        call get_errvar to get the variances of the velocity errors
C        call get_Whro to Build the BC smoothing with functions 
C        call initial_conditions to set the initial conditions
C                                (identically zero as of Arpil 1999) 
C        call statistics to compute the initial statistics  

      call getdata_moody
      call measure_func
      call get_errvar
      call get_Wrho
      call initial_conditions
      call statistics(iter)
      if (Nlags.gt.0 )call timeseries_offset 

      converged = .false. 

C repeat until done 
100   continue 
         iter = iter + 1
	 write(*,*)'moody_saco iteration: ',iter

C the output of moody is DomegaDrho (unscaled)
         call moody_run(iter)

c compute the trial bcs for the saco run. 
         if(method_descent.eq.1) then 
	    call trial_bcs(iter)
	 elseif(method_descent.eq.2) then
	    if (nlags.gt.0)then 
 	     call trial_bcs_cgm_W(iter)
            else
	     call trial_bcs_cgm(iter)
	    endif
	 else
	    write(*,*)'Unknown method of descent: ',method_descent
	    stop 1111
	 endif

C the output of saco is the new error field
	 call saco_run(iter)

C compute the correct step size, adjust the errors and the bcs. 
C compute the cost function and other statistics 
         call adjust_soln(iter)
         call statistics(iter)

C test for convergence
         converged = convergence(iter)
	 write(*,*)'CONVERGENCE: ', iter,converged

      if((iter.lt.itermax).and.(.not.converged)) goto 100 
101   continue

      write(*,*)'CASCO. finished iteration loop '

C                       write the final bcs to a file (*.cbc)
C                   and write the final error to a file
      call write_cbc
c
c in lagged convolution version output convolution coefficients
c
      if (nlags.gt.0)call write_lcc
	
      call write_error

C make a final run of saco to get best time series and check final error
      call dbc_best_to_dbc
      iter = -1 
      call saco_run(iter)
      call error_to_best
      
      call write_error
      stop
      end

