      PROGRAM INITOZ5                                                                                                               
C     PROGRAM INITOZ5 (ICTL,       OZLX,       OZONE,       OUTPUT,     )       I2                                                  
C    1          TAPE99=ICTL, TAPE1=OZLX, TAPE2=OZONE, TAPE6=OUTPUT)                                                                 
C     -------------------------------------------------------------             I2                                                  
C                                                                               I2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       I2                                                  
C     JUN 12/96 - F. MAJAESS (REPLACE FSEEK CALL BY FORTRAN CODE)               I2                                                  
C     NOV 06/95 - F. MAJAESS (REVISE DEFAULT PACKING DENSITY VALUE)                                                                 
C     DEC 22/94 - F. MAJAESS(COMPLETE CONVERSION TO VARIABLE DIMENSIONS)                                                            
C     JUL 13/92 - E. CHAN   (DIMENSION SELECTED VARIABLES AS REAL*8)                                                                
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII AND                                                                
C                            REPLACE CALL TO SKIPF WITH A CALL TO THE                                                               
C                            C ROUTINE FSEEK).                                                                                      
C     JUN 07/88 - M.LAZARE. (FIX BUG WHERE IDAY NOT MONTH BOUNDARY)                                                                 
C     NOV 25/86 - M.LAZARE. (BASED ON INITOZC)                                                                                      
C                                                                               I2                                                  
CINITOZ5 - PRODUCES MONTHLY CROSS-SECTIONS OF ZONAL MEAN OZONE                  I1                                                  
C          FOR GCM5                                                     2  1    I1                                                  
C                                                                               I3                                                  
CAUTHORS - J.P.BLANCHET, M.LAZARE                                               I3                                                  
C                                                                               I3                                                  
CPURPOSE - INTERPOLATES INITIAL MONTHLY ZONAL CROSS-SECTIONS OF PPMV            I3                                                  
C          OZONE FROM EQUALLY-SPACED LATITUDES TO THE GAUSSIAN GRID,            I3                                                  
C          THEN INTEGRATES IN THE VERTICAL TO OBTAIN THE VERTICALLY-            I3                                                  
C          INTEGRATED AMOUNT OF OZONE (IN CMS STP) FROM THE TOP OF THE          I3                                                  
C          ATMOSPHERE TO EACH REFERENCE PRESSURE LEVEL IN THE ORIGINAL          I3                                                  
C          OZONE DATA. THE MODEL ROUTINE OZON5 USES THIS RESULT, ALONG          I3                                                  
C          WITH THE SURFACE PRESSURE, TO CALCULATE THE OZONE AMOUNT             I3                                                  
C          BETWEEN FULL MODEL SIGMA SURFACES.                                   I3                                                  
C          IF THE MODEL INITIALIZATION DATE/TIME IS NOT ON A MONTH              I3                                                  
C          BOUNDARY, THE PROGRAM INTERPOLATES BETWEEN ADJACENT MONTH            I3                                                  
C          BOUNDARIES AND STORES THE RESULT ON THE OZONE FILE AS WELL           I3                                                  
C          WITH VALUE OF IDAY IN IBUF(2).                                       I3                                                  
C                                                                               I3                                                  
CINPUT FILES...                                                                 I3                                                  
C                                                                               I3                                                  
C      ICTL  = INITIALIZATION CONTROL DATASET (SEE ICNTRL5).                    I3                                                  
C      OZLX  = MONTHLY LATITUDINAL GLOBAL CROSS-SECTIONS OF OZONE IN            I3                                                  
C              PPMV (ON EQUALLY-SPACED LATITUDES AND A RELATIVELY DENSE         I3                                                  
C              PRESSURE GRID) OBTAINED FROM THE METEOROLOGICAL RESEARCH         I3                                                  
C              REPORT "REFERENCE OZONE MODELS FOR MIDDLE ATMOSPHERE" BY         I3                                                  
C              KITA AND SUMI FROM THE GEOPHYSICAL INSTITUTE, UNIVERSITY         I3                                                  
C              OF TOKYO.                                                        I3                                                  
C                                                                               I3                                                  
COUTPUT FILE...                                                                 I3                                                  
C                                                                               I3                                                  
C      OZONE = MONTHLY GAUSSIAN LATITUDE CROSS-SECTIONS OF VERTICALLY-          I3                                                  
C              INTEGRATED OZONE AMOUNT (CMS STP) FROM THE TOP OF THE            I3                                                  
C              ATMOSPHERE TO EACH INPUT PRESSURE LEVEL,(ON LEVOZ LEVELS,        I3                                                  
C              WHERE LEVOZ IS THE NUMBER OF PRESSURE LEVELS IN THE              I3                                                  
C              ORIGINAL OZONE DATASET IN THE LLPHYS).                           I3                                                  
C              THERE ARE LEVOZ (=37) PRESSURE LEVELS AND ILAT LATITUDES.        I3                                                  
C---------------------------------------------------------------------------                                                        
C                                                                                                                                   
C                                                                                                                                   
C     *    ILAT = SOUTH - NORTH NUMBER OF GRID POINTS.                                                                              
C     *    P    = VECTOR OF PRESSURE LEVELS IN ORIGINAL OZONE DATA (THE VALUES                                                      
C     *           AT P=2000 WERE ADDED TO EMULATE A "SURFACE" FIELD OF OZONE                                                        
C     *           FOR USE IN INTERPOLATING TO MODEL SIGMA LEVELS BELOW 750 MB).                                                     
C     *           (THEY ACTUALLY ARE A COPY OF THE 750 MB DATA).                                                                    
C     *           THE SIZE OF THIS VECTOR IS 37. ITS VALUES ARE GIVEN IN THE                                                        
C     *           DATASET BELOW FOR REFERENCE, ALTHOUGH THEY ARE NOT USED                                                           
C     *           EXPLICITLY IN THE PROGRAM.                                                                                        
C     *    ALAT = VECTOR OF EQUALLY-SPACED LATITUDES (S.P. TO N.P.) IN ORIGINAL                                                     
C     *           OZONE DATA (THE VALUES AT -90. AND 90. WERE ADDED FOR                                                             
C     *           COMPLETENESS BY COPYING THE ADJACENT VALUES AT -80. AND 80.,                                                      
C     *           RESPECTIVELY).                                                                                                    
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
      REAL*8 SL(96),CL(96),WL(96),WOSSL(96),RADL(96)                                                                                
      REAL   DLAT(96)                                                                                                               
      REAL ALAT(19),P(37),CDZPOT(37),OZZX(19,37),                                                                                   
     1     O3(96,37),O3F(96,37)                                                                                                     
      INTEGER NFDM(12),LBL(96)                                                                                                      
C                                                                                                                                   
      COMMON/AOLD/G(65341)                                                                                                          
      COMMON/ANEW/H(96),HM(96),HP(96)                                                                                               
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(65341)                                                                                               
C                                                                                                                                   
      EQUIVALENCE (G(1),OZZX(1,1))                                                                                                  
C                                                                                                                                   
      DATA NFDM/1,32,60,91,121,152,182,213,244,274,305,335/                                                                         
      DATA ALAT/-90.E0,-80.E0,-70.E0,-60.E0,-50.E0,                                                                                 
     &   -40.E0,-30.E0,-20.E0,-10.E0,  0.E0,                                                                                        
     &   10.E0, 20.E0, 30.E0, 40.E0, 50.E0, 60.E0,                                                                                  
     &   70.E0, 80.E0, 90.E0/                                                                                                       
      DATA    P/0.003E0, 0.004E0, 0.005E0, 0.007E0,                                                                                 
     &  0.010E0, 0.015E0, 0.020E0, 0.030E0, 0.040E0,                                                                                
     &  0.050E0, 0.070E0, 0.100E0, 0.150E0, 0.200E0,                                                                                
     &  0.300E0, 0.400E0, 0.500E0, 0.700E0,                                                                                         
     &  1.00E0, 1.50E0, 2.00E0, 3.00E0, 4.00E0,                                                                                     
     &  5.00E0, 7.00E0, 10.0E0, 15.0E0,                                                                                             
     &  20.0E0, 31.0E0, 44.0E0, 88.0E0, 176.E0,                                                                                     
     &  260.E0, 378.E0, 537.E0, 750.E0, 2000.E0/                                                                                    
      DATA C DZ P O T /6.79E-7, 9.52E-7, 1.56E-6, 2.32E-6, 2.95E-6,                                                                 
     1        3.81E-6, 1.02E-5, 7.14E-6, 8.78E-6, 1.60E-5, 2.35E-5,                                                                 
     2        2.88E-5, 3.88E-5, 1.01E-4, 7.24E-5, 8.88E-5, 1.57E-4,                                                                 
     3        2.37E-4, 2.80E-4, 3.86E-4, 1.04E-3, 7.21E-4, 9.03E-4,                                                                 
     4        1.54E-3, 2.35E-3, 2.20E-3, 2.97E-3, 9.11E-3, 1.35E-2,                                                                 
     5        3.66E-2, 6.18E-2, 7.07E-2, 1.05E-1, 1.56E-1, 1.52E-1,                                                                 
     6        1.76E-1, 0.00E0   /                                                                                                   
      DATA MAXX/65341/, MAXLAT/96/                                                                                                  
C-----------------------------------------------------------------------                                                            
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,99,1,2,6)                                                                                                     
C                                                                                                                                   
C     * GET PARAMETERS FROM CONTROL FILE.                                                                                           
C                                                                                                                                   
      REWIND 99                                                                                                                     
      READ(99,END=902) LABL,ILEV,(XXX,L=1,ILEV)                                                                                     
      READ(99,END=903) LABL,IXX,IXX,IXX,ILAT,IXX,IDAY,GMT                                                                           
      WRITE(6,6010)  ILEV,ILAT,IDAY,GMT                                                                                             
C                                                                                                                                   
C     * GET LAT-HGT CROSS-SECTION FOR EACH MONTH (ARRAY NLAT X LEVOZ),                                                              
C     * WHERE NLAT IS THE NUMBER OF EQUALLY-SPACED LATITUDES AND LEVOZ                                                              
C     * IS THE NUMBER OF PRESSURE LEVELS.                                                                                           
C                                                                                                                                   
      DO 290 MONTH=1,12                                                                                                             
      NDAY=NFDM(MONTH)                                                                                                              
      CALL GETFLD2(-1,G,NC4TO8("ZONL"),NDAY,NC4TO8("OZLX"),0,                                                                       
     +                                          IBUF,MAXX,OK)                                                                       
      IF(.NOT.OK) CALL                             XIT('INITOZ5',-1)                                                                
      WRITE(6,6020) IBUF                                                                                                            
      NLAT  = IBUF(5)                                                                                                               
      LEVOZ = IBUF(6)                                                                                                               
C                                                                                                                                   
      IF(ILAT.GT.MAXLAT.OR.NLAT.GT.MAXLAT) CALL    XIT('INITOZ5',-2)                                                                
C                                                                                                                                   
C     * SET DLAT TO GAUSSIAN LATITUDES (DEG).                                                                                       
C                                                                                                                                   
      ILATH=ILAT/2                                                                                                                  
      CALL GAUSSG(ILATH,SL,WL,CL,RADL,WOSSL)                                                                                        
      CALL  TRIGL(ILATH,SL,WL,CL,RADL,WOSSL)                                                                                        
      DO 110 J=1,ILAT                                                                                                               
  110 DLAT(J)=RADL(J)*180.E0/3.14159E0                                                                                              
C                                                                                                                                   
C     * PERFORM INTERPOLATION ON PRESSURE SURFACES FROM EQUALLY-SPACED                                                              
C     * LATITUDES TO GAUSSIAN LATITUDES.                                                                                            
C     * THE ARRAY LBL HOLDS THE POSITION IN THE ARRAY ALAT WHOSE CORRESPONDING                                                      
C     * LATITUDE IS CLOSEST TO, BUT SOUTH OF, (OR POSSIBLY EQUAL TO) THE                                                            
C     * GAUSSIAN LATITUDE DLAT(J); HENCE ONE VALUE FOR EACH GAUSSIAN LATITUDE.                                                      
C     * THE RESULTING ARRAY IS O3(ILAT,LEVOZ).                                                                                      
C                                                                                                                                   
      DO 145 J=1,ILAT                                                                                                               
        IB=0                                                                                                                        
        DO 140 I=1,NLAT                                                                                                             
          IF(DLAT(J).LT.ALAT(I).AND.IB.EQ.0) IB=I                                                                                   
  140   CONTINUE                                                                                                                    
        IF(IB.EQ.0) THEN                                                                                                            
          WRITE(6,6030) DLAT(J)                                                                                                     
          CALL                                     XIT('INITOZ5',-3)                                                                
        ENDIF                                                                                                                       
        LBL(J)=IB-1                                                                                                                 
  145 CONTINUE                                                                                                                      
C                                                                                                                                   
      DO 160 L=1,LEVOZ                                                                                                              
      DO 150 J=1,ILAT                                                                                                               
        K=LBL(J)                                                                                                                    
        CALL LININT(ALAT(K),OZZX(K,L),ALAT(K+1),OZZX(K+1,L),DLAT(J),                                                                
     1              O3(J,L))                                                                                                        
  150 CONTINUE                                                                                                                      
  160 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * EVALUATE THE AMOUNT OF OZONE (IN CMS STP) FROM THE TOP (P(1)=0.003 MBS)                                                     
C     * TO EACH REFERENCE PRESSURE LEVEL. THE RESULTING ARRAY IS                                                                    
C     * O3F(ILAT,LEVOZ).                                                                                                            
C     * USE A FOUR-POINT QUADRATURE SCHEME WITH CONSTANTS C1,C2,C3 AND C4.                                                          
C                                                                                                                                   
      C1=-1.E0/24.E0                                                                                                                
      C2=13.E0/24.E0                                                                                                                
      C3=C2                                                                                                                         
      C4=C1                                                                                                                         
      LEVOZM1=LEVOZ-1                                                                                                               
      DO 175 J=1,ILAT                                                                                                               
        O3F(J,1)=0.E0                                                                                                               
        O3F(J,2)=(C2*O3(J,1)+C3*O3(J,2)+C4*O3(J,3))*C DZ POT(1)                                                                     
        DO 170 L=3,LEVOZM1                                                                                                          
          O3F(J,L)=O3F(J,L-1)+(C1*O3(J,L-2)+C2*O3(J,L-1)+C3*O3(J,L)+                                                                
     1                         C4*O3(J,L+1))*C DZ POT(L-1)                                                                          
  170   CONTINUE                                                                                                                    
        O3F(J,LEVOZ)=O3F(J,LEVOZM1)+(C1*O3(J,LEVOZ-2)+C2*O3(J,LEVOZ-1)+                                                             
     1                               C3*O3(J,LEVOZ))*C DZ POT(LEVOZM1)                                                              
  175 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * OUTPUT CROSS SECTION IN STANDARD FORMAT.                                                                                    
C                                                                                                                                   
      CALL SETLAB(JBUF,NC4TO8("ZONL"),NDAY,NC4TO8("  OZ"),0,ILAT,1,0,1)                                                             
      DO 200 L=1,LEVOZ                                                                                                              
        DO 180 J=1,ILAT                                                                                                             
          H(J)=O3F(J,L)                                                                                                             
  180   CONTINUE                                                                                                                    
        JBUF(4)=L                                                                                                                   
        CALL PUTFLD2(2,H,JBUF,MAXX)                                                                                                 
        WRITE(6,6020) JBUF                                                                                                          
  200 CONTINUE                                                                                                                      
  290 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * ADD EXTRA DATA FOR STARTING DATE (IDAY) OF A NEW MODEL RUN                                                                  
C     * WHICH IS NOT A MONTH BOUNDARY. TO DO THIS, INTERPOLATE BETWEEN                                                              
C     * THE TWO ADJACENT MONTH-BOUNDARY VALUES.                                                                                     
C                                                                                                                                   
      PARTDAY=GMT/24.E0                                                                                                             
      RDAY=FLOAT(IDAY)+PARTDAY                                                                                                      
      LL=0                                                                                                                          
      MM=0                                                                                                                          
      DO 300 N=1,12                                                                                                                 
        IF(IDAY.GE.NFDM(N)) LL=N                                                                                                    
        IF(IDAY.EQ.NFDM(N)) MM=N                                                                                                    
  300 CONTINUE                                                                                                                      
      IF(LL.EQ.0) CALL                             XIT('INITOZ5',-4)                                                                
                                                                                                                                    
      IF(MM.EQ.0) THEN                                                                                                              
                                                                                                                                    
        NDAYM=NFDM(LL)                                                                                                              
        RFDM=FLOAT(NDAYM)                                                                                                           
        IF(LL.LT.12) THEN                                                                                                           
          RFDP=FLOAT(NFDM(LL+1))                                                                                                    
          NDAYP=NFDM(LL+1)                                                                                                          
        ELSE                                                                                                                        
          RFDP=FLOAT(NFDM(1)+365)                                                                                                   
          NDAYP=NFDM(1)                                                                                                             
        ENDIF                                                                                                                       
        WRITE(6,6040) RDAY,NDAYM,NDAYP                                                                                              
C                                                                                                                                   
        DO 500 L=1,LEVOZ                                                                                                            
                                                                                                                                    
          CALL GETFLD2(-2,HM,NC4TO8("ZONL"),NDAYM,NC4TO8("  OZ"),                                                                   
     +                                            L,JBUF,MAXX,OK)                                                                   
          IF(.NOT.OK) CALL                         XIT('INITOZ5',-5)                                                                
          IF(L.EQ.1) THEN                                                                                                           
            WRITE(6,6020) JBUF                                                                                                      
            NPACK=JBUF(8)                                                                                                           
          ENDIF                                                                                                                     
          CALL GETFLD2(-2,HP,NC4TO8("ZONL"),NDAYP,NC4TO8("  OZ"),                                                                   
     +                                            L,JBUF,MAXX,OK)                                                                   
          IF(.NOT.OK) CALL                         XIT('INITOZ5',-6)                                                                
          IF(L.EQ.1) THEN                                                                                                           
            WRITE(6,6020) JBUF                                                                                                      
            NPACK=MIN(NPACK,JBUF(8))                                                                                                
          ENDIF                                                                                                                     
          NWDS=JBUF(5)*JBUF(6)                                                                                                      
          DO 400 I=1,NWDS                                                                                                           
            CALL LININT(RFDM,HM(I),RFDP,HP(I),RDAY,H(I))                                                                            
  400     CONTINUE                                                                                                                  
          JBUF(2)=IDAY                                                                                                              
C                                                                                                                                   
C         * POSITION FILE POINTER TO THE END OF THE FILE AND APPEND DATA.                                                           
C                                                                                                                                   
C         ISTAT=FSEEK(2,0,2)                                                                                                        
  450     READ(2,END=460)                                                                                                           
          GOTO 450                                                                                                                  
  460     BACKSPACE 2                                                                                                               
          JBUF(8)=NPACK                                                                                                             
          CALL PUTFLD2(2,H,JBUF,MAXX)                                                                                               
          WRITE(6,6020) JBUF                                                                                                        
                                                                                                                                    
  500   CONTINUE                                                                                                                    
                                                                                                                                    
      ENDIF                                                                                                                         
C                                                                                                                                   
      CALL                                         XIT('INITOZ5',0)                                                                 
C                                                                                                                                   
C     * E.O.F. ON FILE ICTL.                                                                                                        
C                                                                                                                                   
  902 CALL                                         XIT('INITOZ5',-7)                                                                
  903 CALL                                         XIT('INITOZ5',-8)                                                                
C--------------------------------------------------------------------                                                               
 6010 FORMAT('0 ILEV,ILAT,IDAY,GMT =',3I5,2X,F5.2)                                                                                  
 6020 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
 6030 FORMAT(////,5X,'**** INTERPOLATION ATTEMPTED OUTSIDE THE RANGE OF                                                             
     1THE LATITUDE GRID (-90 TO 90 DEG OF LAT) FOR LAT= ',F6.2,////)                                                                
 6040 FORMAT('0INTERPOLATING FOR RDAY= ',F6.2,' BETWEEN ',I3,' AND ',                                                               
     1           I3,')')                                                                                                            
      END                                                                                                                           
