!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
***S/P OZOREF2
*

      SUBROUTINE OZOREF2(O3F,LREF,DLAT,NP,NMAX,LBL,NLAT,ALAT, 1
     x     PREF,F)
#include "impnone.cdk"
      INTEGER LREF,NP,NMAX,IB,I,J,K,L,NLAT
      REAL SLOPE,XLATI,PREF(LREF)
      REAL DL,B,F1,F2,A1,A2,OZON,FAC,CONV
      REAL O3F(NMAX,LREF)
      INTEGER LBL(NP)
      REAL DLAT(NP)
      REAL F(NLAT,LREF),ALAT(NLAT)
*
*Author
*          L.Garand (1997)
*          rewritten from original CCRN code
*
*Revisions
* 001      B. Bilodeau (Jan 2000) - Exit if latitudes out of bounds
*
*Object
*          to calculate the ozone mixing ratio (kg/kg) at ozone
*          climatological levels for desired array of latitudes
*
*Arguments
*
*          - Output -
* O3F      ozone (kg O3/kg air) at  each reference
*          climatological level for DLAT latitudes
*
*          - Input -
* PREF     reference ozone pressure (N/m2) levels of field F
* LREF     number of climatological ozone level
* DLAT     latitude in radians of model points
* NP       number of points to process
* NMAX     number of maximum points permitted
* LBL      work field
* NLAT     number of latitude climatological bands
* ALAT     climatological ozone latitudes in degrees
* F        climatological ozone field in PPMV
*
**
*
C 
#include "consphy.cdk"
C
      FAC = 180./PI
C
      DO 145 J=1,NP
        IB=0
        XLATI= FLOAT( NINT(DLAT(J)*FAC) )
C
        DO 140 I=1,NLAT
          IF( XLATI.LT.ALAT(I) .AND. IB.EQ.0 ) IB=I
  140   CONTINUE
C
        IF ( XLATI.EQ. 90.0 ) IB=NLAT
C
        IF(IB.LE.1) THEN
          WRITE(6,6030) XLATI
          WRITE(6,*) (ALAT(I),i=1,NLAT)
          CALL QQEXIT(1)
          RETURN
 6030     FORMAT(1X,' O3 INPOL OUT BOUNDS IN LATITUDE:',E12.4)
        ENDIF
        LBL(J)=IB-1
  145 CONTINUE
C
C  interpolate to desired latitudes
c  and transform into kg/kg using CONV:
C  1.E-6 converts PPMV to PPV , 48.0=M(O3), 28.964= M(dry air)
c
      CONV = 1.E-6*48./28.964

c
      DO 160 L=1,LREF
      DO 160 J=1,NP
      K=LBL(J)
      A1=ALAT(K)
      A2=ALAT(K+1)
      F1=F(K,L)
      F2=F(K+1,L)
      DL=DLAT(J)  *FAC
      SLOPE = (F2-F1)/(A2-A1)
      B = F1 - SLOPE*A1
      OZON = SLOPE*DL + B
      O3F(J,L)=OZON*CONV
  160 CONTINUE
C
C
      RETURN
      END