module intno2mod

!$$$ module documentation block
!           .      .    .                                       .
! module:   intno2mod    module for intno2 and its tangent linear intno2_tl
!   prgmmr:
!
! abstract: module for intno2 and its tangent linear intno2_tl
!
! program history log:
!   2005-05-13  Yanqiu zhu - wrap intno2 and its tangent linear intno2_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-11-26  Todling - remove intno2_tl; add interface back
!   2009-08-13  lueken - update documentation
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - implemented obs adjoint test  
!
! subroutines included:
!   sub intno2_
!   sub intno2lay_
!   sub intno2lev_
!
! variable definitions:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC intno2

interface intno2; module procedure &
          intno2_
end interface

contains

!subroutine intno2_(no2head,no2lhead,rval,sval)
subroutine intno2_(no2head,rval,sval)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intno2       call individual no2 obs operators
!   prgmmr: todling       org: np23                date: 2008-11-28
!
! abstract:  This routine calls the individual components of the 
!            no2 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
!
!   input argument list:
!     no2head  - layer no2 obs type pointer to obs structure
!     no2lhead - level no2 obs type pointer to obs structure
!     sno2     - no2 increment in grid space
!
!   output argument list:
!     rno2    - no2 results from observation operator 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use obsmod, only: no2_ob_type,o3l_ob_type
  use gsi_bundlemod, only: gsi_bundle
  implicit none

! Declare passed variables
  type( no2_ob_type),pointer,intent(in   ) :: no2head
!  type(o3l_ob_type),pointer,intent(in   ) :: no2lhead
  type(gsi_bundle),intent(in   ) :: sval
  type(gsi_bundle),intent(inout) :: rval

!  If obs exist call int routines
  if(associated(no2head))call intno2lay_( no2head,rval,sval)
!  if(associated(no2lhead))call intno2lev_(no2lhead,rval,sval)

end subroutine intno2_

subroutine intno2lay_(no2head,rval,sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intno2       apply nonlin qc obs operator for no2 
!   prgmmr: derber           org: np23                date: 1995-07-11
!
! abstract:  This routine applies the observation operator (forward
!            model) and adjoint of this operator for no2 observations
!            with the addition of nonlinear qc.
!
! program history log:
!   1995-07-11  derber
!   1999-03-01  wu - port cray90 code to ibm-sp (mpi version)
!   2004-06-16  treadon - update documentation
!   2004-08-02  treadon - add only to module use, add intent in/out
!   2004-10-08  parrish - add nonlinear qc
!   2005-03-01  parrish - nonlinear qc change to account for inflated obs error
!   2005-04-11  treadon - merge intno2 and intno2_qc into single routine
!   2005-06-14  wu      - add OMI total ozone
!   2005-09-28  derber  - consolidate location and weight arrays
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!                       - unify NL qc
!   2007-02-15  rancic - add foto
!   2007-02-16  sienkiewicz - add call to routine for level ozone contrib.
!   2007-03-19  tremolet - binning of observations
!   2007-05-30  h.liu   - move interpolation weights w1-w4 inside k loop
!   2007-06-05  tremolet - use observation diagnostics structure
!   2007-07-09  tremolet - observation sensitivity
!   2008-01-04  tremolet - Don't apply H^T if l_do_adjoint is false
!   2008-??-??  ??????   - remove nonlinear qc gradient; folded OMI within layer O3
!   2008-11-28  todling  - turn FOTO optional; changed ptr%time handle
!   2009-01-18  todling  - treat val in quad precision (revisit later)
!   2010-05-13  todling  - update to use gsi_bundle; update interface
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - introduced ladtest_obs         
!
!   input argument list:
!     no2head  - layer no2 obs type pointer to obs structure
!     sno2     - no2 increment in grid space
!     rno2
!
!   output argument list:
!     rno2    - no2 results from observation operator 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use kinds, only: r_kind,i_kind,r_quad
  use obsmod, only: no2_ob_type,lsaveobsens,l_do_adjoint
  use gridmod, only: lat2,lon2,nsig
  use jfunc, only: jiter,l_foto,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 ! ajl
  implicit none

! Declare passed variables
  type( no2_ob_type),pointer,intent(in   ) :: no2head
  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,j1,j2,j3,j4,kk,iz1,iz2,j1x,j2x,j3x,j4x
  real(r_kind) dz1,pob,delz,time_no2
  real(r_quad) val1,valx
  real(r_kind) w1,w2,w3,w4
  real(r_kind),pointer,dimension(:) :: xhat_dt_no2
  real(r_kind),pointer,dimension(:) :: dhat_dt_no2
  real(r_kind),pointer,dimension(:,:,:)  :: sno2p
  real(r_kind),pointer,dimension(:,:,:)  :: rno2p
  real(r_kind),allocatable,dimension(:,:) :: sno2
  real(r_kind),allocatable,dimension(:,:) :: rno2
  real(r_kind),dimension(nsig) :: sno2max,sno2min,rno2max,rno2min ! ajl
  integer(i_kind) :: ntropmax ! ajl
  type(no2_ob_type), pointer :: no2ptr
  ntropmax=1 ! ajl

! Retrieve pointers
! Simply return if any pointer not found
  ier=0 
!  write(860+mype,*)'top intno2'
!  call flush(860+mype)
  call gsi_bundlegetpointer(sval,'no2',sno2p,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'no2',rno2p,istatus);ier=istatus+ier
!  write(860+mype,*)nsig,'sno2p',maxval(sno2p),minval(sno2p)
!  write(860+mype,*)nsig,'rno2p',maxval(rno2p),minval(rno2p)
  if(l_foto) then
     call gsi_bundlegetpointer(xhat_dt,'no2',xhat_dt_no2,istatus);ier=istatus+ier
     call gsi_bundlegetpointer(dhat_dt,'no2',dhat_dt_no2,istatus);ier=istatus+ier
  endif
  if(ier/=0)return

! Can't do rank-2 pointer into rank-2, therefore, allocate work space
  allocate(sno2(lat2*lon2,nsig),rno2(lat2*lon2,nsig))
  do k=1,nsig
!    sno2max(k)=maxval(sno2p(:,:,k))
!    sno2min(k)=minval(sno2p(:,:,k))
!    rno2max(k)=maxval(rno2p(:,:,k))
!    rno2min(k)=minval(rno2p(:,:,k))
!    write(860+mype,*)k,'sno2max',sno2max(k),sno2min(k)
!    write(860+mype,*)rno2max(k),rno2min(k)
     ij=0
     do j=1,lon2
        do i=1,lat2
           ij=ij+1
           sno2(ij,k) = sno2p(i,j,k)
           rno2(ij,k) = rno2p(i,j,k)
        enddo
     enddo
  enddo
!
! SBUV OZONE: LAYER O3 and TOTAL O3
!
! Loop over no2 observations.
  no2ptr => no2head
  do while (associated(no2ptr))

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


!    Accumulate contribution from layer observations
     dz1=nsig+1
!     write(550+mype,*)'nlno2',no2ptr%nlno2
     if ( no2ptr%nlno2 >= 1 ) then

        if(l_foto) time_no2 = no2ptr%time*r3600
        do k=1,no2ptr%nlno2
           val1= zero_quad
           pob = no2ptr%prs(k)
           iz1=dz1
           if (iz1 > nsig) iz1=nsig
           iz2=pob
           do kk=iz1,iz2,-1
              delz=one
              if (kk==iz1) delz=dz1-iz1
              if (kk==iz2) delz=delz-pob+iz2
              w1=no2ptr%wij(1,kk)
              w2=no2ptr%wij(2,kk)
              w3=no2ptr%wij(3,kk)
              w4=no2ptr%wij(4,kk)
              val1=val1 + ( &
                   w1* sno2(j1,kk)+ &
                   w2* sno2(j2,kk)+ &
                   w3* sno2(j3,kk)+ &
                   w4* sno2(j4,kk))*delz
              if (l_foto) then
                 j1x=w1+(kk-1)*lat2*lon2
                 j2x=w2+(kk-1)*lat2*lon2
                 j3x=w3+(kk-1)*lat2*lon2
                 j4x=w4+(kk-1)*lat2*lon2
                 val1=val1 + ( &
                     (w1*xhat_dt_no2(j1x)+ &
                      w2*xhat_dt_no2(j2x)+ &
                      w3*xhat_dt_no2(j3x)+ &
                      w4*xhat_dt_no2(j4x))*time_no2)*delz
              endif
           enddo

           if (lsaveobsens) then
              no2ptr%diags(k)%ptr%obssen(jiter)=val1*no2ptr%err2(k)*no2ptr%raterr2(k)
           else
              if (no2ptr%luse) no2ptr%diags(k)%ptr%tldepart(jiter)=val1
           endif

           if (l_do_adjoint) then
              if (lsaveobsens) then
                 valx = no2ptr%diags(k)%ptr%obssen(jiter)

              else
                 if(ladtest_obs) then
                    valx     = val1
                 else
                    val1=val1-no2ptr%res(k)

                    valx     = val1*no2ptr%err2(k) 
                    valx     = valx*no2ptr%raterr2(k)
                 end if
              endif

              do kk=iz1,iz2,-1
                 delz=one
                 if(kk==iz1)delz=dz1-iz1
                 if(kk==iz2)delz=delz-pob+iz2
                 w1=no2ptr%wij(1,kk)
                 w2=no2ptr%wij(2,kk)
                 w3=no2ptr%wij(3,kk)
                 w4=no2ptr%wij(4,kk)
                 rno2(j1,kk)  =  rno2(j1,kk) + valx*w1*delz
                 rno2(j2,kk)  =  rno2(j2,kk) + valx*w2*delz
                 rno2(j3,kk)  =  rno2(j3,kk) + valx*w3*delz
                 rno2(j4,kk)  =  rno2(j4,kk) + valx*w4*delz
              enddo
              if (l_foto) then
                 do kk=iz1,iz2,-1
                    delz=one
                    if(kk==iz1)delz=dz1-iz1
                    if(kk==iz2)delz=delz-pob+iz2
                    w1=no2ptr%wij(1,kk)
                    w2=no2ptr%wij(2,kk)
                    w3=no2ptr%wij(3,kk)
                    w4=no2ptr%wij(4,kk)
                    j1x=w1+(kk-1)*lat2*lon2
                    j2x=w2+(kk-1)*lat2*lon2
                    j3x=w3+(kk-1)*lat2*lon2
                    j4x=w4+(kk-1)*lat2*lon2
                    dhat_dt_no2(j1x) = dhat_dt_no2(j1x) + valx*w1*delz*time_no2
                    dhat_dt_no2(j2x) = dhat_dt_no2(j2x) + valx*w2*delz*time_no2
                    dhat_dt_no2(j3x) = dhat_dt_no2(j3x) + valx*w3*delz*time_no2
                    dhat_dt_no2(j4x) = dhat_dt_no2(j4x) + valx*w4*delz*time_no2
                 enddo
              endif
              dz1=pob
           endif
        end do

     end if   ! (no2ptr%nlno2 >= 1)

!    Add contribution from total column observation
     k=no2ptr%nlno2+1
     val1= zero
!     do kk=nsig,1,-1
     ntropmax=max(no2ptr%nltrop,ntropmax)
     do kk=no2ptr%nltrop,1,-1
        w1=no2ptr%wij(1,kk)
        w2=no2ptr%wij(2,kk)
        w3=no2ptr%wij(3,kk)
        w4=no2ptr%wij(4,kk)
        val1=val1 + &
             w1* sno2(j1,kk)+ &
             w2* sno2(j2,kk)+ &
             w3* sno2(j3,kk)+ &
             w4* sno2(j4,kk)
     enddo
!     write(550+mype,*)'val1',val1,'nsig',nsig
     if (l_foto) then
        do kk=nsig,1,-1
           w1=no2ptr%wij(1,kk)
           w2=no2ptr%wij(2,kk)
           w3=no2ptr%wij(3,kk)
           w4=no2ptr%wij(4,kk)
           j1x=w1+(kk-1)*lat2*lon2
           j2x=w2+(kk-1)*lat2*lon2
           j3x=w3+(kk-1)*lat2*lon2
           j4x=w4+(kk-1)*lat2*lon2
           val1=val1 + &
               (w1*xhat_dt_no2(j1x)+ &
                w2*xhat_dt_no2(j2x)+ &
                w3*xhat_dt_no2(j3x)+ &
                w4*xhat_dt_no2(j4x))*time_no2
        enddo
     endif

     if (lsaveobsens) then
        no2ptr%diags(k)%ptr%obssen(jiter)=val1*no2ptr%err2(k)*no2ptr%raterr2(k)
     else
        if (no2ptr%luse) no2ptr%diags(k)%ptr%tldepart(jiter)=val1
     endif

     if (l_do_adjoint) then
        if (lsaveobsens) then
           valx = no2ptr%diags(k)%ptr%obssen(jiter)

        else
           if(ladtest_obs) then
              valx     = val1
           else
              val1=val1-no2ptr%res(k)

              valx     = val1*no2ptr%err2(k)
              valx     = valx*no2ptr%raterr2(k)
           end if

        endif

!        do kk=nsig,1,-1
        do kk=no2ptr%nltrop,1,-1
           w1=no2ptr%wij(1,kk)
           w2=no2ptr%wij(2,kk)
           w3=no2ptr%wij(3,kk)
           w4=no2ptr%wij(4,kk)
           rno2(j1,kk)  = rno2(j1,kk) + valx*w1
           rno2(j2,kk)  = rno2(j2,kk) + valx*w2
           rno2(j3,kk)  = rno2(j3,kk) + valx*w3
           rno2(j4,kk)  = rno2(j4,kk) + valx*w4
        enddo
        if (l_foto) then
           do kk=nsig,1,-1
              w1=no2ptr%wij(1,kk)
              w2=no2ptr%wij(2,kk)
              w3=no2ptr%wij(3,kk)
              w4=no2ptr%wij(4,kk)
              j1x=w1+(kk-1)*lat2*lon2
              j2x=w2+(kk-1)*lat2*lon2
              j3x=w3+(kk-1)*lat2*lon2
              j4x=w4+(kk-1)*lat2*lon2
              dhat_dt_no2(j1x) =dhat_dt_no2(j1x) + valx*w1*time_no2
              dhat_dt_no2(j2x) =dhat_dt_no2(j2x) + valx*w2*time_no2
              dhat_dt_no2(j3x) =dhat_dt_no2(j3x) + valx*w3*time_no2
              dhat_dt_no2(j4x) =dhat_dt_no2(j4x) + valx*w4*time_no2
           enddo
        endif
     endif

     no2ptr => no2ptr%llpoint

! End loop over observations
  enddo
!     write(550+mype,*)'rno2',maxval(rno2),minval(rno2)

! Copy output and clean up 
  do k=1,nsig
     ij=0
     do j=1,lon2
        do i=1,lat2
           ij=ij+1
           rno2p(i,j,k) = rno2(ij,k)
        enddo
     enddo
!    rno2max(k)=maxval(rno2p(:,:,k))
!    rno2min(k)=minval(rno2p(:,:,k))
!    write(860+mype,*)'rno2end',k,rno2max(k),rno2min(k)
  enddo
  deallocate(sno2,rno2)
!  write(860+mype,*)'ntropmax',ntropmax

! End of routine
  return
end subroutine intno2lay_

subroutine intno2lev_(no2lhead,rval,sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    into3l       apply nonlin qc obs operator for o3 level
!   prgmmr: sienkiewicz      org: GMAO                date: 2006-09-14
!
! abstract:  This routine applies the observation operator (forward
!            model) and adjoint of this operator for ozone level 
!            observations with the addition of nonlinear qc.
!
! to do: add time derivatives correctly (Todling) 
!
! program history log:
!   2006-09-14  sienkiewicz - add level ozone
!   2007-01-02  sienkiewicz - separate from intno2 (again)
!   2007-02-16  sienkiewicz - changes for new inner loop obs data structure
!   2009-01-08  todling - remove nonlinear qc
!   2009-01-22  sienkiewicz - add time derivative
!   2010-05-13  todling  - update to use gsi_bundle; update interface
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - introduced ladtest_obs         
!
!   input argument list:
!     no2lhead - level ozone obs type pointer to obs structure
!     sno21d   - ozone increment in grid space
!     rno21d
!
!   output argument list:
!     rno21d   - results from observation operator (0 for no data)
!
! attributes:
!   language: f90
!   machine:  -
!
!$$$
!--------

  use kinds, only: r_kind,i_kind
  use obsmod, only: o3l_ob_type,lsaveobsens, l_do_adjoint
  use gridmod, only: latlon1n
  use constants, only: r3600
  use jfunc, only: jiter,l_foto,xhat_dt,dhat_dt
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_4dvar, only: ladtest_obs
  implicit none

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

! Declare local variables
  integer(i_kind) ier,istatus
  integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8
  real(r_kind) val,grad
  real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8,time_o3l
  real(r_kind),pointer,dimension(:) :: xhat_dt_no2
  real(r_kind),pointer,dimension(:) :: dhat_dt_no2
  real(r_kind),pointer,dimension(:) :: sno21d
  real(r_kind),pointer,dimension(:) :: rno21d
  type(o3l_ob_type), pointer :: o3lptr

! Retrieve pointers
! Simply return if any pointer not found
  ier=0 
  call gsi_bundlegetpointer(sval,'oz',sno21d,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'oz',rno21d,istatus);ier=istatus+ier
  if(l_foto) then
     call gsi_bundlegetpointer(xhat_dt,'oz',xhat_dt_no2,istatus);ier=istatus+ier
     call gsi_bundlegetpointer(dhat_dt,'oz',dhat_dt_no2,istatus);ier=istatus+ier
  endif
  if(ier/=0)return

! LEVEL-OZONE OBSERVATIONS

! Loop over ozone observations.


  o3lptr => no2lhead

  do while (associated(o3lptr))
     j1=o3lptr%ij(1)
     j2=o3lptr%ij(2)
     j3=o3lptr%ij(3)
     j4=o3lptr%ij(4)
     j5=o3lptr%ij(5)
     j6=o3lptr%ij(6)
     j7=o3lptr%ij(7)
     j8=o3lptr%ij(8)
     w1=o3lptr%wij(1)
     w2=o3lptr%wij(2)
     w3=o3lptr%wij(3)
     w4=o3lptr%wij(4)
     w5=o3lptr%wij(5)
     w6=o3lptr%wij(6)
     w7=o3lptr%wij(7)
     w8=o3lptr%wij(8)


!    Forward model
     val=w1*sno21d(j1)+w2*sno21d(j2)+w3*sno21d(j3)+w4*sno21d(j4)+ &
          w5*sno21d(j5)+w6*sno21d(j6)+w7*sno21d(j7)+w8*sno21d(j8)

     if ( l_foto ) then
        time_o3l=o3lptr%time*r3600
        val=val+&
            (w1*xhat_dt_no2(j1)+w2*xhat_dt_no2(j2)+ &
             w3*xhat_dt_no2(j3)+w4*xhat_dt_no2(j4)+ &
             w5*xhat_dt_no2(j5)+w6*xhat_dt_no2(j6)+ &
             w7*xhat_dt_no2(j7)+w8*xhat_dt_no2(j8))*time_o3l
     endif

     if (lsaveobsens) then
        o3lptr%diags%obssen(jiter) = val*o3lptr%raterr2*o3lptr%err2
     else
        if (o3lptr%luse) o3lptr%diags%tldepart(jiter)=val
     endif

     if (l_do_adjoint) then
        if (lsaveobsens) then
           grad = o3lptr%diags%obssen(jiter)

        else
           if( ladtest_obs ) then
              grad = val
           else
              val=val-o3lptr%res

              grad = val*o3lptr%raterr2*o3lptr%err2
           end if
        endif

!    Adjoint
        rno21d(j1)=rno21d(j1)+w1*grad
        rno21d(j2)=rno21d(j2)+w2*grad
        rno21d(j3)=rno21d(j3)+w3*grad
        rno21d(j4)=rno21d(j4)+w4*grad
        rno21d(j5)=rno21d(j5)+w5*grad
        rno21d(j6)=rno21d(j6)+w6*grad
        rno21d(j7)=rno21d(j7)+w7*grad
        rno21d(j8)=rno21d(j8)+w8*grad

        if ( l_foto ) then
           grad=grad*time_o3l
           dhat_dt_no2(j1)=dhat_dt_no2(j1)+w1*grad
           dhat_dt_no2(j2)=dhat_dt_no2(j2)+w2*grad
           dhat_dt_no2(j3)=dhat_dt_no2(j3)+w3*grad
           dhat_dt_no2(j4)=dhat_dt_no2(j4)+w4*grad
           dhat_dt_no2(j5)=dhat_dt_no2(j5)+w5*grad
           dhat_dt_no2(j6)=dhat_dt_no2(j6)+w6*grad
           dhat_dt_no2(j7)=dhat_dt_no2(j7)+w7*grad
           dhat_dt_no2(j8)=dhat_dt_no2(j8)+w8*grad
        endif

     endif

     o3lptr => o3lptr%llpoint

  end do

! End of routine
  return

end subroutine intno2lev_

end module intno2mod
