SUBROUTINE ozon2d(imt,lm,pozn,xlat,rsin1,rcos1,rcos2,qo3)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! subprogram:    ozon2d      compute ozone in model columns
  !   prgrmmr: katz/campana    org: w/np2      date: ??-??-??
  ! abstract:
  !     ozon2d computes the ozone mixing ratio in each grid box
  !     of columns within the model domain
  ! program history log:
  !   ??-??-??  katz/kc    - originator
  !   96-07-26  black      - modified for eta model
  ! usage: call ozon2d from subroutine radtn
  !   input argument list:
  !     imt:   number of columns in the slab
  !     lm:    number of layers in columns
  !     xlat:  geodetic latitude of grid columns in degrees
  !     rsin1: information relating to position of earth in its orbit
  !     rcos1: information relating to position of earth in its orbit
  !     rcos2: information relating to position of earth in its orbit
  !   output argument list:
  !     qo3:  ozone mixing ratio at midlayers of model columns

  IMPLICIT NONE
  INTEGER, PARAMETER :: nl=81,nlp1=nl+1,lngth=37*nl
  REAL, PARAMETER :: rtd=57.2957795

  INTEGER, INTENT(in) :: imt,lm
  REAL, INTENT(in) :: rsin1,rcos1,rcos2
  REAL, DIMENSION(imt,lm), INTENT (out) :: qo3
  REAL, DIMENSION(imt,lm),INTENT(in) ::pozn
  REAL, DIMENSION(imt),INTENT(in)   ::  xlat

  REAL, DIMENSION(imt)   :: tthan
  INTEGER,DIMENSION(imt) :: jjrow
  REAL, DIMENSION(imt,nl) :: qo3o3

  INTEGER :: i,k,numitr,ilog,nhalf,it
  REAL :: th2, do3v,do3vp,aplo,aphi

  INCLUDE "seaso3.f90"

  DO i=1,imt
     th2=0.2*xlat(i)
     jjrow(i)=19.001-th2
     tthan(i)=(19-jjrow(i))-th2
  END DO

  !seasonal and spatial interpolation done below.
  DO k=1,nl
     DO i=1,imt
        do3v=xduo3n(jjrow(i),k)+rsin1*xdo3n2(jjrow(i),k)  &
             +rcos1*xdo3n3(jjrow(i),k) +rcos2*xdo3n4(jjrow(i),k)
        do3vp=xduo3n(jjrow(i)+1,k)+rsin1*xdo3n2(jjrow(i)+1,k)  &
             +rcos1*xdo3n3(jjrow(i)+1,k) +rcos2*xdo3n4(jjrow(i)+1,k)
        !now latitudinal interpolation
        !and convert o3 into mass mixing ratio (orig data mpy by 1.e4)
        qo3o3(i,k)=1.e-4*(do3v+tthan(i)*(do3vp-do3v))
     END DO
  END DO

  !vertical interpolation for each gridpoint (linear in ln p)
  numitr=0
  ilog=nl
  DO
     ilog=(ilog+1)/2
     IF(ilog == 1) EXIT
     numitr=numitr+1
  END DO
  DO  k=1,lm
     nhalf=(nl+1)/2
     DO i=1,imt
        jjrow(i)=nhalf
     END DO
     DO  it=1,numitr
        nhalf=(nhalf+1)/2
        DO i=1,imt
           IF(pozn(i,k) < prgfdl(jjrow(i)-1))THEN
              jjrow(i)=jjrow(i)-nhalf
           ELSE IF(pozn(i,k) >= prgfdl(jjrow(i)))THEN
              jjrow(i)=jjrow(i)+nhalf
           END IF
           jjrow(i)=MIN(jjrow(i),nl)
           jjrow(i)=MAX(jjrow(i),2)
        END DO
     END DO
     DO  i=1,imt
        IF(pozn(i,k) < prgfdl(1))THEN
           qo3(i,k)=qo3o3(i,1)
        ELSE IF(pozn(i,k) > prgfdl(nl))THEN
           qo3(i,k)=qo3o3(i,nl)
        ELSE
           aplo=alog(prgfdl(jjrow(i)-1))
           aphi=alog(prgfdl(jjrow(i)))
           qo3(i,k)=qo3o3(i,jjrow(i))+(alog(pozn(i,k))-aphi)/ (aplo-aphi)*  &
                (qo3o3(i,jjrow(i)-1)-qo3o3(i,jjrow(i)))
        END IF
     END DO
  END DO
END SUBROUTINE ozon2d
