module intcomod

!$$$ module documentation block
!           .      .    .                                       .
! module:   intcomod    module for intco and its tangent linear intco_tl
!   prgmmr:
!
! abstract: module for intco and its tangent linear intco_tl
!
! program history log:
!   2005-05-13  Yanqiu zhu - wrap intoz and its tangent linear intoz_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-11-26  Todling - remove intoz_tl; add interface back
!   2009-08-13  lueken - update documentation
!   2010-06-02  tangborn - converted intoz into intco 
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - implemented obs adjoint test  
!
! subroutines included:
!   sub intco_
!   sub intcolev_
!
! variable definitions:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC intco

interface intco; module procedure &
          intco_
end interface

contains

subroutine intco_(colvkhead,rval,sval)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intco       call individual carbon monoxide obs operators
!   prgmmr: todling       org: np23                date: 2008-11-28
!
! abstract:  This routine calls the individual components of the 
!            carbon monoxide observation operator.
!
! program history log:
!   2008-11-28  todling
!   2009-01-08  todling - remove reference to ozohead
!   2010-05-13  todling - update to use gsi_bundle
!   2010-06-02  tangborn - made version for carbon monoxide
!
!   input argument list:
!     colvkhead  - level carbon monoxide obs type pointer to obs structure for MOPITT
!     sco     - carbon monoxide increment in grid space
!
!   output argument list:
!     rco    - carbon monoxide results from observation operator 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use obsmod, only: colvk_ob_type
  use gsi_bundlemod, only: gsi_bundle
  implicit none

! Declare passed variables
  type(colvk_ob_type),pointer,intent(in   ) :: colvkhead
  type(gsi_bundle),intent(in   ) :: sval
  type(gsi_bundle),intent(inout) :: rval

  call intcolev_(colvkhead,rval,sval)

end subroutine intco_

subroutine intcolev_(colvkhead,rval,sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intco       apply nonlin qc obs operator for carbon monoxide
!   prgmmr: derber           org: np23                date: 1995-07-11
!
! abstract:  This routine applies the observation operator (forward
!            model) and adjoint of this operator for ozone observations
!            with the addition of nonlinear qc.
!
! program history log:
!   1995-07-11  derber
!   2010-06-07  tangborn - carbon monoxide based on ozone code
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - introduced ladtest_obs         
!
!   input argument list:
!     colvkhead  - level carbon monoxide obs type pointer to obs structure
!     sco     - carbon monoxide increment in grid space
!
!   output argument list:
!     rco     - carbon monoxide results from observation operator 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use kinds, only: r_kind,i_kind,r_quad,r_double
  use obsmod, only: colvk_ob_type,lsaveobsens,l_do_adjoint
  use gridmod, only: lat2,lon2,nsig
  use jfunc, only: jiter,xhat_dt,dhat_dt
  use constants, only: one,zero,r3600,zero_quad
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_4dvar, only: ladtest_obs
  use mpimod, only : mype
  use modairsco, only :dpairs,logicalairsco,debugco,usewdpij
  use modairsco, only : nllpointmin
  implicit none

! Declare passed variables
  type(colvk_ob_type),pointer,intent(in   ) :: colvkhead
  type(gsi_bundle)          ,intent(in   ) :: sval
  type(gsi_bundle)          ,intent(inout) :: rval

! Declare local variables
  integer(i_kind) i,j,ij,ier,istatus
  integer(i_kind) k,k1,k2,j1,j2,j3,j4,kk,iz1,iz2,j1x,j2x,j3x,j4x,nllpoint
  integer j1save,j2save,j3save,j4save
  real(r_kind) pob,time_co,maxper,per
!  !real(r_quad) val1,valx,val1dp
!  real(r_quad) val1q,valxq,val1dpq
  real(r_double) val1,valx,val1dp
  real(r_double) val1q,valxq,val1dpq
  real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8
  real(r_kind),pointer,dimension(:) :: xhat_dt_co
  real(r_kind),pointer,dimension(:) :: dhat_dt_co
  real(r_kind),pointer,dimension(:,:,:)  :: scop
  real(r_kind),pointer,dimension(:,:,:)  :: rcop
  real(r_kind),allocatable,dimension(:,:) :: sco
  real(r_kind),allocatable,dimension(:,:) :: rco
  real(r_kind),allocatable,dimension(:,:) :: rcodp
  real(r_kind),allocatable,dimension(:,:,:) :: rcoinc
!  real(r_kind),allocatable,dimension(:)   :: coak
!  real(r_quad)                             :: coak
  real(r_double)                             :: coak
  real(r_kind),allocatable,dimension(:)   :: vali,valxk
  real(r_kind),allocatable,dimension(:)   :: val_ret
  real(r_kind),allocatable,dimension(:,:,:) :: rcodpinc
! put in choice to do error before ak backwards like original intco
  type(colvk_ob_type), pointer :: colvkptr
  logical docohmatr
  character*10 ccohmatr
  save docohmatr
  integer :: entry=0
!  real*8 wall1,wall2,wallco
!  save entry,wallco
!  data wallco/0.0/


!  If no co observations return
  if(.not. associated(colvkhead))return
  entry=entry+1
!  call getwall(wall1)
  if(entry.eq.1)then
    ccohmatr=' '
    call getenv('COHMATR',ccohmatr)
    if(ccohmatr.eq.'YES')then
      docohmatr=.true.
    else
      docohmatr=.false.
    endif
    write(6,*)'docohmatr',docohmatr
  endif
    
  j1save=-9999
  j2save=-9999
  j3save=-9999
  j4save=-9999
! Retrieve pointers
! Simply return if any pointer not found
  ier=0 
  nllpoint=0
  if(nllpointmin>0)print *,mype,'top intco nllpointmin ',nllpointmin
  call gsi_bundlegetpointer(sval,'co',scop,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'co',rcop,istatus);ier=istatus+ier
  if(ier/=0)return
!   write(6,*)mype,'top of intco'
!   call flush(6)
!   write(980+mype,*)'top of intco'
!  call flush(980+mype)

! Can't do rank-2 pointer into rank-2, therefore, allocate work space
!    write(6,*)'usewdpij',usewdpij
!    call flush(6)
  allocate(sco(lat2*lon2,nsig),rco(lat2*lon2,nsig))
  allocate(rcoinc(lat2,lon2,nsig))
  allocate(rcodpinc(lat2,lon2,nsig)) 
  allocate(rcodp(lon2*lat2,nsig))
  rcoinc=0.0
  rcodpinc=0.0
  do k=1,nsig
     ij=0
     do j=1,lon2
        do i=1,lat2
           ij=ij+1
           sco(ij,k) = scop(i,j,k)
           rco(ij,k) = rcop(i,j,k)
           rcodp(ij,k) = rcop(i,j,k)
        enddo
     enddo
  enddo
  if(debugco)then
    write(800+mype,*)'entry',entry
    write(800+mype,*)'top of intco sval',maxval(sco),minval(sco),' rval ',maxval(rco),minval(rco)
    call flush(800+mype)
  endif
!  write(800+mype,*)entry,'top of intco ',maxval(sco),minval(sco),' rval ',maxval(rco),minval(rco)
!  call flush(800+mype)
!  write(6,*)mype,'top of intco sval ',maxval(sco),minval(sco),' rval ',maxval(rco),minval(rco)
!  call flush(6)
!
! MOPITT CARBON MONOXIDE: LAYER CO 
!
! Loop over carbon monoxide observations.
  colvkptr => colvkhead
!  write(800+mype,*)'associated colvkptr ',associated(colvkptr)
!  write(980+mype,*)'associated colvkptr ',associated(colvkptr)
  maxper=0.0
  do while (associated(colvkptr))
    nllpoint=nllpoint+1

!   Set location
    j1=colvkptr%ij(1)
    j2=colvkptr%ij(2)
    j3=colvkptr%ij(3)
    j4=colvkptr%ij(4)


!   Accumulate contribution from layer observations
!   This repeats the algorithm inside intrp3co
!   with several of the terms already calculated 


    allocate(vali(colvkptr%nlco))
    allocate(valxk(colvkptr%nlco))
    allocate(val_ret(colvkptr%nlco))
    if(usewdpij)then
      do k=1,colvkptr%nlco   ! loop over MOPITT ave. ker. contribution levels 
        val1dp=0.0
        do kk=1,colvkptr%nsub(k)
          k1=colvkptr%ksub(k,kk)
          w1=colvkptr%wdpij(k,kk,1)
          w2=colvkptr%wdpij(k,kk,2)
          w3=colvkptr%wdpij(k,kk,3)
          w4=colvkptr%wdpij(k,kk,4)
          val1dp=val1dp+ w1* sco(j1,k1)+ &
               w2* sco(j2,k1)+ &
               w3* sco(j3,k1)+ &
               w4* sco(j4,k1)
        end do
        vali(k)=val1dp
      end do
    else

      do k=1,colvkptr%nlco   ! loop over MOPITT ave. ker. contribution levels 
         
        pob = colvkptr%prs(k)
      
        k1=int(pob)
        k2=min(k1+1,nsig)
        w1=colvkptr%wij(1,k)
        w2=colvkptr%wij(2,k)
        w3=colvkptr%wij(3,k)
        w4=colvkptr%wij(4,k)
        w5=colvkptr%wij(5,k)
        w6=colvkptr%wij(6,k)
        w7=colvkptr%wij(7,k)
        w8=colvkptr%wij(8,k)
        val1=  w1* sco(j1,k1)+ &
               w2* sco(j2,k1)+ &
               w3* sco(j3,k1)+ &
               w4* sco(j4,k1)+ &
               w5* sco(j1,k2)+ &
               w6* sco(j2,k2)+ &
               w7* sco(j3,k2)+ & 
               w8* sco(j4,k2)
       vali(k)=val1
      enddo
      if(logicalairsco)then
        do k=colvkptr%nlco,2, -1
          vali(k)=.5*(vali(k)+vali(k-1))*dpairs(k)
        end do
        vali(1)=vali(2)
      endif
    endif
!   for airs need val1 bar *dp
!   new guess
!      alt=matmul(colvkptr%cohmat,vali)
    do k=1,colvkptr%nlco
!      val1=zero_quad
      val1=zero
      do j=1,colvkptr%nlco
        val1=val1+colvkptr%cohmat(k,j)*vali(j)
!        if(.not.val1<=0.0.and..not.val1>=0.0)then
!          write(6,*)mype,'val1',k,j,val1,'cohmat',colvkptr%cohmat(k,j),'vali',vali(j)
!          call flush(6)
!          call stop2(888)
!        endif
      end do
      if (lsaveobsens) then
        colvkptr%diags(k)%ptr%obssen(jiter)=val1*colvkptr%err2(k)*colvkptr%raterr2(k)
      else
        if (colvkptr%luse) colvkptr%diags(k)%ptr%tldepart(jiter)=val1
      endif
      if (l_do_adjoint) then
        if (lsaveobsens) then
          valx = colvkptr%diags(k)%ptr%obssen(jiter)
        else
          if( ladtest_obs ) then
            valx = val1
          else
!            per=abs(val1)/max(1.e-30,abs(colvkptr%res(k)))
!            maxper=max(per,maxper)
!            if(nllpoint.eq.1)then
!              write(6,*)mype,'val1',k,val1,'res',colvkptr%res(k)
!              call flush(6)
!            endif
            val1=val1-colvkptr%res(k)
            valx=val1*colvkptr%err2(k) 
            valx=valx*colvkptr%raterr2(k)
            valxk(k)=valx
!            if(nllpoint.eq.1)then
!              write(6,*)mype,'valx',k,valx
!            endif
          endif
        endif
      endif
    end do
    if(l_do_adjoint)then
      do k=1,colvkptr%nlco
!        coak=zero_quad
        coak=zero
        if(docohmatr)then
          do j=1,colvkptr%nlco  !loop over profile levels  
            coak=coak+colvkptr%cohmatr(j,k)*valxk(j)
          enddo
        else
          do j=1,colvkptr%nlco  !loop over profile levels  
            coak=coak+colvkptr%cohmat(j,k)*valxk(j)
          enddo
        endif
        val_ret(k)=coak
!        if(abs(coak)>1.e5)then
!          write(6,*)mype,'coak',k,coak
!          do j=1,colvkptr%nlco
!            write(6,*)mype,k,'valxk',j,valxk(j),'cohmat',colvkptr%cohmat(j,k),'prod',colvkptr%cohmat(j,k)*valxk(j)
            !call flush(6)
!          end do
!        endif
!            if(nllpoint.eq.1)then
!               write(6,*)mype,'val_ret',k,val_ret(k)
!               call flush(6)
!            endif
      end do
    endif
    if(l_do_adjoint)then 
        if(nllpoint.eq.nllpointmin.and.mype.eq.16)then
!          write(6,*)mype,'val_ret',k,val_ret(k)
          j1save=j1
          j2save=j2
          j3save=j3
          j4save=j4
        endif
      do k=1,colvkptr%nlco   ! loop over MOPITT ave. ker. contribution levels 
        do kk=1,colvkptr%nsub(k)
          k1=colvkptr%ksub(k,kk)
          w1=colvkptr%wdpij(k,kk,1)
          w2=colvkptr%wdpij(k,kk,2)
          w3=colvkptr%wdpij(k,kk,3)
          w4=colvkptr%wdpij(k,kk,4)
          rcodp(j1,k1)=rcodp(j1,k1)+val_ret(k)*w1
          rcodp(j2,k1)=rcodp(j2,k1)+val_ret(k)*w2
          rcodp(j3,k1)=rcodp(j3,k1)+val_ret(k)*w3
          rcodp(j4,k1)=rcodp(j4,k1)+val_ret(k)*w4
          go to 801
          if(abs(rcodp(j1,k1))>1.e12)then
            write(6,*)mype,'big rcodp j1',j1,nllpoint,'k',k,kk,'k1',k1,'val_ret',val_ret(k),'w1',w1
            call flush(6)
          endif
          if(abs(rcodp(j2,k1))>1.e12)then
            write(6,*)mype,'big rcodp j2',j2,nllpoint,'k',k,kk,'k1',k1,'val_ret',val_ret(k),'w2',w2
            call flush(6)
          endif
          if(abs(rcodp(j3,k1))>1.e12)then
            write(6,*)mype,'big rcodp j3',j3,nllpoint,'k',k,kk,'k1',k1,'val_ret',val_ret(k),'w3',w3
            call flush(6)
          endif
          if(abs(rcodp(j4,k1))>1.e12)then
            write(6,*)mype,'big rcodp j4',j4,nllpoint,'k',k,kk,'k1',k1,'val_ret',val_ret(k),'w4',w4
            call flush(6)
          endif
801       continue
        end do
      end do
    endif
    deallocate(vali,val_ret,valxk)
    colvkptr => colvkptr%llpoint

! End loop over observations
  enddo

! Copy output and clean up 
!  write(900+mype,*)'bot inco before  max ',maxval(rco),minval(rco)
  if(usewdpij)then
    do k=1,nsig
      ij=0
      do j=1,lon2
        do i=1,lat2
          ij=ij+1
          rcodpinc(i,j,k)=rcodp(ij,k)-rcop(i,j,k)
          rcop(i,j,k) = rcodp(ij,k)
!          if(nllpoint.eq.1)then
!          endif
!          if(abs(rcodpinc(i,j,k))>1.e11)then
!           write(6,*)mype,'big rcodpinc ',i,j,k,'ij',ij,'rcodp',rcodp(ij,k),rcop(i,j,k)
!           call flush(6)
!          endif
!          if(ij.eq.j1save.or.ij.eq.j2save.or.ij.eq.j3save.or.ij.eq.j4save)then
!            if(rcodpinc(i,j,k).ne.0.0)then
!            !write(6,*)k,'i,j',i,j,'rcodpinc ',rcodpinc(i,j,k)
!            endif
!          endif
!          if(.not.rcop(i,j,k)<=0.0.and..not.rcop(i,j,k)>=0.0)then
!            write(6,*)mype,'NaN intco ',i,j,k,'ij',ij
!            call flush(6)
!            call stop2(444)
!          endif
        enddo
      enddo
!      write(900+mype,*)'rocdpinc ',k,maxval(rcodpinc(:,:,k)),minval(rcodpinc(:,:,k))
    enddo
!    write(900+mype,*)'bot incodplev rcoinc',maxval(abs(rcodpinc)),' max ',maxval(rco),minval(rco)
!    write(6,*)mype,'bot incodplev rcoinc',maxval(abs(rcodpinc)),' max ',maxval(rco),minval(rco)
!    call flush(6)
  else
    do k=1,nsig
      ij=0
      do j=1,lon2
        do i=1,lat2
          ij=ij+1
          rcoinc(i,j,k)=rco(ij,k)-rcop(i,j,k)
          rcop(i,j,k) = rco(ij,k)
        enddo
      enddo
    enddo
!    write(900+mype,*)'bot incolev rcoinc',maxval(abs(rcoinc)),' max ',maxval(rco),minval(rco)
  endif
!  write(900+mype,*)'bottom intco ',nllpoint
!  call flush(900+mype)
  deallocate(sco,rco)
  deallocate(rcoinc,rcodpinc,rcodp)

! End of routine
!  write(6,*)mype,'bottom intco',nllpoint
  !call flush(6)
!  write(6,*)mype,'maxper',maxper*100.
!  call flush(6)
!  call getwall(wall2)
!  wallco=wallco+wall2-wall1
!  write(500+mype,*)'wallco',entry,wallco/60.
!  call flush(500+mype)
  return
end subroutine intcolev_


end module intcomod
