C** moody_run
C
C Run the adjoint (inverse) model MOODY. 
C The primary input is the velocity errors uerr_best, verr_best.
C The output is DomegaDrho, the gradient of the cost function. 
C
C History: Written Oct/Nov 1998 by CGH and DRL 
C          added DomegaDrho initialization 16 Jan 1999.
C          set DomegaDrho(i,k=1) = 0  16 Jan 1999. 
C          Added Wrho_t  20 Aug 1999 CGH 
C
      subroutine moody_run(casco_iter)
C
      include 'CASCO.DIM'
      include 'mesh.h'
      include 'saco.h'
      include 'startstop.h'
      include 'statarrays.h'
      include 'constants.h'
      include 'weights.h'
      include 'observes.h' 
      include 'errors.h'
      include 'covariance.h' 
      include 'gradcost.h'
      include 'sacotime.h'
      include 'deltabcs.h' 
      include 'wrho.h'

      integer casco_iter

C local variables  
      integer i,j,k
      integer jay
      real    tw
      integer obs_ptr(2) 
      integer eobs_ptr(2) 
      real    seckd
      integer kd
      real    r(NRDIM)

       write(*,*)'MOODY_RUN: casco_iter',casco_iter 

C PREPARATION FOR TIME STEPPING SECTION
C    what needs to be done ? 
        do i = 1, nn
	   zetaold(i) = 0.0
	   zetamid(i) = 0.0
	   zetanew(i) = 0.0
	   umid(i) = 0. 
	   vmid(i) = 0. 
	enddo
	do j = 1, nnv
	   do i = 1, nn
	      uzmid(i,j) = 0.0
	      vzmid(i,j) = 0.0 
	      uznew(i,j) = 0.0 
	      vznew(i,j) = 0.0
	   enddo
	enddo

CGH 16 Jan 1999
C Set DomegaDrho = dbc_best*Wrho as the starting point for 
C the next iteration.  Wrho contains the weight information 
C related the boundary control parameters w0 and w1. 
C Moved here from adjust_soln.  
C 
C This is overkill since  DomegaDrho = 0 in initial_conditions. 
C But I do it here anyway. 
C
      if(casco_iter.eq.1) then 
 	 do j = 1, ntbc
 	   do i = 1, nr
 	      DomegaDrho(i,j) = 0.
  	   enddo
 	 enddo
       elseif(casco_iter.gt.1) then 
          call xWRHO(DomegaDrho,dbc_best,Wrho,Wrho_t,IQrho,JQrho
     $               ,nr,ntbc)
       else
          write(*,*)'MOODY_RUN: casco_iter < 1 ', casco_iter
          write(*,*)'Illegal state.  STOP NOW' 
          stop
       endif
C
      if(casco_iter.le.2) then 
         write(22,*)'casco_iter ', casco_iter
         write(22,*)'BASISOBS '
         do i = nobs,1,-1
 	    write(22,1001)i,(basisobs(j,i),j=1,6)
         enddo
1001  format(1x,i4,1x,f8.1,1x,3e12.4,1x,f8.1,1x,e12.4)
         write(22,*)
         write(22,*)'VWEIGHT'
         do i = nobs,1,-1
 	    do j = 1, nnv
 	       write(22,*)i,j,vweight(i,j)
 	    enddo
         enddo
         write(22,*)
         write(22,*)'W_delta'
         do i = 1, nobs
            write(22,*)i,W_delta(i)
         enddo
      endif 

C Preparation for time stepping 
      write(*,*)'nsteps ', nsteps
      kd = kdf
      seckd = secf
      obs_ptr(1) = nobs+1
      obs_ptr(2) = nobs+1 
      eobs_ptr(1) = n_eobs+1
      eobs_ptr(2) = n_eobs+1 
      k = 0 
      call outputm2(k,kd,seckd,delt,nn,nnv,ne,x,y,in,hdown,
     &     zetamid,umid,vmid,zmid,uzmid,vzmid,zetaold,domegadxi,
     $     domegadrho,lrvs,nr,nsteps,ntbc,casco_iter)

C TIME STEPPING LOOP

      do k = nsteps, 1,-1 
         seckd = seckd - delt
         call up_date2(kd,seckd)

         call verticalm4(nn,nnv,ne,in,delt,sv,cxmat,cymat,
     &        dxmat,dymat,iq,jq,
     &        d1,d2,d3,dz,zetamid,uzmid,vzmid,uznew,vznew,
     &        basisobs,vweight,uerr_best,verr_best,W_delta,k,obs_ptr) 

         call elevationm4(nn,nnv,in,nhbw,theta,
     &        zetaold,zetamid,uzmid,vzmid,uznew,vznew,
     &        sh,iq,jq,G0xmat,G0ymat,Amat,Bmat,dz,
     &        zetanew,ebasisobs,eerr_best,W_eps,k,eobs_ptr)

	 Call EMATpreMULT(EMAT,zetanew,R,LRVS,NN,NHBW,NR,NBWE)
	 jay = pkj(k)
	 tw  = kjweight(k)
	 do i = 1, nr
	    DomegaDxi(i,k) = -r(i)
	    DomegaDrho(i,jay) = DomegaDrho(i,jay) -r(i)*tw
	    DomegaDrho(i,jay-1) = DomegaDrho(i,jay-1) -r(i)*(1.0-tw)
	 enddo



C*** update 
C *OLD<=*MID:  only need to do adjoint elevations
       do i = 1, nn
	  zetaold(i) = zetamid(i)
       enddo
C *MID<=*NEW: 
       do i = 1, nn
	  zetamid(i) = zetanew(i)
	  zetanew(i) = 0.0
       enddo
       do j = 1,nnv
	  do i = 1,nn
	     uzmid(i,j) = uznew(i,j)
	     vzmid(i,j) = vznew(i,j)
	     uznew(i,j) = 0.0
	     vznew(i,j) = 0.0
	  enddo
       enddo
C Result: Arrays at the end of the current time step
C                                      are now stored as *MID
C
C-----------------------------------------------------------------------
C compute vertically averaged velocities at current time (*MID)
C
          call vertavg(uzmid,umid,zmid,nn,nnv,nndim)
          call vertavg(vzmid,vmid,zmid,nn,nnv,nndim)

      call outputm2(k,kd,seckd,delt,nn,nnv,ne,x,y,in,hdown,
     &     zetamid,umid,vmid,zmid,uzmid,vzmid,zetaold,domegadxi,
     $     domegadrho,lrvs,nr,nsteps,ntbc,casco_iter)
      
       enddo
C end of time stepping loop
C
CGH set DomegaDrho(i,k=1) = 0 
C   This ensures that the initial elevations are zero and reduces
C   startup shocks.  The spin-up time is BCS_AVG_TSCALE.  
       k = 1 
       do i = 1, nr
          DomegaDrho(i,k) = 0.0
       enddo

C
C  do output stuff 
      k = -1 
      call outputm2(k,kd,seckd,delt,nn,nnv,ne,x,y,in,hdown,
     $     zetamid,umid,vmid,zmid,uzmid,vzmid,zetaold,domegadxi,
     $     domegadrho,lrvs,nr,nsteps,ntbc,casco_iter)

      write(*,*)'MOODY_RUN. done ' 

      RETURN
      END

