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

C VERTICAL COLUMN DENSITIES BASED ON "FINE-GRID" MODEL ATMOSPHERE
C [29Jan02 JBishop]

C -----------------------------------------------------------------------------
C  TCOLM_H  . . . atomic hydrogen column density up to exobase
C  TCOLM_O2 . . . O2 column density up to exobase
C  TCOLM_TOT  . . total "effective opacity" optical depth, i.e., sum of
C                 TCOLM_H and TCOLM_O2 scaled by optical cross sections ratio
C  COLM_EXO . . . atomic hydrogen column density from exobase on
C -----------------------------------------------------------------------------

      SUBROUTINE COLM_DEN(ITHERM,THERMO,
     &                    TCOLM_H,TCOLM_O2,TCOLM_TOT,COLM_EXO)

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

      DIMENSION THERMO(5,61)

      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 /GAUS16/WW16(16),XX16(16)
      COMMON /GAUS32/WW32(32),XX32(32)

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

      TCOLM_H  = 0.0D0
      TCOLM_O2 = 0.0D0

      DO ISUM = 1,ITHERM

        DZ  = THERMO(2,ISUM+1) - THERMO(2,ISUM)

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

        H1H = DZ / LOG(H1N/THERMO(4,ISUM+1))
        O2H = DZ / LOG(O2N/THERMO(5,ISUM+1))

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

        TCOLM_H  = TCOLM_H  + H1C
        TCOLM_O2 = TCOLM_O2 + O2C

      END DO

      TCOLM_TOT = TCOLM_H + (ABSCSX/CENTER) * TCOLM_O2

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

      AUPR  = CONL / TEXO / RC
      ALWR  = CONL / TEXO / RUPR

      COLM_EXO = 0.0D0

      DO ISUM = 1,16

        ALOC = ((AUPR-ALWR) * XX16(ISUM) + AUPR+ALWR) / 2.0D0
        WAL  = (AUPR-ALWR)  * WW16(ISUM) / 2.0D0
        RLOC = CONL / TEXO / ALOC

        CALL CORONA(RLOC,DKIN,TKIN)

        COLM_EXO = COLM_EXO + WAL * DKIN / ALOC**2

      END DO

      COLM_EXO = COLM_EXO * CONL / TEXO

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

      RETURN
      END

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