C =============================================================================

C ROUTINE TO CALL MSIS AND SET UP REFERENCE BACKGROUND ATMOSPHERE ON A
C "FINE" ALTITUDE GRID SUITABLE FOR USE IN lyao_rt-RELATED RADIATIVE
C TRANSPORT CALCULATIONS.

C -----------------------------------------------------------------------------
C NOTES:
C   - version for generating hab_rt "fine grid" thermosphere
C   - IYEAR value not really relevant in MSIS, so 2-digit depiction is kept
C     (i.e., "00" for "2000").
C -----------------------------------------------------------------------------

      SUBROUTINE GET_MSIS(LINE_LABEL)

      IMPLICIT REAL*4 (A-H,O-Z)

      DIMENSION AP(7)
      REAL*8 THERMO_MSIS(7,61)
      REAL*8 TGCM_Z,TGCM_TN,TGCM_H,TGCM_O2,TGCM_N2,TGCM_O1
      CHARACTER(LEN=256) LINE
      
      COMMON/PARA_MSIS/ GLAT,GLONG,IDAY,IYEAR,UT,STL,IAPH,AP,F107,F107A
      COMMON/FINE_GRID/ THERMO_MSIS,BASE,TOP,ITHERM

      DIMENSION TMSIS(2),DMSIS(8),SV(25)
      COMMON/CSW/ SW(25),ISW,SWC(25)
      DATA SV   / 25 * 1.0 /
      DATA ISW  /     1    /

      DATA PLANETR / 6371.0E5 /

      ITHERM = 60

C -----------------------------------------------------------------------------
C "FINE GRID" THERMOSPHERE MODEL & RELATED QUANTITIES
C -----------------------------------------------------------------------------

C ALTITUDE GRID SETTINGS [see get_msis_lyao.f]

      IF (LINE_LABEL.EQ.2) THEN
        BASE = 102.0                              ! Lyman beta
      ELSE IF (LINE_LABEL.EQ.3) THEN
        BASE = 160.0                              ! Lyman gamma
      END IF

      QFCTR  = 0.075128
      QPIVOT = BASE

      TOP    = BASE + 392.0                       ! exobase

C CONSTRUCTING THERMOSPHERE REFERENCE TABLE

C _____ PREPARING TO CALL MSIS

      MASS48 = 48
      IYEAR  = MOD(IYEAR,100)                     ! 2-digit "year" value
      IYD    = IDAY + 1000*IYEAR
      SEC    = 3600.0 * UT
      IF (IAPH.EQ.1) THEN
        CALL TSELEC(SV)
        SW(9) = -1.0
      END IF

C _____ THE BIG CALL

      ALT = BASE
      
      
C     Open the file and read the header, also set up the format statement for the data
 12   FORMAT(f7.3,f9.3,4e14.7)
      OPEN(UNIT=15,FILE='tgcm_JunSMAX.dat')
      DO I=1,5
        READ(15,*) line
      ENDDO
      
      DO IBIN = 1,ITHERM+1

        IF (IBIN.EQ.1) THEN
          QUAD = 0.0
        ELSE
          QUAD = QFCTR * ABS(ALT-QPIVOT)
          IF (QUAD.LT.1.0) QUAD = 1.0
        END IF

        ALT = ALT + QUAD
        IF (IBIN.EQ.ITHERM+1) ALT = TOP

C        CALL GTD6(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS48,
C     &            DMSIS,TMSIS)
        THERMO_MSIS(1,IBIN) = DBLE(ALT)                   ! ALTITUDE
        THERMO_MSIS(2,IBIN) = DBLE(PLANETR + ALT*1.0E5)   ! RADIUS
C        THERMO_MSIS(3,IBIN) = DBLE(TMSIS(2))              ! TEMPERATURE
C        THERMO_MSIS(4,IBIN) = DBLE(DMSIS(7))              ! ATOMIC HYDROGEN
C        THERMO_MSIS(5,IBIN) = DBLE(DMSIS(4))              ! MOLECULAR OXYGEN
C        THERMO_MSIS(6,IBIN) = DBLE(DMSIS(3))              ! MOLECULAR NITROGEN
C        THERMO_MSIS(7,IBIN) = DBLE(DMSIS(2))              ! ATOMIC OXYGEN
        
        READ (15,12) TGCM_Z,TGCM_TN,TGCM_H,TGCM_O2,TGCM_N2,TGCM_O1
        THERMO_MSIS(3,IBIN) = TGCM_TN
        THERMO_MSIS(4,IBIN) = TGCM_H
        THERMO_MSIS(5,IBIN) = TGCM_O2
        THERMO_MSIS(6,IBIN) = TGCM_N2
        THERMO_MSIS(7,IBIN) = TGCM_O1

      END DO

      CLOSE(15)

C -----------------------------------------------------------------------------

      RETURN
      END

C =============================================================================