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

C SET UP GRIDDED ATMOSPHERE DENSITIES & ATMOSPHERIC TEMPERATURE FOR
C MULTIPLE SCATTERING CALCULATIONS
C [29Jan02 JBishop]

C -----------------------------------------------------------------------------
C ZONE DENSITIES & TEMPERATURES:  SPHERICALLY SYMMETRIC
C
C   HDNS  . . . ATOMIC HYDROGEN DENSITIES
C   O2DNS . . . O2 DENSITIES
C   TATM  . . . ATMOSPHERIC TEMPERATURE
C -----------------------------------------------------------------------------

      SUBROUTINE RT_GRID(ITHERM,THERMO,IEXO,IKNT,BRAD,HDNS,O2DNS,TATM)

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

      DIMENSION THERMO(5,61)
      DIMENSION BRAD(IKNT+1)
      DIMENSION HDNS(IKNT), O2DNS(IKNT), TATM(IKNT)

      COMMON /EXOS/  GM,CONL,PLANETR,RBASE,RC,RUPR,RP,TEXO,DEXO,VELT,
     &               RADPF,CENTER,ABSCSX,FTSAT,FDSAT,IGEO
      COMMON /NMBR/  PI,RTPI,PID2,OFFSET
      COMMON /GAUS08/WW08(08),XX08(08)
      COMMON /GAUS16/WW16(16),XX16(16)
      COMMON /GAUS32/WW32(32),XX32(32)

C -----------------------------------------------------------------------------
C THERMOSPHERE:  PLANE PARALLEL, QUASI-ISOTHERMAL TREATMENT
C -----------------------------------------------------------------------------

      ZSOFAR = BRAD(1)
      HCOLM  = 0.0D0
      O2COLM = 0.0D0
      TWHYD  = 0.0D0

      II = 1

      DO IBIN = 1,ITHERM

        INXT  = IBIN + 1

        DZ    = THERMO(2,INXT) - THERMO(2,IBIN)

        H1N   = THERMO(4,IBIN)
        O2N   = THERMO(5,IBIN)

        H1H   = DZ / LOG(H1N/THERMO(4,INXT))
        O2H   = DZ / LOG(O2N/THERMO(5,INXT))

        H1C   = H1N*H1H * (1.0D0 - EXP(-DZ/H1H))
        O2C   = O2N*O2H * (1.0D0 - EXP(-DZ/O2H))

        TPL   = THERMO(3,IBIN)
        TGRAD = (THERMO(3,INXT) - TPL) / DZ

        ZTEST = ZSOFAR + DZ

        IF ( (ZTEST.GT.BRAD(II+1)) .OR. (IBIN.EQ.ITHERM) ) THEN

 401      CONTINUE

          DR  = BRAD(II+1) - ZSOFAR

          H1R = H1N*H1H * (1.0D0 - EXP(-DR/H1H))
          O2R = O2N*O2H * (1.0D0 - EXP(-DR/O2H))
          TPR = TPL + DR * TGRAD

          HCOLM  = HCOLM  + H1R
          O2COLM = O2COLM + O2R
          TWHYD  = TWHYD  + H1R * (TPL + TPR) / 2.0D0

          HDNS(II)  = HCOLM  / (BRAD(II+1)-BRAD(II))
          O2DNS(II) = O2COLM / (BRAD(II+1)-BRAD(II))
          TATM(II)  = TWHYD  / HCOLM

          IF (II.EQ.IEXO) GO TO 402

          ZSOFAR = BRAD(II+1)

          II = II + 1

          IF (ZTEST.GT.BRAD(II+1)) THEN

            HCOLM  = 0.0D0
            O2COLM = 0.0D0
            TWHYD  = 0.0D0

            H1N    = H1N * EXP(-DR/H1H)
            O2N    = O2N * EXP(-DR/O2H)
            TPL    = TPL + DR * TGRAD

            H1C    = H1C - H1R
            O2C    = O2C - O2R

            GO TO 401

          END IF

          HCOLM  = H1C - H1R
          O2COLM = O2C - O2R
          TWHYD  = HCOLM * (TPR + THERMO(3,INXT)) / 2.0D0

          ZSOFAR = ZTEST

        ELSE

          HCOLM  = HCOLM  + H1C
          O2COLM = O2COLM + O2C
          TWHYD  = TWHYD  + H1C * (TPL + THERMO(3,INXT)) / 2.0D0

          ZSOFAR = ZTEST

        END IF

      END DO

 402  CONTINUE

C -----------------------------------------------------------------------------
C LOCAL EXOSPHERE:  ANALYTIC EXTENSION OF THERMOSPHERE PROFILE
C -----------------------------------------------------------------------------

      DO II = IEXO+1,IKNT

        AUPR = CONL / TEXO / BRAD(II)
        ALWR = CONL / TEXO / BRAD(II+1)

        SUMD = 0.0D0
        SUMT = 0.0D0

        DO ISUM = 1,8

          ALOC = ((AUPR-ALWR) * XX08(ISUM) + AUPR+ALWR) / 2.0D0
          WAL  = (AUPR-ALWR)  * WW08(ISUM) / 2.0D0
          RLOC = CONL / TEXO / ALOC
          CALL CORONA(RLOC,DKIN,TKIN)
          SUMD = SUMD + WAL * DKIN * (RLOC**3) / ALOC
          SUMT = SUMT + WAL * DKIN * TKIN * (RLOC**3) / ALOC

        END DO

        DNOM = BRAD(II+1)**3 - BRAD(II)**3
        HDNS(II)  = 3.0D0 * SUMD/DNOM
        O2DNS(II) = 0.0D0
        TATM(II)  = SUMT / SUMD                                !! weighted mean

      END DO

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

      RETURN
      END

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