!-------------------------------------------------------------------------
!     NASA/GSFC, Global Modeling and Assimilation Office, Code 601.1     !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE: model_ad:  
!
! !INTERFACE:

subroutine model_ad(xini,xobs,ldprt)

! !USES:

use kinds, only: r_kind,i_kind
use gsi_4dvar, only: nsubwin,nobs_bins,winlen,winsub,hr_obsbin
use gsi_4dvar, only: iadateend,idmodel
use constants, only: zero,r3600
use state_vectors
use geos_pertmod, only: ndtpert
use m_tick, only: tick
use timermod, only: timer_ini,timer_fnl

#ifdef GEOS_PERT
use prognostics, only: dyn_prog
use prognostics, only: prognostics_initial
use prognostics, only: prognostics_final
use geos_pertmod, only: gsi2pgcm
use geos_pertmod, only: pgcm2gsi
use m_model_ad, only: initial_ad
use m_model_ad, only: amodel_ad
use m_model_ad, only: final_ad
#endif /* GEOS_PERT */

use lag_fields, only: nlocal_orig_lag, ntotal_orig_lag
use lag_fields, only: lag_ad_vec,lag_tl_spec_i,lag_tl_spec_r
use lag_fields, only: lag_u_full,lag_v_full
use lag_fields, only: lag_ADscatter_stateuv
use lag_traj, only: lag_rk2iter_ad
! use lag_traj, only: lag_rk4iter_ad

implicit none

! !INPUT PARAMETERS:

type(state_vector), intent(in   ) :: xobs(nobs_bins) ! Adjoint state variable at observations times
logical           , intent(in   ) :: ldprt           ! Print-out flag

! !INPUT/OUTPUT PARAMETERS:

type(state_vector), intent(inout) :: xini(nsubwin)   ! Adjoint state variable at control times

! !DESCRIPTION: Run AGCM adjoint model.
!
! !REVISION HISTORY:
!
!  19Apr2007  tremolet - initial code
!  29May2007  todling  - add actual calls to interface and AGCM AD model
!  29Jun2007  todling  - adm verified against tlm
!  30Sep2007  todling  - add timer
!  30Apr2009  meunier  - add trajectory model for lagrangian data
!
!EOP
!-----------------------------------------------------------------------

! Declare local variables
character(len=*), parameter :: myname = 'model_ad'

type(state_vector) :: xx
integer(i_kind)    :: nstep,istep,nfrctl,nfrobs,ii,jj,ierr
integer(i_kind)    :: nymdi,nhmsi,ndt,dt
real(r_kind)       :: d0,tstep

! Temporary vector for lagrangian backward integration
real(r_kind),dimension(3):: ad_tmp_locvect

#ifdef GEOS_PERT
type(dyn_prog) :: xpert
#endif /* GEOS_PERT */

!******************************************************************************

! Initialize timer
call timer_ini('model_ad')

! Initialize variables
if (idmodel) then
   tstep  = r3600
   dt     = tstep
else
   ndt    = NINT(hr_obsbin*r3600/ndtpert)
   dt     = ndt*ndtpert
   tstep  = dt
endif
nstep  = NINT(winlen*r3600/tstep)
nfrctl = NINT(winsub*r3600/tstep)
nfrobs = NINT(hr_obsbin*r3600/tstep)
nymdi  =  iadateend/100
nhmsi  = (iadateend-100*nymdi)*10000

if (ldprt) write(6,'(a,3(1x,i4))')'model_ad: nstep,nfrctl,nfrobs=',nstep,nfrctl,nfrobs

call allocate_state(xx)
xx = zero
d0 = zero

! Initialize GCM ADM and its perturbation vector
#ifdef GEOS_PERT
if (.not.idmodel ) then
   call initial_ad
   call prognostics_initial(xpert)
endif
#endif /* GEOS_PERT */

ii=nobs_bins

! Post-process final state
if (nobs_bins>1) then
   if (ldprt) write(6,'(a,i8.8,1x,i6.6,2(1x,i4))')'model_ad: retrieving state nymdi,nhmsi,nobs_bins=',nymdi,nhmsi,nobs_bins
   call self_add(xx,xobs(nobs_bins))
   if (.not.idmodel) then
      d0=dot_product(xobs(1),xobs(1))
   endif
endif
! Run AD model
do istep=nstep-1,0,-1
   call tick(nymdi,nhmsi,-dt)

!  Apply AD model
#ifdef GEOS_PERT
   if (.not.idmodel) then
      call gsi2pgcm  (xx,xpert,'adm',ierr)                          ! T'(-1)
      call amodel_ad (xpert,nymdi,nhmsi,ntsteps=ndt,g5pert=.true.)  ! M'
      call pgcm2gsi  (xpert,xx,'adm',ierr)                          ! T'
   endif
#endif /* GEOS_PERT */

! Apply AD trajectory model (same time steps as obsbin)
   if (MOD(istep,nfrobs)==0 .and. ntotal_orig_lag>0) then
      ii=istep/nfrobs+1
      if (ldprt) write(6,'(a,i8.8,1x,i6.6,2(1x,i4))')'model_ad: trajectory model nymd,nhms,istep,ii=',nymdi,nhmsi,istep,ii
      if (ii<1.or.ii>nobs_bins) call abor1('model_ad: error xobs')
      ! Execute AD model for each balloon (loop step insensitive)
      do jj=1,nlocal_orig_lag
         ad_tmp_locvect = lag_ad_vec(jj,ii+1,:)
         ! if (.not.idmodel) then
         call lag_rk2iter_ad(lag_tl_spec_i(jj,ii,:),lag_tl_spec_r(jj,ii,:),&
           &ad_tmp_locvect(1),ad_tmp_locvect(2),ad_tmp_locvect(3),&
           &lag_u_full(:,:,ii),lag_v_full(:,:,ii))
         print '(A,I3,A,F16.6,F16.6)',"ADiter: ",ii," location",lag_ad_vec(jj,ii,1),lag_ad_vec(jj,ii,2)
         ! end if
         lag_ad_vec(jj,ii,:)=lag_ad_vec(jj,ii,:)+ad_tmp_locvect
      end do
      ! Give the sensitivity back to the GCM
      call lag_ADscatter_stateuv(xx%u,xx%v,ii)
      ! To not add the contribution 2 times
      lag_u_full(:,:,ii)=zero; lag_v_full(:,:,ii)=zero;
   endif

!  Post-process x_{istep}
   if (MOD(istep,nfrobs)==0) then
      ii=istep/nfrobs+1
      if (ldprt) write(6,'(2a,i8.8,1x,i6.6,2(1x,i4))')myname,': retrieving state nymd,nhms,istep,ii=',nymdi,nhmsi,istep,ii
      call self_add(xx,xobs(ii))
   endif

!  Apply control vector to x_{istep}
   if (MOD(istep,nfrctl)==0) then
      ii=istep/nfrctl+1
      if (ldprt) write(6,'(2a,i8.8,1x,i6.6,2(1x,i4))')myname,': adj adding WC nymd,nhms,istep,ii=',nymdi,nhmsi,istep,ii
      xini(ii)=xx
      d0=dot_product(xobs(ii),xx)
   endif

enddo

if(ldprt) print *, myname, ': total (gsi) dot product ', d0
#ifdef GEOS_PERT
if (.not.idmodel ) then
   call prognostics_final(xpert)
   call final_ad
endif
#endif /* GEOS_PERT */

call deallocate_state(xx)

! Finalize timer
call timer_fnl('model_ad')

return
end subroutine model_ad