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
  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,fulllin,debugco,aknorm,usewdpij
  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
  real(r_kind) pob,time_co
  real(r_quad) val1,valx,val1dp
  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_kind),allocatable,dimension(:)   :: vali,validp
  real(r_kind),allocatable,dimension(:)   :: val_ret
  real(r_kind),allocatable,dimension(:)   :: conew,dfga,alt,cogesnew,resnew
  real(r_kind),allocatable,dimension(:,:,:) :: rcodpinc
  type(colvk_ob_type), pointer :: colvkptr
  integer :: entry=0
  save entry
  entry=entry+1

!  If no co observations return
  if(.not. associated(colvkhead))return
! Retrieve pointers
! Simply return if any pointer not found
  ier=0 
  nllpoint=0
  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
  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
!
! 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)
  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)
!      write(500+mype,*)mype,'j1',j1,j2,j3,j4
!      call flush(500+mype)
!     if(j1<1.or.j2<1.or.j3<1.or.j4<1)then
!       call stop2(432)
!     endif


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


        allocate(vali(colvkptr%nlco))
        allocate(validp(colvkptr%nlco))
        allocate(coak(colvkptr%nlco))
        allocate(conew(colvkptr%nlco))
        allocate(cogesnew(colvkptr%nlco))
        allocate(resnew(colvkptr%nlco))
        allocate(dfga(colvkptr%nlco))
        allocate(alt(colvkptr%nlco))
        allocate(val_ret(colvkptr%nlco))
        if(usewdpij)then
          do k=1,colvkptr%nlco   ! loop over MOPITT ave. ker. contribution levels 
!           if(mype.eq.0)then
!             write(980,*)k,'nsub',colvkptr%nsub(k)
!           endif
            val1dp=0.0
            do kk=1,colvkptr%nsub(k)
!             if(mype.eq.0)then
!               write(980,*)'kk',kk,' ksub ',colvkptr%ksub(k,kk),colvkptr%wdpij(k,kk,1)
!             endif
              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
           validp(k)=val1dp
          end do
          do k=1,colvkptr%nlco
            vali(k)=validp(k)
          enddo
        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
        if(fulllin)then
!         keep fully linear no log or exp
          do k=1,colvkptr%nlco
!            conew(k)=colvkptr%coges(k)+vali(k)
!            if(conew(k)<0.0)then
!               write(400+mype,*)entry,nllpoint,'conewmod neg',k,conew(k),'ges',colvkptr%coges(k),'vali',vali(k)
!            endif
            dfga(k)=vali(k)/colvkptr%comodges(k)
 
            if(debugco)then
              if(dfga(k).ne.0)then
                if(abs(dfga(k))>.6)then
                  write(500+mype,*)entry,nllpoint,' big fraction ',k,dfga(k),'vali',vali(k),' ges ',colvkptr%comodges(k)
                endif
                write(380+mype,*)entry,'inc scaled ',k,dfga(k),' vali ',vali(k),' modges ',colvkptr%comodges(k)
              endif
            endif
          end do
          alt=matmul(colvkptr%ak,dfga)
!         need to put back to non fractional units for change of res
          do k=1,colvkptr%nlco
             alt(k)=alt(k)*colvkptr%coges(k)
             if(debugco)then
               if(alt(k).ne.0)then
                write(380+mype,*)entry,'alt descaled ',k,alt(k)
               endif
             endif
             conew(k)=colvkptr%coges(k)+alt(k)
!             if(conew(k)<0.0)then
!               write(400+mype,*)entry,nllpoint,'conew neg',k,conew(k),'ges',colvkptr%coges(k),'val1',alt(k)
!               conew(k)=1.e-20
!             endif
             val1=alt(k)
             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
!   ajl should have consistent units of val1 here to be same as res
!             if(abs(val1)>.01*abs(colvkptr%res(k)))then
!             write(800+mype,*)'val1',val1,'res',colvkptr%res(k)
!             call flush(800+mype)
!             endif
                    val1=alt(k)-colvkptr%res(k)
                    if(debugco)write(380+mype,*)entry,' new -res ',val1,' alt ',alt(k),' old res ',colvkptr%res(k)
                    val1=val1/colvkptr%coges(k) ! make fractional because of log
                    if(debugco)write(420+mype,*)entry,nllpoint,' val 1 scaled ',k,val1
!                    val1=-resnew(k)/colvkptr%coges(k) ! make fractional because of log

                    valx=val1*colvkptr%err2(k) 
                    valx=valx*colvkptr%raterr2(k)
!                    if(nllpoint<4)then
!             write(400+mype,*)'k',k,'res',colvkptr%res(k) ,'err2',colvkptr%err2(k),'raterr2',colvkptr%raterr2(k)
!                   endif
                 end if
              endif
              val_ret(k)=valx  
              if(debugco)write(420+mype,*)entry,nllpoint,' valret',k,val_ret(k)
            endif 
          end do
        else
          do k=1,colvkptr%nlco
            conew(k)=colvkptr%comodges(k)+vali(k)
            if(conew(k)<0.0)then
!               write(400+mype,*)entry,nllpoint,' conew neg',k,conew(k),'ges',colvkptr%comodges(k),'vali',vali(k)
               conew(k)=1.e-16
            endif
!           write(400+mype,*)'conew',k,conew(k),'ap',colvkptr%ap(k),'vali',vali(k),'ges',colvkptr%comodges(k)
!           call flush(400+mype)
            dfga(k)=log(conew(k))-log(colvkptr%ap(k))
!           write(400+mype,*)'dfga',k,dfga(k),' new ',conew(k),' ap ',colvkptr%ap(k)
!           call flush(400+mype)
          end do
          alt=matmul(colvkptr%ak,dfga)
          do k=1,colvkptr%nlco
            cogesnew(k)=colvkptr%ap(k)*exp(alt(k))
            resnew(k)=colvkptr%coobs(k)-cogesnew(k)
!            if(nllpoint.eq.19)then
!              write(500+mype,*)entry,'find resnew',k,resnew(k),' cogesnew ',cogesnew(k),' alt ',alt(k), &
!               ' obs ',colvkptr%coobs(k)
!            endif
!           write(400+mype,*)'alt',alt(k),' coges new ',cogesnew(k),' old ',colvkptr%coges(k),' diff ', &
!           cogesnew(k)-colvkptr%coges(k)
!           write(400+mype,*)'resnew ',k,resnew(k),' old ',colvkptr%res(k),' diff ', &
!           resnew(k)-colvkptr%res(k)
          enddo
!         val1 is convolved increment on airs model layers

!         Averaging kernel  

          do k=1,colvkptr%nlco   ! loop over MOPITT retrieval levels
            val1=cogesnew(k)-colvkptr%coges(k)

            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
!   ajl should have consistent units of val1 here to be same as res
!             if(abs(val1)>.01*abs(colvkptr%res(k)))then
!             write(800+mype,*)'val1',val1,'res',colvkptr%res(k)
!             call flush(800+mype)
!             endif
!                    val1=val1-colvkptr%res(k)
                    !val1=resnew(k)/colvkptr%coges(k) ! make fractional because of log
                    val1=-resnew(k)/colvkptr%coges(k) ! make fractional because of log

                    valx=val1*colvkptr%err2(k) 
                    valx=valx*colvkptr%raterr2(k)
!                    if(nllpoint<4)then
!              write(400+mype,*)'k',k,'res',colvkptr%res(k) ,'err2',colvkptr%err2(k),'raterr2',colvkptr%raterr2(k), &
!                'ges',colvkptr%coges(k) 
!                   endif
                 end if
              endif
              val_ret(k)=valx  
!              if(nllpoint.eq.19)then
!                write(500+mype,*)entry,'val_ret',k,val_ret(k),' resnew ',resnew(k),' coges ',colvkptr%coges(k)
!                write(500+mype,*)' err2 ',colvkptr%err2(k),' raterr2 ',colvkptr%raterr2(k)
!              endif
            endif 
          enddo ! k
        endif

!  Averaging kernel First - spread values to ak contribution levels 

        if(l_do_adjoint)then 
!          if(colvkptr%ak(1,1)<-100.)then
!            coak=val_ret
!            if(logicalairsco)then
              !do k=1,colvkptr%nlco
!                coak(k)=coak(k)/dpairs(k)
!               take mass back out
!              end do
!            endif
!          else
              do k=1,colvkptr%nlco  !loop over ak levels 
                coak(k)=zero_quad
!                if(nllpoint<20)then
!                write(800+mype,*)'ak',k,maxval(colvkptr%ak(:,k)),minval(colvkptr%ak(:,k))
!                endif
                do j=1,colvkptr%nlco  !loop over profile levels  
                  coak(k)=coak(k)+colvkptr%ak(j,k)*val_ret(j) ! Contribution to kth ak level from jth retrieval level
               enddo
               if(debugco)write(420+mype,*)entry,nllpoint,'coak',k,coak(k)
!               if(abs(coak(k))>.6)then
!                 write(500+mype,*)entry,nllpoint,' big coak ',k,coak(k),' ges ',colvkptr%comodges(k),' res ',colvkptr%res(k)
!                 write(500+mype,*)'sum ak ',sum(colvkptr%ak(:,k))
!               endif
               coak(k)=coak(k)*colvkptr%comodges(k) ! defractional back to mass of co
               if(debugco)write(420+mype,*)entry,nllpoint,' descale ',coak(k)
               coak(k)=coak(k)/dpairs(k) ! put back to vmr
               if(debugco)write(420+mype,*)entry,nllpoint,' vmr ',coak(k),' modges/dp ',colvkptr%comodges(k)/dpairs(k)
               if(aknorm>0.)then
                 coak(k)=coak(k)*aknorm
               endif
!                if(nllpoint<20)then
!               write(800+mype,*)'coak',k,coak(k),'ges',colvkptr%comodges(k),' aknorm ',colvkptr%aknorm
!               endif
             enddo 
!          endif
!
! Adjoint of interpolation - spreads each ave. kernel level to interpolant gridpoints  
            if(.not.usewdpij)then

              do kk=colvkptr%nlco,1,-1    !loop over averaging kernel levels 
                 pob = colvkptr%prs(kk)
                 k1=int(pob)
                 k2=min(k1+1,nsig)
                 w1=colvkptr%wij(1,kk)
                 w2=colvkptr%wij(2,kk)
                 w3=colvkptr%wij(3,kk)
                 w4=colvkptr%wij(4,kk)
                 w5=colvkptr%wij(5,kk)
                 w6=colvkptr%wij(6,kk)
                 w7=colvkptr%wij(7,kk) 
                 w8=colvkptr%wij(8,kk) 
                 rco(j1,k1)  =  rco(j1,k1) + coak(kk)*w1
                 rco(j2,k1)  =  rco(j2,k1) + coak(kk)*w2
                 rco(j3,k1)  =  rco(j3,k1) + coak(kk)*w3
                 rco(j4,k1)  =  rco(j4,k1) + coak(kk)*w4
                 rco(j1,k2)  =  rco(j1,k2) + coak(kk)*w5
                 rco(j2,k2)  =  rco(j2,k2) + coak(kk)*w6
                 rco(j3,k2)  =  rco(j3,k2) + coak(kk)*w7
                 rco(j4,k2)  =  rco(j4,k2) + coak(kk)*w8
                 if(nllpoint<10)then
!                   write(900+mype,*)entry,'k',kk,'tco inc1 ',j1,k1,coak(kk)*w1,' rco ',rco(j1,k1),' sco ',sco(j1,k1)
!                   write(900+mype,*)kk,'tco inc2 ',j1,k1,coak(kk)*w2
                   !write(900+mype,*)kk,'tco inc3 ',j1,k1,coak(kk)*w3
                   !write(900+mype,*)kk,'tco inc4 ',j1,k1,coak(kk)*w4
                   !write(900+mype,*)kk,'tco inc5 ',j1,k2,coak(kk)*w5
                   !write(900+mype,*)kk,'tco inc6 ',j1,k2,coak(kk)*w6
                   !write(900+mype,*)kk,'tco inc7 ',j1,k2,coak(kk)*w7
                   !write(900+mype,*)kk,'tco inc8 ',j1,k1,coak(kk)*w8
                   !call flush(900+mype)
                 endif
              enddo  ! k2
           else
              do k=1,colvkptr%nlco   ! loop over MOPITT ave. ker. contribution levels 
!               if(mype.eq.0)then
!                 write(980,*)k,'nsub',colvkptr%nsub(k)
!               endif
                do kk=1,colvkptr%nsub(k)
!                 if(mype.eq.0)then
!                   write(980,*)'kk',kk,' ksub ',colvkptr%ksub(k,kk),colvkptr%wdpij(k,kk,1)
!                 endif
                  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)+coak(kk)*w1
                  rcodp(j2,k1)=rcodp(j2,k1)+coak(kk)*w2
                  rcodp(j3,k1)=rcodp(j3,k1)+coak(kk)*w3
                  rcodp(j4,k1)=rcodp(j4,k1)+coak(kk)*w4
                end do
              end do
            endif


        deallocate(coak,vali,val_ret,validp)
        deallocate(alt,dfga,conew,cogesnew,resnew)

        endif ! l_do_adjoint
        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(.not.rcop(i,j,k)<=0.0.and..not.rcop(i,j,k)>=0.0)then
            write(6,*)'NaN intco ',i,j,k,'ij',ij
            call flush(6)
            call stop2(444)
          endif
        enddo
      enddo
    enddo
!    write(900+mype,*)'bot incodplev rcoinc',maxval(abs(rcodpinc)),' max ',maxval(rco),minval(rco)
  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)
  return
end subroutine intcolev_


end module intcomod
