C**********************************************************                  O3INTERP.2     
C             SUBROUTINE O3INTERP                                            O3INTERP.3     
C**********************************************************                  O3INTERP.4     
      subroutine O3INTERP(o3clim,o3mod,ps)                                   O3INTERP.5     
c                                                                            O3INTERP.6     
c This subroutine takes in a profile from climatology (o3clim) and           O3INTERP.7     
c interpolates to the model vertical grid (nl.gt.15), computing the          O3INTERP.8     
c model profile by constraining its column to agree with climatology.        O3INTERP.9     
c The model profile is passed back in o3mod.                                 O3INTERP.10    
c N.B. The o3clim passed to the routine must be on the pressure levels       O3INTERP.11    
c defined in the routine (ocliml).                                           O3INTERP.12    
c                                                                            O3INTERP.13    
C                                                                            PARAM1.2     
C     Determines model resolution                                            PARAM1.3     
C                                                                            PARAM1.4     
      PARAMETER(NN=21,MM=21,NHEM=2,NL=5,MOCT=1,MG=64,JG=16,NWJ2=121          NETCDF.35    
     P         ,NCRAY=64,JGL=JG,NTRAC=1,NLEVRF=1)                            NETCDF.36    
                                                                             PARAM1.7     
C                                                                            PARAM1.8     
C                                                                            PARAM2.2     
C     Sets basic constants, especially those needed for array dimensions     PARAM2.3     
C                                                                            PARAM2.4     
      PARAMETER(MH=2,PI=3.14159265359,PI2=2.0*PI                             PARAM2.5     
     +,NNP=NN+1,MGPP=MG+2,JGP=JG+1,JGG=JG*NHEM,JGGP=JGG+1,MJP=NWJ2+NWJ2      PARAM2.6     
     +,NLM=NL-1,NLP=NL+1,NLPP=NL+2,NLA=NL+3,NLB=NL+4,NL2=NL*NL               PARAM2.7     
     +,IDA=(MG+MG+MG)/2+1,IDB=NWJ2*NL,IDC=IDB+IDB,IDD=MGPP*NL                PARAM2.8     
     +,IDE=NL2*NN,IDF=NCRAY*(MG+1),IDG=JG*NL,IDH=JG*MG                       PARAM2.9     
     +,IDI=NNP/2,IDJ=IDI*IDI,IDK=NL*IDI,IDL=MGPP/2,IDM=NNP/2,IDN=IDM*NL      PARAM2.10    
     +,NWW=1+(MM-1)/MOCT)                                                    PARAM2.11    
      PARAMETER(IGA=NWJ2*NHEM,IGB=IDB*NHEM,IGC=MGPP*NHEM,IGD=IDD*NHEM        PARAM2.12    
     +,IGG=IDG*NHEM,IGL=IDL*NHEM,IGM=IDM*NHEM,IGN=IDN*NHEM                   PARAM2.13    
     +,IGO=IGA+IGA,IGP=IGB+IGB,NFTWG=(5+NTRAC)*NL+3                          PARAM2.14    
     +,NFTGW=(6+3*NTRAC)*NL+2,NFTGD=(3+NTRAC)*NL,NLTR=NL*NTRAC)              PARAM2.15    
C     Number of 2D (surface) output fields. This value is                    PARAM2.16    
C     Doubled due to averaged and instantaneous fields.                      PARAM2.17    
      PARAMETER (N2DFLD=21,NGRPAD=N2DFLD*2*IGC)                              PARAM2.18    
C                                                                            PARAM2.19    
C                                                                            BLANK.2     
C     Basic planetary parameters for run plus information about              BLANK.3     
C     vertical grid structure                                                BLANK.4     
C                                                                            BLANK.5     
C     Note that RD and GASCON are identical and CPD is set from RD,AKAP.     BLANK.6     
      COMMON        SQ(NNP),RSQ(NNP),SIGMAH(NLM),SIGMA(NL)                   BLANK.7     
     +              ,T01S2(NLM),T0(NL),ALPHA(NL),DSIGMA(NL),RDSIG(NL)        BLANK.8     
     +              ,TKP(NL),C(NL2),SQH(NNP)                                 BLANK.9     
     +              ,MF,MFP,JZF,NF,NFP                                       BLANK.10    
     +              ,AKAP,GA,GASCON,RADEA,WW,P0,PFAC,EZ,AIOCT                BLANK.11    
     +              ,RD,RV,CPD,CLATNT                                        BLANK.12    
     +              ,LRSTRT,LSHORT,LTVEC,LSTRETCH                            BLANK.13    
     +              ,LFLUX                                                   BLANK.14    
     +              ,LBALAN,LRESTIJ                                          BLANK.15    
     +              ,LCLIM, LPERPET, L22L,LOROG ,LCSFCT                      BLANK.16    
     +              ,LNOISE                                                  BLANK.17    
      COMPLEX EZ,AIOCT                                                       BLANK.18    
      LOGICAL LRSTRT,LSHORT,LTVEC,LSTRETCH,LBALAN,LRESTIJ                    BLANK.19    
     +       ,LFLUX,LNOISE                                                   BLANK.20    
     +       ,LCLIM, LPERPET, L22L,LOROG,LCSFCT                              BLANK.21    
C                                                                            BLANK.22    
C                                                                            SPECTR.2     
C     Array ordering in SPECTR must correspond to that in GRIDP.             SPECTR.3     
C                                                                            SPECTR.4     
      COMMON/SPECTR/Z(IGB),D(IGB),T(IGB),TRA(IGB,NTRAC),SP(IGA),GS(IGA)      SPECTR.5     
     :              ,SPA(IGA),VP(IGA),DTE(IGB),TT(IGB)                       SPECTR.6     
     :              ,TRAT(IGB,NTRAC),DT(IGB),ZT(IGB)                         SPECTR.7     
     :              ,ZMI(IGB),DMI(IGB),TMI(IGB)                              SPECTR.8     
     :              ,TRAMI(IGB,NTRAC),SPMI(IGA)                              SPECTR.9     
      COMPLEX Z,D,T,TRA,SP,GS,SPA,VP,DTE,TT,TRAT,DT,ZT                       SPECTR.10    
     :       ,ZMI,DMI,TMI,TRAMI,SPMI                                         SPECTR.11    
C                                                                            SPECTR.12    
c                                                                            O3INTERP.18    
      real o3clim(15),ocliml(15),halflev(14),deltap(15),ohdcol(15)           O3INTERP.19    
      real o3mod(nl),dpmod(nl),modcol(nl)                                    O3INTERP.20    
      real ps, pressu                                                        O3INTERP.21    
      logical ifirst                                                         O3INTERP.22    
      save                                                                   O3INTERP.23    
      data ocliml/1.0,3.0,10.0,30.0,50.0,70.0,100.0,150.0,200.0,             O3INTERP.24    
     &              250.0,300.0,400.0,500.0,700.0,850.0/                     O3INTERP.25    
      data ifirst/.true./                                                    O3INTERP.26    
c                                                                            O3INTERP.27    
      rhoo3=2.14 ! density of O3 at STP                                      O3INTERP.28    
      GR=GA*rhoo3 !g*density of O3 at STP                                    O3INTERP.29    
      if (ifirst) then                                                       O3INTERP.30    
        ifirst=.false.                                                       O3INTERP.31    
        do l=1,14                                                            O3INTERP.32    
          halflev(l)=(ocliml(l+1)+ocliml(l))/2.0    ! in mb                  O3INTERP.33    
        enddo                                                                O3INTERP.34    
        deltap(1)=halflev(1)                        ! Top                    O3INTERP.35    
        do l=2,14                                                            O3INTERP.36    
          deltap(l)=halflev(l)-halflev(l-1)         ! In between             O3INTERP.37    
        enddo                                                                O3INTERP.38    
        deltap(15)=1.0e3-halflev(14)                ! Bottom                 O3INTERP.39    
        do l=1,15                                                            O3INTERP.40    
          deltap(l)=deltap(l)*100.0          ! Convert all to Pa             O3INTERP.41    
        enddo                                                                O3INTERP.42    
      endif                                                                  O3INTERP.43    
                                                                             O3INTERP.44    
                               ! Compute overhead column of climatology      O3INTERP.45    
      col=0.0                  ! for this particular profile.                O3INTERP.46    
      do l=1,15                                                              O3INTERP.47    
        col=col+o3clim(l)*deltap(l)                                          O3INTERP.48    
        ohdcol(l)=col*1.0e5/(GR)    ! In DU                                  O3INTERP.49    
      enddo                                                                  O3INTERP.50    
                                                                             O3INTERP.51    
      do l=1,nl                    ! model delta pressure                    O3INTERP.52    
        dpmod(l)=dsigma(l)*ps      ! in Pa                                   O3INTERP.53    
      enddo                                                                  O3INTERP.54    
                                                                             O3INTERP.55    
                                                                             O3INTERP.56    
      sumrdp=0.0     ! MMR*dp summed down column                             O3INTERP.57    
                                                                             O3INTERP.58    
      do 269 l=1,nl                                                          O3INTERP.59    
        IF (l.LT.NL)THEN                                                     O3INTERP.60    
           pressu=sigmah(l)*ps/100.0   ! in mb                               O3INTERP.61    
        ELSE                                                                 O3INTERP.62    
           pressu=ps/100.0                                                   O3INTERP.63    
        ENDIF                                                                O3INTERP.64    
        do 270 lo=1,14                                                       O3INTERP.65    
          if (pressu.lt.halflev(1)) then                                     O3INTERP.66    
                                                                             O3INTERP.67    
            modcol(l)=(pressu/halflev(1))*ohdcol(1)  ! Model column (DU)     O3INTERP.68    
                                               ! forced from climatology     O3INTERP.69    
            o3mod(l)=(modcol(l)*GR/1.0e5-sumrdp)/dpmod(l)                    O3INTERP.70    
            goto 271                                                         O3INTERP.71    
                                                                             O3INTERP.72    
          elseif (pressu.eq.halflev(lo)) then                                O3INTERP.73    
                                                                             O3INTERP.74    
            modcol(l)=ohdcol(lo)                                             O3INTERP.75    
            o3mod(l)=(modcol(l)*GR/1.0e5-sumrdp)/dpmod(l)                    O3INTERP.76    
            goto 271                                                         O3INTERP.77    
                                                                             O3INTERP.78    
          elseif (pressu.gt.halflev(lo).and.pressu.lt.halflev(lo+1))         O3INTERP.79    
     &         then                                                          O3INTERP.80    
                                                                             O3INTERP.81    
            fraction=(pressu-halflev(lo))/(halflev(lo+1)-halflev(lo))        O3INTERP.82    
            modcol(l)=ohdcol(lo)+fraction*(ohdcol(lo+1)-ohdcol(lo))          O3INTERP.83    
            o3mod(l)=(modcol(l)*GR/1.0e5-sumrdp)/dpmod(l)                    O3INTERP.84    
            goto 271                                                         O3INTERP.85    
                                                                             O3INTERP.86    
          elseif (pressu.ge.halflev(14)) then                                O3INTERP.87    
                                                                             O3INTERP.88    
            fraction=(pressu-halflev(14))/(1.0e3-halflev(14))                O3INTERP.89    
            modcol(l)=ohdcol(14)+fraction*(ohdcol(15)-ohdcol(14))            O3INTERP.90    
            o3mod(l)=(modcol(l)*GR/1.0e5-sumrdp)/dpmod(l)                    O3INTERP.91    
            goto 271                                                         O3INTERP.92    
                                                                             O3INTERP.93    
          endif                                                              O3INTERP.94    
                                                                             O3INTERP.95    
 270  continue                                                               O3INTERP.96    
 271  sumrdp=sumrdp+o3mod(l)*dpmod(l)                                        O3INTERP.97    
 269  continue                                                               O3INTERP.98    
      RETURN                                                                 O3INTERP.99    
      END                                                                    O3INTERP.100   
