SUBROUTINE READSABERATM(LPR,ATMOS,RE,XLAT,MR,NR,NGG,NAM,IGS,ZR,PR,
     &TR,TGR,Q,QGR,ATMFILE,ATMPATH,IHYDRO)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C VERSION 1.0 (c.j. mertens, november 1998)
C
C READSABERATM READS IN MODEL ATMOSPHERIC DATA FOR SABER.
C
C PARAMETER LIST:
C
C   LPR---- UNIT NUMBER OF FILE CONTAINING INPUT MODEL ATMOSPHERIC DATA.
C   RE----- EARTH RADIUS.
C   XLAT--- GEOCENTRIC LATITUDE OF TANGENT POINT.
C   MR----- MAXIMUM NUMBER OF TANGENT ALTITUDES IN MODEL ATMOSPHERE.
C   NR----- ACTUAL NUMBER OF TANGENT ALTITUDES IN MODEL ATMOSPHERE.
C   NGG---- NUMBER AS GASES IN MODEL ATMOSPHERE.
C   NAM---- GASNAMES IN MODEL ATMOSPHERE.
C   IGS---- AFGL GASID's FOR GASES IN MODEL ATMOSPHERE.
C   ZR----- TANGENT ALTITUDES IN MODEL ATMOSPHERE.
C   PR----- PRESSURE AT TANGENT ALTITUDES IN MODEL ATMOSPHERE.
C   TR----- TEMPERATURE AT TANGENT ALTITUDES IN MODEL ATMOSPHERE.
C   TGR---- TEMPERATURE GRADIENTS AT TANGENT ALTITUDES IN MODEL.
C   Q------ VOLUME MIXING RATIOS FOR GASES IN MODEL ATMOSPHERE.
C   QGR---- MIXING RATIO GRADIENTS AT TANGENT ALTITUDES.
C   ATMPATH DIRECTORY PATH FOR MODEL ATMOSPHERIC DATA.
C
C SOFTWARE WRITTEN BY Christopher J. Mertens
C                     NASA Langley Research Center
C                     Atmospheric Sciences Competency
C                     21 Langley Blvd., Mail Stop 401B
C                     Hampton, VA 23681-2199
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PARAMETER(NFILE=9,MTMP=300)
      CHARACTER*(*) ATMOS,ATMFILE,ATMPATH,FILEN(NFILE)*10,BLANK*1
      CHARACTER*8 NAM(*),NAME(NFILE)
      LOGICAL IHYDRO
      REAL*8 REARTH
      DIMENSION ZR(MR),PR(MR),TR(MR),TGR(2,MR),IGS(*)
      DIMENSION Q(MR,*),QGR(2,MR,*),ZTMP(MTMP),QTMP(MTMP),PTMP(MTMP)
      DIMENSION ZTMP2(MTMP),QTMP2(MTMP),IGAS(NFILE)
      DATA BLANK/' '/
      DATA FACTOR/1.0E-6/
      DATA FILEN/ 'H2O.ALT', 'CO2.ALT', 'O3.ALT', 'O2.ALT', 'N2.ALT',
     &'O.ALT', 'O1D.ALT', 'TEMP.ALT', 'PRESS.ALT'/ 
      DATA NAME/ 'H2O', 'CO2', 'O3', 'O2', 'N2', 'O(3P)','O(1D)','TEMP',
     &'PRESS'/
      DATA IGAS/100,200,300,700,2200,3400,6100,0,0/
      NGG=NFILE-2
C
C ..... Find directory path to atmospheric data.
C
      NCHAR=LEN(ATMFILE)
      DO IC=1,NCHAR
        IF(ATMFILE(IC:IC).EQ.BLANK) GO TO 5
        ICC=IC
      ENDDO
  5   ATMOS=ATMFILE(1:ICC)
      ATMPATH='./atmos/'//ATMOS//'/'
      DO IC=1,NCHAR
	ICC1=IC
        IF(ATMPATH(IC:IC).NE.BLANK) GO TO 1
      ENDDO	
  1   DO IC=ICC1,NCHAR
	IF(ATMPATH(IC:IC).EQ.BLANK) GO TO 2
	ICC2=IC
      ENDDO
      ATMOS=ATMPATH(ICC1:ICC2)
C
C ..... Read latitude.
C
  2   ATMFILE=ATMPATH(ICC1:ICC2)//'LATITUDE.DAT'
      OPEN(LPR,FILE=ATMFILE,ERR=980)
      READ(LPR,'(F7.2)') XLAT
      CLOSE(LPR)
C
C ..... Read pressure profile.
C
      ATMFILE=ATMPATH(ICC1:ICC2)//FILEN(NFILE)
      OPEN(LPR,FILE=ATMFILE,ERR=980)
      READ(LPR,*)
      I=0
  10  I=I+1
      READ(LPR,*,END=15) ZTMP(I),PTMP(I)
      GO TO 10
  15  NR=I-1
      IF(NR.EQ.0) THEN
        WRITE(*,*) '!!!! COULD NOT FIND ',ATMFILE
        STOP 'ERROR IN READSABERATM'
      ENDIF
      IF(NR.GT.MR) THEN
        WRITE(*,*) '!!!! DIMENSION ERROR. INCREASE MR TO AT LEAST ',NR
	STOP 'ERROR IN READSABERATM'
      ENDIF
      IF(NR.GT.MTMP) THEN 
    	WRITE(*,900) NR
        STOP 'ERROR IN READSABERATM'
      ENDIF 
      DO I=1,NR
        II=NR-I+1
        ZR(II)=ZTMP(I)
        PR(II)=PTMP(I)
      ENDDO
      CLOSE(LPR)           
C
C ..... Read temperature profile.
C
      ATMFILE=ATMPATH(ICC1:ICC2)//FILEN(NFILE-1)
      OPEN(LPR,FILE=ATMFILE,ERR=980)
      READ(LPR,*)
      I=0
  20  I=I+1
      READ(LPR,*,END=25) ZTMP(I),QTMP(I)
      GO TO 20
  25  NTMP=I-1
      IF(NTMP.EQ.0) THEN
        WRITE(*,*) '!!!! COULD NOT FIND ',ATMFILE
        STOP 'ERROR IN READSABERATM'
      ENDIF
      IF(NTMP.GT.MTMP) THEN
        WRITE(*,900) NTMP
        STOP 'ERROR IN READSABERATM'
      ENDIF
      DO I=1,NTMP
        II=NTMP-I+1
        ZTMP2(II)=ZTMP(I)
        QTMP2(II)=QTMP(I)
      ENDDO
      CLOSE(LPR)      
      CALL KP_INTPROF(NR,NTMP,1,NR,ZTMP2,ZR,QTMP2,TR)
C
C ..... Initialize gradients to zero.
C      
      DO I=1,NR
        DO IS=1,2
          TGR(IS,I)=0.0
          DO IG=1,NGG
            QGR(IS,I,IG)=0.0
          ENDDO
        ENDDO
      ENDDO
C
C ..... Read constituent profiles.
C
      DO IG=1,NGG
        IGS(IG)=IGAS(IG)
        NAM(IG)=NAME(IG)
        ATMFILE=ATMPATH(ICC1:ICC2)//FILEN(IG)
        OPEN(LPR,FILE=ATMFILE)
        READ(LPR,*)
        I=0
  30    I=I+1
        READ(LPR,*,END=35) ZTMP(I),QTMP(I)
        GO TO 30
  35    NTMP=I-1
        IF(NTMP.EQ.0) THEN
          WRITE(*,*)
          WRITE(*,910) ATMFILE,NAM(IG)
          DO IL=1,NR
            Q(IL,IG)=0.0
          ENDDO
          GO TO 40
        ENDIF
        IF(NTMP.GT.MTMP) THEN
          WRITE(*,900) NTMP
          STOP 'ERROR IN READSABERATM'
        ENDIF
        DO I=1,NTMP
          II=NTMP-I+1
          ZTMP2(II)=ZTMP(I)
          QTMP2(II)=QTMP(I)*FACTOR
        ENDDO
        CLOSE(LPR)           
        CALL KP_INTPROF(NR,NTMP,1,NR,ZTMP2,ZR,QTMP2,Q(1,IG))
   40   CONTINUE
      ENDDO
C      WRITE(88,*) NR,XLAT,RE
C      DO I=1,NR
C        WRITE(88,999) I,ZR(I),TR(I),PR(I),Q(I,1),Q(I,2),Q(I,3),Q(I,4),
C     &  Q(I,5)
C      ENDDO
C      RE=REARTH(XLAT)
      IPTIE=NR
      CALL RADGRAV(ZR,XLAT,G,RE)
C      CALL BARO(MR,NR,1,NR,IPTIE,NGG,IGS,RE,XLAT,ZR,PR,TR,Q,.TRUE.)
      IF(IHYDRO) CALL HYDRO(NR,1,NR,IPTIE,RE,XLAT,ZR,PR,TR)
C      DO I=1,NR
C        TR(I)=TR(I)+2.0
C      ENDDO
      RETURN
  900 FORMAT(1X,'!!!! DIMENSION ERROR IN READSABERATM. INCREASE MTMP',
     &       1X,'TO AT LEAST ',I5)
  910 FORMAT(1X,'!!!! WARNING: COULD NOT FIND ',A80,/,' SETTING MIXING',
     &       1X,'RATIOS EQUAL TO ZERO FOR GAS: ',A8)
  980 WRITE(*,990) ATMFILE
  990 FORMAT(1X,'PROBLEM FINDING THE FOLLOWING FILE IN READSABERATM:',/,
     &       1X,A80)
 999  FORMAT(1X,I2,2F9.3,6E12.5)
      END