      SUBROUTINE WACCMC2_ATMSTATE_INTERP(NTIMES,NLATITUDES,NLEVELS,
     &Z_WACCM,P_WACCM,ZGEO_WACCM,T_WACCM,
     &H2O_WACCM,CO2_WACCM,O3_WACCM,N2O_WACCM,CH4_WACCM,O_WACCM,
     &O2_WACCM,NO_WACCM,NO2_WACCM,OH_WACCM,N2_WACCM,HO2_WACCM,
     &CFCL3_WACCM,CF2CL2_WACCM,H_WACCM,O1D_WACCM,
     &MR,MGM,IDG,RE,XLAT,ZR,PR,TR,Q,CO2,CO2_FIX)
************************************************************************
* WACCM_ATMSTATE CREATES ATMOSPHERIC STATE VECTORS FROM THE WACCM
* MODEL DATABASE ARRAYS. THERE IS ONE FILE PER MONTH THROUGHOUT THE
* SABER MISSION.
*
* ZGEO_WACCM IS GEOPOTENTIAL ALTITUDE IN METERS WE MUST CONVERT THAT
* TO GEOMETRIC ALTITUDE IN KM. WE TAKE ALL PROFILES 
* TO O KM BY INSISTING THAT THE VALUE AT O KM FOR A PARTICULAR 
* PROFILE IS SET EQUAL TO THE CORRESPONDING VALUE AT 2 KM. THERE 
* ARE 46 LATITUDE DATA POINTS FOR EACH ATMOSPHERIC CONDITION.
* LATITUDE COVERAGE IS FROM -90 to +90 AT 4 DEGREE INCREMENTS.
* 
************************************************************************
************************************************************************
*
*      (c) Copyright 1987-2001 by GATS, Inc.
*          11864 Canon Blvd. Suite 101
*          Newport News, VA 23606
*          Phone: (757) 873-5920
*
*  All Rights Reserved. No part of this software or publication may be
*  reproduced, stored in a retrieval system, or transmitted, in any form
*  or by any means, electronic, mechanical, photocopying, recording, or
*  otherwise without the prior written permission of GATS, Inc.
*
************************************************************************
      IMPLICIT   NONE

      REAL*4     R,DZ,RI,WGHT,WGHT2,GREF
      PARAMETER  (R = 287.0, GREF=9.8)
      INTEGER*4  MAXLEV,CO2NUM
      PARAMETER  (MAXLEV = 1000)
      LOGICAL    IVAR,IQUIT
      INTEGER*4  NTIMES,LEVELS,MR,MGM,IDG(MGM)
      INTEGER*4  J,IFIX,IPTIE
      INTEGER*4  NLATITUDES,NLEVELS,I,II,III,IG,II2
      REAL*4    G,RE,REARTH,XLAT,CO2I,CO2MAX,ZMBEG,ZMEND,FACT,XLATINC,ZG
      REAL*4     Q(MR,MGM),ZR(MR),PR(MR),TR(MR),CO2(MR),CO2_FIX,CO2_FACT
      REAL*4     Z_WACCM(NLEVELS),P_WACCM(NLEVELS)
      REAL*4     ZGEO_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     T_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     H2O_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     CO2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     O3_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     N2O_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     CH4_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     O_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     O2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     NO_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     NO2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     OH_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     N2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     HO2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     CFCL3_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     CF2CL2_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     H_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     O1D_WACCM(NLEVELS,NLATITUDES,NTIMES)
      REAL*4     EVAR(MAXLEV)
      REAL*4     CO2PERB(MAXLEV)
      REAL*4     CO2TEMP(MAXLEV)

      IF( MR .LE. NLEVELS )THEN
         WRITE(*,*) 'MR must be at least ',NLEVELS+1
         STOP 'ERROR in ATMSTATE'
      ENDIF
      IF( MAXLEV .LT. MR )THEN
         WRITE(*,*) 'MAXLEV must be at least ',MR
         STOP 'ERROR in ATMSTATE'
      ENDIF

      IF( MGM .LT. 16 )THEN
         WRITE(*,*) 'MGM must be at least 16'
         STOP 'ERROR in ATMSTATE'
      ENDIF

*
* ... Determine index for tangent-point latitude
*
      I=1
      IF(XLAT.LT.-90.0.OR.XLAT.GT.90.0) THEN
        WRITE(*,*) 'ERROR: XLAT not in range'
      ENDIF
      XLATINC=180.0/(NLATITUDES-1)
      RI=1.0+(XLAT+90.0)/XLATINC
      II=INT(RI)
      IF(II.LT.1) THEN
        II=1
        WGHT=1.0
      ELSEIF(II.GT.NLATITUDES-1) THEN
        II=NLATITUDES-1
        WGHT=0.0
      ELSE
        WGHT=1.0-(RI-REAL(II))
      ENDIF
      IF(WGHT.LT.0.0.OR.WGHT.GT.1.0) THEN
        WRITE(*,*) 'ERROR: XLAT not in range'
      ENDIF        
      II2=II+1
      WGHT2=1.0-WGHT
CTM Set the q values. Not all are available now. 
      DO III = 1, NLEVELS
        PR(III)   = P_WACCM(III)
CTM Set ZR at this point to Geopotential altitude
        ZR(III)  = WGHT*ZGEO_WACCM(III,II,I)+WGHT2*ZGEO_WACCM(III,II2,I)
        TR(III)   = WGHT*T_WACCM(III,II,I)+WGHT2*T_WACCM(III,II2,I)
        Q(III,1)  = WGHT*O2_WACCM(III,II,I)+WGHT2*O2_WACCM(III,II2,I)
        Q(III,2)  = WGHT*CO2_WACCM(III,II,I)+WGHT2*CO2_WACCM(III,II2,I)
        Q(III,3)  = WGHT*O3_WACCM(III,II,I)+WGHT2*O3_WACCM(III,II2,I)
        Q(III,4)  = WGHT*H2O_WACCM(III,II,I)+WGHT2*H2O_WACCM(III,II2,I)
        Q(III,5)  = WGHT*N2O_WACCM(III,II,I)+WGHT2*N2O_WACCM(III,II2,I)
        Q(III,6)  = WGHT*CH4_WACCM(III,II,I)+WGHT2*CH4_WACCM(III,II2,I)
        Q(III,7)  = WGHT*NO_WACCM(III,II,I)+WGHT2*NO_WACCM(III,II2,I)
        Q(III,8)  = WGHT*NO2_WACCM(III,II,I)+WGHT2*NO2_WACCM(III,II2,I)
        Q(III,9)  = WGHT*OH_WACCM(III,II,I)+WGHT2*OH_WACCM(III,II2,I)
        Q(III,10) = WGHT*N2_WACCM(III,II,I)+WGHT2*N2_WACCM(III,II2,I)
        Q(III,11) = WGHT*HO2_WACCM(III,II,I)+WGHT2*HO2_WACCM(III,II2,I)
        Q(III,12) = WGHT*O_WACCM(III,II,I)+WGHT2*O_WACCM(III,II2,I)
        Q(III,13) = WGHT*H_WACCM(III,II,I)+WGHT2*H_WACCM(III,II2,I)
        Q(III,14) = WGHT*CFCL3_WACCM(III,II,I)+
     &              WGHT2*CFCL3_WACCM(III,II2,I)
        Q(III,15 )= WGHT*CF2CL2_WACCM(III,II,I)+
     &              WGHT2*CF2CL2_WACCM(III,II2,I)
        Q(III,16) = WGHT*O1D_WACCM(III,II,I)+WGHT2*O1D_WACCM(III,II2,I)
      ENDDO
      CALL RADGRAV(0.0,XLAT,G,RE)
      CO2MAX=0.0
      CO2NUM=0
      DO III=1,NLEVELS
        CO2I=Q(III,2)
CTM calculate ZR from Geopotential Altitude, 0.001 factor converts to km
        ZG=0.001*ZR(III)*GREF/G
C        write(99,*) III,RE,ZG
        ZR(III)=RE*ZG/(RE-ZG)
C        write(99,*) III,ZR(III)
        IF(PR(III).GT.1.0.AND.PR(III).LT.10.0) THEN
          CO2MAX=CO2MAX+CO2I
          CO2NUM=CO2NUM+1
        ENDIF
      ENDDO
      CO2MAX=CO2MAX/CO2NUM
      WRITE(*,*) 'CO2: TREND, WACCM_REFC2',CO2_FIX,CO2MAX
CTM USE WACCM as is
C FIX CO2 PROFILE
CTM First check whether CO2 VMR exceeds CO2_FIX. If so then scale upper
CTM altitudes using last value .LE. CO2_FIX
C      CO2MAX=0
C      DO III=1,NLEVELS
C        CO2I=Q(III,2)
C        IF(PR(III).GT.0.05.AND.PR(III).LT.0.5) THEN
C          CO2MAX=MAX(CO2MAX,CO2I)
C          IF(CO2I.GT.CO2_FIX) THEN
C            CO2_FACT=CO2_FIX/CO2I
C            IFIX=III
C            GOTO 11
C          ENDIF
C        ENDIF
C      ENDDO
CTM IF the above condition isn't met then scale upper altitudes using
CTM maximum value in array as long as it's close enough
C      CO2_FACT=CO2_FIX/CO2MAX
C      IF(CO2_FACT.LT.1.02) THEN
C        DO III=1,NLEVELS
C          CO2I=Q(III,2)
C          IF(PR(III).GT.0.05.AND.PR(III).LT.0.5.AND.CO2I.EQ.CO2MAX) THEN
C            IFIX=III
C            GOTO 11
C          ENDIF
C        ENDDO
C      ENDIF
C      WRITE(*,*) 'ERROR FIXING CO2 PROFILE'
C      STOP 'ERROR IN ATMSTATE'
C 11   DO III=1,NLEVELS
C         IF(III.LT.IFIX) THEN
C           Q(III,2)=CO2_FACT*Q(III,2)
C         ELSE
C           Q(III,2)=CO2_FIX
C         ENDIF
C         CO2(III)=Q(III,2)
C      ENDDO
CTM Perturb CO2 profile above 70km. Start with 0% change at 70km and then
CTM slowly increase to 15% at 80km and keep it at 15% above 80km.
CTM Do this by first creating CO2PERB=CO2*1.15 and splicing with CO2
CTM using ZMBEG=80km and ZMEND=70km. Write out warning to standard out
CTM so it is obvious this is being done.
CTM Set up parameters for splicing perturbed CO2 to WACCM for test
CTM TURN THIS OFF!!!
C      IVAR = .FALSE.
C      ZMBEG = 80.0
C      ZMEND = 70.0
C      FACT=1.0
C      WRITE(*,*) "!!!!WARNING: Perturbing CO2 for test run :WARNING!!!!"
C      DO III=1,NLEVELS
C         CO2PERB(III)=CO2(III)*1.15
C      ENDDO
C      CALL FUSE(MR,NLEVELS,ZMBEG,ZMEND,ZR,CO2,CO2PERB,FACT,IVAR,
C     &             EVAR,EVAR,CO2TEMP,IQUIT)
C      DO III=1,NLEVELS
C         WRITE(*,*) ZR(III),CO2PERB(III),CO2(III),CO2TEMP(III)
C         CO2(III)=CO2TEMP(III)
C      ENDDO
CTM END CO2 Perturb 
*
      ZR(MR) = 0.0
      DZ = 1000.0*ZR(NLEVELS)
      TR(MR) = TR(NLEVELS)
      PR(MR) = PR(NLEVELS)*EXP(G*DZ/R/TR(NLEVELS))
      CO2(MR)= CO2(NLEVELS)
      DO IG = 1, MGM
         Q(MR,IG) = Q(NLEVELS,IG)
      ENDDO
CTM Do not call Hypsometric, Pressures are not hydrostatic at this point.
CTM HYPSOMETRIC recalculates altitudes for hydrostatic balance
CTM Hysometric integrates from top need to integrate from bottom
C      CALL HYPSOMETRIC(1,MR,MR,RE,XLAT,MR,MGM,IDG,ZR,PR,TR,Q)
C      IPTIE=MR
C      DO III=1,MR
C        IF(PR(III).LT.10.0) IPTIE=III
C      ENDDO
C      CALL HYDROZ(1,MR,IPTIE,RE,XLAT,ZR,PR,TR)
*
*      DO III = 1, MR
*         WRITE(45,101) III, II, IREC, DAY, ZR(III), PR(III), TR(III)
*      ENDDO
*
*101   FORMAT(I3,2X,I2,2X,I3,2X,I3,2X,1PE12.5,2X,1PE12.5,2X,1PE12.5)
*      DO I = 1, MR
*         WRITE(46,102) I,ZR(I),Q(I,1),Q(I,2),Q(I,3),Q(I,4)
*      ENDDO
*102   FORMAT(I2,5(1X,1PE12.5))
*
      RETURN
      END                                                                           
