subroutine intrp3trace1(f,g,dx,dy,dz,obstime,nlevs,thpresatobs,conversion,wdpij,maxsub,nsub,ksub,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intrp3trace    space-time linear interpolation for
! trace
!   prgmmr: parrish          org: np22                date: 2013-01-26
!
! abstract: same as intrp3oz but for special case n=1, with n argument removed.
!            This has been created to solve problem of type mismatch debug
!            compile
!            error on WCOSS.
!
! program history log:
!   2013-01-26  parrish
!
!   input argument list:
!     f        - input interpolator (gridded guess trace gas fields)
!     dx,dy,dz - input x,y,z-coords of interpolation points (grid units)
!     obstime  - observation times
!     nlevs    - number of observational layers 
!     mype     - mpi task id
!
!   output argument list:
!     g        - output interpolatees (guess trace gas at observation location)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use kinds, only: r_kind,i_kind
  use guess_grids, only: nfldsig,hrdifsig,ges_prsi
  use gridmod, only: lat2,lon2,nlat,nlon,nsig,lon1,istart,jstart
!  use constants, only: zero, one, rozcon
  use constants, only: zero, one
  use modairsco, only : pobsairsco,dpairs
  implicit none

! Declare passed variables
  integer(i_kind)                               ,intent(in   ) :: mype,nlevs
  real(r_kind)                                  ,intent(in   ) :: dx,dy,obstime
  real(r_kind),dimension(nlevs)               ,intent(in   ) :: dz
  real(r_kind),dimension(nsig+1)              ,intent(in   ) :: thpresatobs
  real(r_kind),dimension(lat2,lon2,nsig,nfldsig),intent(in   ) :: f
  real(r_kind),dimension(nlevs)                 ,intent(  out) :: g
  real(r_kind),intent(in)                                      :: conversion
  real(r_kind),dimension(nlevs,maxsub,8)      ,intent(out)   :: wdpij
  integer(i_kind),                             intent(in)    :: maxsub
  integer(i_kind),dimension(nlevs),intent(inout)              :: nsub
  integer(i_kind),dimension(nlevs,maxsub)      ,intent(out)  ::ksub

! Declare local variables
  integer(i_kind) j,k,ix,ix1,iy,iy1,kk,itsig,itsigp,iz1,iz2
  integer(i_kind) ixp,iyp,mm1
  real(r_kind) w00,w01,w10,w11,delx,dely,delx1,dely1
  real(r_kind) delz,dz1,dtsig,dtsigp,pob
  real(r_kind) delp1,delp2,delp3,delp4,delp5,delp6,delp7,delp8,delpav
  real(r_kind) delpavk(nlevs)
  real(r_kind) :: r10=10._r_kind


!*************************************************************************
! Initialize variables
  g=zero
  mm1=mype+1


! Loop over number of observations.

!    Get horizontal interpolation information.  This information includes
!    the (i,j) indices of the grid points surrounding the observation,
!    plus the corresponding interpolation weights between these points
!    and the observation.

   ix1=dx; iy1=dy
   ix1=max(1,min(ix1,nlat))
   delx=dx-ix1; dely=dy-iy1; delx=max(zero,min(delx,one))
   ix=ix1-istart(mm1)+2; iy=iy1-jstart(mm1)+2
   if(iy<1) then
      iy1=iy1+nlon
      iy=iy1-jstart(mm1)+2
   end if
   if(iy>lon1+1) then
      iy1=iy1-nlon
      iy=iy1-jstart(mm1)+2
   end if
   ixp=ix+1; iyp=iy+1
   if(ix1==nlat) then
      ixp=ix
   end if
   delx1=one-delx; dely1=one-dely
   w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely
        
!  Set the weights for linear time iterpolation from the guess to 
!  the observation time.

   if(obstime > hrdifsig(1) .and. obstime < hrdifsig(nfldsig))then
     do j=1,nfldsig-1
       if(obstime > hrdifsig(j) .and. obstime <= hrdifsig(j+1))then
         itsig=j
         itsigp=j+1
         dtsig=((hrdifsig(j+1)-obstime)/(hrdifsig(j+1)-hrdifsig(j)))
       end if
     end do
   else if(obstime <=hrdifsig(1))then
     itsig=1
     itsigp=1
     dtsig=one
   else
     itsig=nfldsig
     itsigp=nfldsig
     dtsig=one
   end if
   dtsigp=one-dtsig
     
!  Given horizontal (spatial) and temporal interpolate weights, loop 
!  over the number of layered trace gas observations at the given location

   dz1=nsig+1
!  write(6,*)mype,'nlevs',nlevs
!  call flush(6)
!   write(400+mype,*)'conversion',conversion,'nlevs',nlevs,'w00',w00,w10,w01,w11
   if(dtsigp<=0.0)then
     delpavk=0.0
     do k=1,nlevs
       pob = dz(k)
       !write(400+mype,*)'pob',k,pob,'nsub',nsub(k)
       if(pob>float(nsig+1))then
         kk=nsig
         nsub(k)=1
         ksub(k,1)=kk
         delz=dpairs(k)
         delpavk(k)=delz
         !write(400+mype,*)'above top k',k,'pob',pob,'kk',kk,'delz',delz
         g(k)=g(k)+ &
              ((f(ix ,iy ,kk,itsig )*w00 &
              + f(ixp,iy ,kk,itsig )*w10 &
              + f(ix ,iyp,kk,itsig )*w01 &
              + f(ixp,iyp,kk,itsig )*w11)*delz)
         wdpij(k,1,1)=w00*delz
         wdpij(k,1,2)=w10*delz
         wdpij(k,1,3)=w01*delz
         wdpij(k,1,4)=w11*delz
         !write(400+mype,*)'kout',k,'g',g(k),'delz',delz
       elseif(pob<1.)then
         if(k.eq.1)then
           write(6,*)'how can pob be < 1. ',pob,'k',k,dz(k)
           write(400+mype,*)'how can pob be < 1. ',pob,'k',k,'dz',dz(k)
           call flush(400+mype)
           call flush(6)
!           call stop2(333)
         elseif(dz(k-1)<1.)then
!          whole layer below sfc
           kk=1
           nsub(k)=1
           ksub(k,1)=kk
           delz=dpairs(k)
           delpavk(k)=delz
           g(k)=g(k)+ &
              ((f(ix ,iy ,kk,itsig )*w00 &
              + f(ixp,iy ,kk,itsig )*w10 &
              + f(ix ,iyp,kk,itsig )*w01 &
              + f(ixp,iyp,kk,itsig )*w11)*delz)
           wdpij(k,1,1)=w00*delz
           wdpij(k,1,2)=w10*delz
           wdpij(k,1,3)=w01*delz
           wdpij(k,1,4)=w11*delz
!           write(400+mype,*)'layer below surface stop',k,dz(k),'dzkj-1',k-1,dz(k-1),'pob',pob
!           call flush(400+mype)
!           write(6,*)'layer below surface stop',k,dz(k)
!           call flush(6)
!           call stop2(444)
         else
!          layer partly below surface
           delz=pobsairsco(k)-thpresatobs(1)*r10
           delpavk(k)=delz
           kk=1
           ksub(k,1)=kk
           nsub(k)=1
           !write(400+mype,*)'layer partly below surface ',k,'kk',kk,'delz',delz
!           write(400+mype,*)'pobsairsco',pobsairsco(k),'thpresbot',thpresatobs(1)*r10,' dp ',delz
           g(k)=g(k)+ &
             ((f(ix ,iy ,kk,itsig )*w00 &
             + f(ixp,iy ,kk,itsig )*w10 &
             + f(ix ,iyp,kk,itsig )*w01 &
             + f(ixp,iyp,kk,itsig )*w11)*delz)
           wdpij(k,1,1)=w00*delz
           wdpij(k,1,2)=w10*delz
           wdpij(k,1,3)=w01*delz
           wdpij(k,1,4)=w11*delz
           iz1=dz(k-1) 
           dz1=dz(k-1)
           do kk=1,iz1
             delz=one
             if (kk==iz1) delz=dz1-iz1
             nsub(k)=nsub(k)+1
             ksub(k,nsub(k))=kk
             !write(400+mype,*)' do part above surface delz',delz,'dz1',dz1,'iz1',iz1
             delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
             delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
             delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
             delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
             delpav=w00*delp1+w10*delp2+w01*delp3+w11*delp4
!             write(400+mype,*)'delpav above ground ',delpav,'dthpres',thpresatobs(1)-thpresatobs(2)
             delz=delz*r10 ! make mb
             delpavk(k)=delpavk(k)+delpav*delz
             g(k)=g(k) + &
                ((f(ix ,iy ,kk,itsig )*w00*delp1 &
                + f(ixp,iy ,kk,itsig )*w10*delp2 &
                + f(ix ,iyp,kk,itsig )*w01*delp3 &
                + f(ixp,iyp,kk,itsig )*w11*delp4)*delz)
             wdpij(k,kk+1,1)=w00*delp1*delz
             wdpij(k,kk+1,2)=w10*delp2*delz
             wdpij(k,kk+1,3)=w01*delp3*delz
             wdpij(k,kk+1,4)=w11*delp4*delz
           end do
         endif
       elseif(k.ne.1.and.dz(k-1)>float(nsig+1))then
!        layer partly above the top
         delz=thpresatobs(nsig+1)*r10-pobsairsco(k-1)
         delpavk(k)=delz
         kk=nsig
!         write(400+mype,*)'layer partly above top k',k,'delz',delz,'kk',kk
         g(k)=g(k)+ &
             ((f(ix ,iy ,kk,itsig )*w00 &
             + f(ixp,iy ,kk,itsig )*w10 &
             + f(ix ,iyp,kk,itsig )*w01 &
             + f(ixp,iyp,kk,itsig )*w11)*delz)
         wdpij(k,1,1)=w00*delz
         wdpij(k,1,2)=w10*delz
         wdpij(k,1,3)=w01*delz
         wdpij(k,1,4)=w11*delz
         dz1=nsig
         iz1=dz1
         iz2=pob
         nsub(k)=1
         ksub(k,1)=kk
!         write(400+mype,*)'layer part above top iz1',iz1,'iz2',iz2,'pob',pob,'k',k,'kk',kk
         do kk=iz1,iz2,-1
           delz=one
           if (kk==iz2) delz=delz-pob+iz2
           nsub(k)=nsub(k)+1
           ksub(k,nsub(k))=kk
!           write(400+mype,*)'part of top layer in raqms k',k,'kk',kk,'delz',delz
           delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
           delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
           delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
           delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
           delpav=w00*delp1+w10*delp2+w01*delp3+w11*delp4
           delz=delz*r10 ! make mb
           delpavk(k)=delpavk(k)+delpav*delz
           g(k)=g(k) + &
                ((f(ix ,iy ,kk,itsig )*w00*delp1 &
                + f(ixp,iy ,kk,itsig )*w10*delp2 &
                + f(ix ,iyp,kk,itsig )*w01*delp3 &
                + f(ixp,iyp,kk,itsig )*w11*delp4)*delz)
           wdpij(k,kk-iz2+2,1)=w00*delp1*delz
           wdpij(k,kk-iz2+2,2)=w10*delp2*delz
           wdpij(k,kk-iz2+2,3)=w01*delp3*delz
           wdpij(k,kk-iz2+2,4)=w11*delp4*delz
         end do
!         if(nsub(k)>1)then
!           write(400+mype,*)'nsub interp aa',k,nsub(k)
!           do kk=1,nsub(k)
!             write(400+mype,*)'kk ',kk,ksub(k,kk)
!           end do
!         endif
       else
         if(k.eq.1)then
           dz1=nsig+1
         else
           dz1=dz(k-1)
         endif
         iz1 = dz1
         if (iz1>nsig) iz1=nsig
         iz2 = pob
!        write(6,*)mype,'k',k,'iz1',iz1,iz2,'pob',pob
!        call flush(6)
!         write(400+mype,*)'k',k,'iz1',iz1,iz2
         nsub(k)=0
!         write(400+mype,*)' layer inside ',k,'dz1',dz1,'iz1',iz1,'pob',pob,'iz2',iz2 
         do kk=iz1,iz2,-1
           delz=one
           if (kk==iz1) delz=dz1-iz1
           if (kk==iz2) delz=delz-pob+iz2
           nsub(k)=nsub(k)+1
           ksub(k,nsub(k))=kk
!           write(400+mype,*)k,'inside nsub',nsub(k),'kk',kk
!           call flush(400+mype)
!           write(400+mype,*)'k',k,'kk',kk,'delz',delz
           delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
           delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
           delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
           delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
           delpav=w00*delp1+w10*delp2+w01*delp3+w11*delp4
           delz=delz*r10 ! make mb
           g(k)=g(k) + &
                ((f(ix ,iy ,kk,itsig )*w00*delp1 &
                + f(ixp,iy ,kk,itsig )*w10*delp2 &
                + f(ix ,iyp,kk,itsig )*w01*delp3 &
                + f(ixp,iyp,kk,itsig )*w11*delp4)*delz)
           wdpij(k,kk-iz2+1,1)=w00*delp1*delz
           wdpij(k,kk-iz2+1,2)=w10*delp2*delz
           wdpij(k,kk-iz2+1,3)=w01*delp3*delz
           wdpij(k,kk-iz2+1,4)=w11*delp4*delz
           if(wdpij(k,kk-iz2+1,1)<=0.0)then
              write(400+mype,*)'weight neg w00',w00,'delz',delz,'k',k,'iz1',iz1,iz2, &
                'dz1',dz1,'pob',pob 
              call flush(400+mype)
           endif
           !write(400+mype,*)'kk',kk,'k',k,'g(k)',g(k),'delp1',delp1,'delpav',delpav*delz*10.
           delpavk(k)=delpavk(k)+delpav*delz
!          write(6,*)mype,'g',k,g(k),'delz',delz,'kk',kk,'delp',delp1,'f',f(ix,iy,kk,itsig)
!          call flush(6)
         enddo
!         write(400+mype,*)'k',k,'g',g(k),'wdpij',wdpij(k,1,1),'delpavk',delpavk(k)*10.
         dz1=pob
       endif
!       if(nsub(k)>1)then
!         write(400+mype,*)k,'did inside bb set nsub',nsub(k)
!         do kk=1,nsub(k)
!           write(400+mype,*)kk,'ksub',ksub(k,kk)
!         end do
!       endif
     enddo
   else
     do k=1,nlevs
       pob = dz(k)
       iz1 = dz1
       if (iz1>nsig) iz1=nsig
       iz2 = pob
!        write(6,*)mype,'k',k,'iz1',iz1,iz2,'pob',pob
!        call flush(6)
        do kk=iz1,iz2,-1
           delz=one
           if (kk==iz1) delz=dz1-iz1
           if (kk==iz2) delz=delz-pob+iz2
           delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
           delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
           delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
           delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
           delp5=ges_prsi(ix ,iy ,kk,itsigp)-ges_prsi(ix ,iy ,kk+1,itsigp)
           delp6=ges_prsi(ixp,iy ,kk,itsigp)-ges_prsi(ixp,iy ,kk+1,itsigp)
           delp7=ges_prsi(ix ,iyp,kk,itsigp)-ges_prsi(ix ,iyp,kk+1,itsigp)
           delp8=ges_prsi(ixp,iyp,kk,itsigp)-ges_prsi(ixp,iyp,kk+1,itsigp)
           g(k)=g(k) + &
                ((f(ix ,iy ,kk,itsig )*w00*delp1 &
                + f(ixp,iy ,kk,itsig )*w10*delp2 &
                + f(ix ,iyp,kk,itsig )*w01*delp3 &
                + f(ixp,iyp,kk,itsig )*w11*delp4)*delz)*dtsig + &
                ((f(ix ,iy ,kk,itsigp)*w00*delp5 &
                + f(ixp,iy ,kk,itsigp)*w10*delp6 &
                + f(ix ,iyp,kk,itsigp)*w01*delp7 &
                + f(ixp,iyp,kk,itsigp)*w11*delp8)*delz)*dtsigp  
           wdpij(k,kk-iz1+1,1)=w00*delp1*delz*dtsig
           wdpij(k,kk-iz1+1,2)=w10*delp2*delz*dtsig
           wdpij(k,kk-iz1+1,3)=w01*delp3*delz*dtsig
           wdpij(k,kk-iz1+1,4)=w11*delp4*delz*dtsig
           wdpij(k,kk-iz1+1,5)=w00*delp1*delz*dtsigp
           wdpij(k,kk-iz1+1,6)=w10*delp2*delz*dtsigp
           wdpij(k,kk-iz1+1,7)=w01*delp3*delz*dtsigp
           wdpij(k,kk-iz1+1,8)=w11*delp4*delz*dtsigp
!          write(6,*)mype,'g',k,g(k),'delz',delz,'kk',kk,'delp',delp1,'f',f(ix,iy,kk,itsig)
!          call flush(6)
        enddo
        dz1=pob
      enddo
     endif
     g=g*conversion
!     do k=1,nlevs
!       if(nsub(k)>1)then
        !write(400+mype,*)'delpavk ',k,delpavk(k),' dp airs ',dpairs(k)
!        write(400+mype,*)k,' nsub ',nsub(k),' g ',g(k)
!         do kk=1,nsub(k)
!           write(400+mype,*)k,'kk',kk,'ksub',ksub(k,kk),' wdjpij ',wdpij(k,kk,1:2)
         !end do
!       endif
!     end do

! End of routine
  return
end subroutine intrp3trace1
