module stpno2mod

!$$$ module documentation block
!           .      .    .                                       .
! module:   stpno2mod    module for stpno2 and its tangent linear stpno2_tl
!  prgmmr:
!
! abstract: module for stpno2 and its tangent linear stpno2_tl
!
! program history log:
!   2005-05-17  Yanqiu zhu - wrap stpno2 and its tangent linear stpno2_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-12-02  Todling - remove stpno2_tl
!   2009-01-21  Sienkiewicz - add stpo3l (level ozone) again
!   2009-08-12  lueken - update documentation
!   2010-05-13  todling - uniform interface across stp routines
!
! subroutines included:
!   sub stpno2
!   sub stpno2lay_
!   sub stpno2lev_
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC stpno2

contains

!subroutine stpno2(no2head,o3lhead,rval,sval,out,sges,nstep)
subroutine stpno2(no2head,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stpno2       call components to calculate contrib. to
!                            penalty and stepsize for no2
!   prgmmr: sienkiewicz     org: GMAO                 date: 2009-01-22
!
! abstract: The routine calls individual components that calculate 
!           contribution to the penalty and step size from layer 
!           and level no2 measurements
!
! program history log:
!   2009-01-22  Sienkiewicz - incorporation of level ozone routine
!   2010-01-04  zhang,b - bug fix: accumulate penalty for multiple obs bins
!   2010-05-13  todling - udpate interface; gsi_bundle use
!
!   input argument list:
!     no2head
!     o3lhead
!     rno2  - search direction for no2
!     sno2  - input no2 correction field
!     sges - step size estimates (nstep)
!     nstep- number of stepsize estimates (==0 means use outer iteration value)
!
!   output argument list:
!     out(1:nstep) - contribution of no2 data to penalty sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:
!
!$$$  
  use kinds, only: r_kind,r_quad,i_kind
  use obsmod, only: no2_ob_type,o3l_ob_type
  use gridmod, only: latlon1n
  use constants, only: zero_quad,zero
  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   ) :: o3lhead
  integer(i_kind)                     ,intent(in   ) :: nstep
  type(gsi_bundle)                    ,intent(in   ) :: sval
  type(gsi_bundle)                    ,intent(in   ) :: rval
  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges
  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out

  out=zero_quad

  if(associated(no2head))call stpno2lay_(no2head, rval,sval,out,sges,nstep)
!  if(associated(o3lhead))call stpno2lev_(o3lhead,rval,sval,out,sges,nstep)

  return

end subroutine stpno2

subroutine stpno2lay_(no2head,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stpno2       compute contribution to penalty and
!                            stepsize for no2, using nonlinear qc
!   prgmmr: derber          org: np23                 date: 1995-07-11
!
! abstract: The routine computes the contribution to the penalty from no2
!           observations.  The routine also computes the contribution of
!           no2 observations to the step size.  This version includes
!           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 documenation
!   2004-07-29  treadon - add only to module use, add intent in/out
!   2004-10-07  parrish - add nonlinear qc option
!   2005-04-11  treadon - merge stpno2 and stpno2_qc into single routine
!   2005-06-14  wu      - add OMI toz
!   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
!   2006-09-18  derber  - modify output values of b1 and b3
!   2007-02-15  rancic  - add foto
!   2007-03-19  tremolet - binning of observations
!   2007-05-30  h.liu   - move interpolation weights w1-w4 inside k loop
!   2007-06-04  derber  - use quad precision to get reproducability over number of processors
!   2008-12-03  todling - update handle of foto
!   2010-05-13  todling - udpate to use gsi_bundle
!
!   input argument list:
!     no2head  - layer no2 obs type pointer to obs structure
!     rno2  - search direction for no2
!     sno2  - input no2 correction field
!     sges - step size estimates (nstep)
!     nstep- number of stepsize estimates (==0 means use outer iteration value)
!
!   output argument list:
!     out(1:nstep) - contribution of no2 data to penalty sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_quad
  use obsmod, only: no2_ob_type
  use constants, only: one,half,two,zero_quad,r3600,zero
  use gridmod, only: lat2,lon2,nsig
  use jfunc, only: l_foto,xhat_dt,dhat_dt
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use omimod, only : nlno2_omi=>nlevscatwt,no2scatterweight
  use mpimod, only : mype
  implicit none

! Declare passed variables
  type( no2_ob_type),pointer           ,intent(in   ) ::  no2head
  integer(i_kind)                     ,intent(in   ) :: nstep
  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
  type(gsi_bundle)                    ,intent(in   ) :: rval,sval
  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges

! Declare local variables
  integer(i_kind) i,j,ij,ier,istatus
  integer(i_kind) k,j1,j2,j3,j4,iz1,iz2,j1x,j2x,j3x,j4x,kk,kl
  real(r_kind) dz1,pob,delz
  real(r_kind) w1,w2,w3,w4,time_no2,no2
  real(r_kind),dimension(max(1,nstep))::pen
  real(r_kind),pointer,dimension(:) :: xhat_dt_no2
  real(r_kind),pointer,dimension(:) :: dhat_dt_no2
  real(r_kind),allocatable,dimension(:,:) :: rno2,sno2
  real(r_kind),pointer,dimension(:,:,:)   :: rno2p,sno2p
  type( no2_ob_type), pointer ::  no2ptr

  real(r_kind) :: val,val1
  real(r_kind),dimension(nlno2_omi) :: val_lay,val_lay1
!  logical no2scatterweight
!  data no2scatterweight/.true./

! Get pointers and return if not found
  ier=0
!  write(700+mype,*)'top stpno2'
  call gsi_bundlegetpointer(sval,'no2',sno2p,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'no2',rno2p,istatus);ier=istatus+ier
  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
!  write(700+mype,*)'ier',ier

  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
     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
!  if(nstep.eq.0)then
!    write(700+mype,*)'stpno2 NSTEP ZERO AAAAAAA'
!    call flush(700+mype)
!  endif
  do while (associated(no2ptr))
!     write(700+mype,*)'luse',no2ptr%luse

     if(no2ptr%luse)then

        if(nstep > 0)then
!          Get location
           j1=no2ptr%ij(1)
           j2=no2ptr%ij(2)
           j3=no2ptr%ij(3)
           j4=no2ptr%ij(4)
           if(l_foto)time_no2=no2ptr%time*r3600

!          Accumulate contribution from layer observations
           dz1=nsig+1
        end if

        if ( no2ptr%nlno2 >= 1 ) then

           do k=1,no2ptr%nlno2
              if(nstep > 0)then
                 val1= -no2ptr%res(k)
                 val = zero_quad
                 pob = no2ptr%prs(k)
                 iz1 = dz1
                 if (iz1>nsig) iz1=nsig
                 iz2 = pob

                 do kk=iz2,iz1
                    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)
                    val=val + ( &
                         w1* rno2(j1,kk)+ &
                         w2* rno2(j2,kk)+ &
                         w3* rno2(j3,kk)+ &
                         w4* rno2(j4,kk))*delz
                    val1=val1 + ( &
                         w1* sno2(j1,kk)+ &
                         w2* sno2(j2,kk)+ &
                         w3* sno2(j3,kk)+ &
                         w4* sno2(j4,kk))*delz
                    if(l_foto) then
                       j1x=j1+(kk-1)*lat2*lon2
                       j2x=j2+(kk-1)*lat2*lon2
                       j3x=j3+(kk-1)*lat2*lon2
                       j4x=j4+(kk-1)*lat2*lon2
                       val=val + ( &
                         (w1*dhat_dt_no2(j1x)+ &
                          w2*dhat_dt_no2(j2x)+ &
                          w3*dhat_dt_no2(j3x)+ &
                          w4*dhat_dt_no2(j4x))*time_no2 )*delz
                       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
                    end if
                 end do
                 do kk=1,nstep
                    no2=val1+sges(kk)*val
                    pen(kk)= no2ptr%err2(k)*no2*no2
                 end do
              else
                 pen(1)=no2ptr%res(k)*no2ptr%res(k)*no2ptr%err2(k)
              end if

              out(1) = out(1)+pen(1)*no2ptr%raterr2(k)
              do kk=2,nstep
                 out(kk) = out(kk)+(pen(kk)-pen(1))*no2ptr%raterr2(k)
              end do
              dz1=pob
           end do
           
        end if

!       Add contribution from total column observation
!        write(700+mype,*)'nstep',nstep
        if(nstep > 0)then
           k   = no2ptr%nlno2+1
           val1= -no2ptr%res(k)
           val  = zero
!           write(700+mype,*)'k',k,'val1',val1
!           call flush(700+mype)
!           do kk=1,nsig
!           if(allocated(no2ptr%avgkernel))then
           if(associated(no2ptr%avgkernel))then
             dz1=no2ptr%prs(1)
             do kl=1,no2ptr%nltrop
               val_lay(kl) = zero
               val_lay1(kl)= zero
               pob = no2ptr%prs(kl+1)
               iz1=dz1
               if (iz1 > nsig) iz1=nsig
               iz2=pob
!               write(700+mype,*)'kl',kl,'iz1',iz1,iz2,'pob',pob
!               call flush(700+mype)
               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)
!                  write(700+mype,*)'kk',kk,'delz',delz,'w1',w1
                  val_lay(kl)=val_lay(kl) + ( & 
                       w1* rno2(j1,kk)+ &
                       w2* rno2(j2,kk)+ &
                       w3* rno2(j3,kk)+ &
                       w4* rno2(j4,kk))*delz   
                  val_lay1(kl)=val_lay1(kl) + ( & 
                       w1* sno2(j1,kk)+ &
                       w2* sno2(j2,kk)+ &
                       w3* sno2(j3,kk)+ &
                       w4* sno2(j4,kk))*delz    
!                 write(750+mype,*)'rno2 ',kk,rno2(j1,kk)
!                 write(500+mype,*)'stpno2 val_lay',kl,val_lay(kl),val_lay1(kl)
               enddo
               dz1=pob 
!               write(700+mype,*)'stpno2 ',kl,'val_lay',val_lay(kl),'scwt',no2ptr%scatterweight(kl)
!               write(700+mype,*)'val_lay1',val_lay1(kl)
               val_lay(kl)=val_lay(kl)*no2ptr%avgkernel(kl)
               val_lay1(kl)=val_lay1(kl)*no2ptr%avgkernel(kl)
               val=val+val_lay(kl)
               val1=val1+val_lay1(kl)
!               write(500+mype,*)'stpno2 val',val,val1,'delz',delz
             enddo
           elseif(no2scatterweight)then
             dz1=no2ptr%prs(1)
             do kl=1,no2ptr%nltrop
               val_lay(kl) = zero
               val_lay1(kl)= zero
               pob = no2ptr%prs(kl+1)
               iz1=dz1
               if (iz1 > nsig) iz1=nsig
               iz2=pob
!               write(700+mype,*)'kl',kl,'iz1',iz1,iz2,'pob',pob
!               call flush(700+mype)
               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)
!                  write(700+mype,*)'kk',kk,'delz',delz,'w1',w1
                  val_lay(kl)=val_lay(kl) + ( & 
                       w1* rno2(j1,kk)+ &
                       w2* rno2(j2,kk)+ &
                       w3* rno2(j3,kk)+ &
                       w4* rno2(j4,kk))*delz   
                  val_lay1(kl)=val_lay1(kl) + ( & 
                       w1* sno2(j1,kk)+ &
                       w2* sno2(j2,kk)+ &
                       w3* sno2(j3,kk)+ &
                       w4* sno2(j4,kk))*delz    
!                 write(750+mype,*)'rno2 ',kk,rno2(j1,kk)
!                 write(700+mype,*)'val_lay',kl,val_lay(kl),val_lay1(kl)
               enddo
               dz1=pob 
!               write(700+mype,*)'stpno2 ',kl,'val_lay',val_lay(kl),'scwt',no2ptr%scatterweight(kl)
!               write(700+mype,*)'val_lay1',val_lay1(kl)
               val_lay(kl)=val_lay(kl)*no2ptr%scatterweight(kl)
               val_lay1(kl)=val_lay1(kl)*no2ptr%scatterweight(kl)
               val=val+val_lay(kl)
               val1=val1+val_lay1(kl)
!               write(700+mype,*)'val',val,val1,'delz',delz
             enddo
           else
!             write(700+mype,*)'nltrop',no2ptr%nltrop
             do kk=1,no2ptr%nltrop
                w1=no2ptr%wij(1,kk)
                w2=no2ptr%wij(2,kk)
                w3=no2ptr%wij(3,kk)
                w4=no2ptr%wij(4,kk)
                val=val+  (          &
                     w1* rno2(j1,kk)+ &
                     w2* rno2(j2,kk)+ &
                     w3* rno2(j3,kk)+ &
                     w4* rno2(j4,kk))
!                write(750+mype,*)'rno2',kk,rno2(j1,kk)
                val1=val1 +  (       &
                     w1* sno2(j1,kk)+ &
                     w2* sno2(j2,kk)+ &
                     w3* sno2(j3,kk)+ & 
                     w4* sno2(j4,kk))
!                write(700+mype,*)'val ',val,' val1 ',val1
                if(l_foto)then
                   j1x=j1+(kk-1)*lat2*lon2
                   j2x=j2+(kk-1)*lat2*lon2
                   j3x=j3+(kk-1)*lat2*lon2
                   j4x=j4+(kk-1)*lat2*lon2
                   val=val+ ( &
                     (w1*xhat_dt_no2(j1x)+ &
                      w2*xhat_dt_no2(j2x)+ &
                      w3*xhat_dt_no2(j3x)+ & 
                      w4*xhat_dt_no2(j4x))*time_no2 )
                   val1=val1 + ( &
                     (w1*dhat_dt_no2(j1x)+ &
                      w2*dhat_dt_no2(j2x)+ &
                      w3*dhat_dt_no2(j3x)+ &
                      w4*dhat_dt_no2(j4x))*time_no2 )
                end if
!                write(500+mype,*)'val',val,val1
             enddo
             
           endif
           do kk=1,nstep
              no2=val1+sges(kk)*val
!              write(700+mype,*)'kk',kk,'no2',no2,'sges',sges(kk)
!              write(750+mype,*)kk,'no2',no2,'sges',sges(kk),'val1',val1,val,'err',no2ptr%err2(k)
              pen(kk)= no2ptr%err2(k)*no2*no2
           end do
        else
!           pen(kk)=no2ptr%res(k)*no2ptr%res(k)*no2ptr%err2(k)
!          fix ajl 3/24/2017 had fixed in ozone before
!           write(700+mype,*)'USE stpno2 Pen1 fix AAAA'
           pen(1)=no2ptr%res(k)*no2ptr%res(k)*no2ptr%err2(k)
        end if

        out(1) = out(1) +pen(1)*no2ptr%raterr2(k)
!        write(750+mype,*)'pen',pen(1:nstep)
        do kk=2,nstep
           out(kk) = out(kk) +(pen(kk)-pen(1))*no2ptr%raterr2(k)
!           write(750+mype,*)'dpen',kk,pen(kk)-pen(1)
        end do
     end if

     no2ptr => no2ptr%llpoint

! End of loop over observations
  enddo
!  write(700+mype,*)'final stpno2lay out',out(1:nstep)
!  write(500+mype,*)'final stpno2lay out',out(1:nstep)
!  call flush(500+mype)

! Clean 
  deallocate(sno2,rno2)

! End of routine.
  return
end subroutine stpno2lay_

subroutine stpno2lev_(o3lhead,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stpno2lev    compute contribution to penalty and
!                            stepsize for o3 level obs, using nonlinear qc
!   prgmmr: sienkiewicz     org: GMAO                 date: 2006-09-14
!
! abstract: The routine computes the contribution to the penalty from ozone
!           observations.  The routine also computes the contribution of
!           ozone observations to the step size.  This version includes
!           nonlinear qc.
!
! program history log:
!   2006-09-14  sienkiewicz - add level ozone obs
!   2007-01-02  sienkiewicz - separate subroutine
!   2007-01-05  sienkiewicz - update to 9/2006 GSI (new obs structure)
!   2009-01-21  sienkiewicz - update to 1/2009 GSI, changes based on stpq & stpno2
!   2010-01-04  zhang,b - bug fix: accumulate penalty for multiple obs bins
!   2010-05-13  todling - udpate to use gsi_bundle
!
!   input argument list:
!     o3lhead - level ozone obs type pointer to obs structure
!     rno21d  - search direction for ozone (as 1d var)
!     sno21d  - input ozone correction field (as 1d var)
!     sges - step size estimates (4)
!     nstep- number of stepsizes (==0 means use outer iteration values)
!
!   output argument list:
!     out(1:nstep) - contribution of ozone data to penalty sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_quad
  use obsmod, only: o3l_ob_type
  use constants, only: zero,one,half,two,r3600
  use gridmod, only: latlon1n
  use jfunc, only: l_foto,xhat_dt,dhat_dt
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  implicit none

! Declare passed variables
  type(o3l_ob_type),pointer           ,intent(in   ) :: o3lhead
  integer(i_kind)                     ,intent(in   ) :: nstep
  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
  type(gsi_bundle)                    ,intent(in   ) :: rval,sval
  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges

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

! Get pointers and return if 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

! Initialize output variables to zero

  time_no2 = zero

  o3lptr => o3lhead

! Loop over level ozone observations
!
  do while (associated(o3lptr))
     if(o3lptr%luse)then
        if(nstep > 0)then
           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)


           val= w1*rno21d(j1)+w2*rno21d(j2)+w3*rno21d(j3)+w4*rno21d(j4)+ &
                w5*rno21d(j5)+w6*rno21d(j6)+w7*rno21d(j7)+w8*rno21d(j8)   
           val2=w1*sno21d(j1)+w2*sno21d(j2)+w3*sno21d(j3)+w4*sno21d(j4)+ &
                w5*sno21d(j5)+w6*sno21d(j6)+w7*sno21d(j7)+w8*sno21d(j8)-o3lptr%res

           if(l_foto) then
              time_no2=o3lptr%time*r3600
              val=val+ (w1*dhat_dt_no2(j1)+w2*dhat_dt_no2(j2)+ &
                        w3*dhat_dt_no2(j3)+w4*dhat_dt_no2(j4)+ &
                        w5*dhat_dt_no2(j5)+w6*dhat_dt_no2(j6)+ &
                        w7*dhat_dt_no2(j7)+w8*dhat_dt_no2(j8))*time_no2
              val2=val2+ (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_no2
           end if

           do kk=1,nstep
              oz=val2+sges(kk)*val
              pen(kk)= oz*oz*o3lptr%err2
           end do
        else
           pen(1) = o3lptr%res*o3lptr%res*o3lptr%err2
        end if

        out(1) = out(1)+pen(1)*o3lptr%raterr2
        do kk=2,nstep
           out(kk) = out(kk)+(pen(kk)-pen(1))*o3lptr%raterr2
        end do
     end if

     o3lptr => o3lptr%llpoint

  end do

! End of routine.
  return
end subroutine stpno2lev_

end module stpno2mod
