      SUBROUTINE THERMDAT2(xvf, pvf, O2vf, CN2vf, Ovf, AMuvf, 
     1                     Cpvf, Gvf, 
     2                     effH, effSRC, tO2, effO2nir,
     3                     X1chem, X2chem, DXchem, WCHEM, K1CHEM,
     4                     WWVS2, WWVS3, WCO2S, WCO2ST,
     5                     WWVL3, WCH4L, WO3L, K1CO2S, K2CO2S,
     6                     ETAG, LAY)

C***********************************************************************
C*                                                                     *
C*              SUBROUTINE THERMDAT2                                   *
C*                                                                     *
C***********************************************************************

C modified from THERMDAT to be used with the CKD scheme
C        V. Fomichev, November, 2004
C
c 1) to define input data on concentrations of O, O2, N2, as well as data
c    on the air molecular weight and specific heat to be used in the
c    radiative scheme for the lower thermosphere (taken from the MSIS
c    global average model;
c 2) to calculate solar heating efficiency in Hartley (Ha) O3 band,
c    Schumann-Runge O2 continuum (SRC) (SRBands efficiency = 1).
c    Method by Mlynczack&Solomon (JGR, vol 98, p 10517, 1993) is used;
c 3) to calculate O2 column number density at the model mid-layer levels
c    used for calculation of the O2 solar heating in SRC and SRB.
c 4) to define parameters used for the chemical heating calculation.
c 5) to define indexes and weights used to smoothly indroduce the
c    upper atmosphere radiation scheme.
c 6) to estimate the solar heating efficiency in O2 bands: .762, .688, .629
c    and 1.27 micron. Approximation based on results shown on Fig. 1 of
c    Mlynczak&Marshall (GRL, 1996, vol. 23, p. 657). Efficiency for .762 um
c    is used for its giving the dominating contribution.
c
c       V. Fomichev, November, 1997
c                    October,  1999 - modified by inclusion (4) - VF.
c                    November, 2004 - (5) and (6) included
c
C Called by model driver outside of parallel region
C Calls A18LIN

      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)

c INPUT:  ETAG(LAY) - model hybrid levels [p/p(surface)] for the mid-layers
c                     levels ORDERED FROM THE TOP (1 - moon layer) DOWN.

      REAL ETAG(LAY)

c OUTPUT: all information is giving at log-pressure grid defined for
c         "mid-layer" levels. ORDERED FROM THE SURFACE UP.
c
c  Output results are passed throught the following common blocks in
c  the driver (fields are passed out to the driver):
c  /THERDAT/:
c         pvf - pressure (mb); xvf - log-pressure height;
c         O2vf,CN2vf,Ovf - vmr's for O2, N2, and O;
c         AMuvf - air molecule weight (g/mol);
c         Cpvf -specific heat (erg/g/K); Gvf - gravity acceleration (m s-2)
c         CO2vf - vmr for CO2 is determined in the vmrCO2 subroutine
c  /THERSH/:
c  effH - solar heating efficiency in the Hartley O3 band
c  effSRC - solar heating efficiency for SRC
c  tO2 - O2 column number density (cm-2) at the model mid-layer levels
c  effO2nir - solar heating efficiency in the "NIR" O2 bands (.629-1.27um)
c  /CHHDAT/:
c  X1chem, X2chem, DXchem - merging region boundaries (to implement chem.heat)
c  K1CHEM, WCHEM - model index for X1chem and weights at x > X2chem
c  /RADIDX/:
c indexes and weights from this block are used to smoothly indroduce the
c upper atmosphere radiation scheme (similar to K1CHEM, WCHEM)
c WWVS2 - weights for H2O heating in 4200-8400 cm-1
c WWVS3 - weights for H2O heating in 2500-4200 cm-1
c K1CO2S,K2CO2S,WCO2S,WCO2ST - indexes and weights for the NIR CO2 bands
c WWVL3 - weights for H2O IR heating in 1400-1900 cm-1
c WCH4L - weights for CH4 IR heating in 1100-1400 cm-1
c WO3L  - weights for O3 IR heating in 980-1100 cm-1 (9.6um)

      real xvf(lay), pvf(lay), O2vf(lay), CN2vf(lay),
     *     Ovf(lay), AMuvf(lay), Cpvf(lay), Gvf(lay)
      real effH(lay), effSRC(lay), tO2(lay), effO2nir(lay)
      real X1chem, X2chem, DXchem, WCHEM(lay)
      real WWVS2(lay), WWVS3(lay), WCO2S(lay), WCO2ST(lay),
     *     WWVL3(lay), WCH4L(lay), WO3L(lay)
      integer K1CHEM, K1CO2S, K2CO2S

c initial input info is taken from the MSIS global average model:
c Muo(g) - air mol. weight; Oo,N2o,O2o - vmr; Cpo(J/kg) - specific heat;
c Go(m s-2) gravity acceleration.
      real xo(45), Muo(45), Oo(45), N2o(45), O2o(45), Cpo(45), Go(45)
      real Muoo(45), Ooo(45), N2oo(45), O2oo(45), Cpoo(45), Goo(45)

C Efficiency for 762 nm O2 (taken from Mlynczak&Marshall, GRL, 1996)
      real xO2nir(10), effO2o(10)

C Coefficients for the Chebyshev polynomial fit for the heating efficiency
      real ChO3l(4), ChO3h(4), ChO2(4)
C ----- for the Ha band:
      data ChO3l/0.926210, 0.133960,-0.076863, 0.006897/
      data ChO3h/0.669650,-0.009682, 0.033093, 0.017938/
C ----- for the SRC:
      data ChO2/0.75349, 0.0036, 0.059468, -0.022795/

C ----- efficiency for 762 nm O2 (taken from Mlynczak&Marshall, GRL, 1996):
      data (xO2nir(i), i=1,10)/
     & 4.44, 5.84, 7.12, 8.44, 9.90,11.50,13.24,15.07,16.62,17.18/
      data (effO2o(i), i=1,10)/
     & 1.,  0.992,0.976,0.919,0.732,0.455,0.301,0.175,0.050,0.0/

c Avogadro number * 10 (only for tO2 calculation)
      data aNa/6.022045E24/

C ----- MSIS data:
      data (xo(i), i=1,45)/
     &  0.0,     0.5,     1.0,     1.5,     2.0,     2.30000, 2.55259,
     &  3.05259, 3.55259, 4.05259, 4.55259, 5.05259, 5.55259, 6.05259,
     &  6.55259, 7.05259, 7.55259, 8.05259, 8.55259, 9.05259, 9.55259,
     & 10.05260,10.55260,11.05260,11.55260,12.05260,12.55260,13.05260,
     & 13.55260,14.05260,14.55260,15.05260,15.55260,16.05260,16.55260,
     & 17.05260,17.55260,18.05260,18.55260,19.05260,19.55260,20.05260,
     & 20.55260,21.05260,21.30000/
      data (Muo(i), i=1,45)/
     & 28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,
     & 28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,28.95,
     & 28.95,28.95,28.94,28.94,28.94,28.93,28.92,28.89,28.83,28.73,
     & 28.58,28.38,28.13,27.80,27.39,26.89,26.30,25.65,24.97,24.27,
     & 23.53,22.73,21.89,21.03,20.60/
      data (Oo(i), i=1,45)/
     & 0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     & 0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     & 0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     & 0.000E+00,0.000E+00,0.000E+00,0.240E-08,0.127E-06,0.176E-05,
     & 0.177E-04,0.125E-03,0.615E-03,0.217E-02,0.575E-02,0.122E-01,
     & 0.219E-01,0.350E-01,0.520E-01,0.736E-01,0.101E+00,0.135E+00,
     & 0.177E+00,0.224E+00,0.275E+00,0.329E+00,0.388E+00,0.452E+00,
     & 0.519E+00,0.589E+00,0.623E+00/
      data (N2o(i), i=1,45)/
     & 0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,
     & 0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,0.781,
     & 0.782,0.782,0.782,0.783,0.783,0.784,0.785,0.786,0.786,0.785,
     & 0.782,0.778,0.774,0.767,0.758,0.745,0.726,0.698,0.662,0.618,
     & 0.568,0.512,0.452,0.389,0.357/
      data (O2o(i), i=1,45)/
     & 0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.210,
     & 0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.210,0.209,
     & 0.209,0.209,0.209,0.208,0.208,0.207,0.205,0.203,0.200,0.195,
     & 0.188,0.179,0.168,0.153,0.136,0.116,0.094,0.076,0.061,0.051,
     & 0.043,0.036,0.029,0.022,0.020/
      data (Cpo(i), i=1,45)/
     & 1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,
     & 1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,1002.5,
     & 1002.5,1002.6,1002.6,1002.7,1002.8,1002.9,1003.0,1003.2,1003.5,
     & 1004.1,1005.2,1006.9,1009.4,1012.8,1017.3,1023.0,1030.3,1039.2,
     & 1049.7,1061.3,1073.3,1085.8,1099.4,1114.8,1132.1,1151.0,1160.9/
      data (Go(i), i=1,45)/
     * 9.8,9.79,9.78,9.77,9.765,9.76,9.75,9.74,9.73,9.72,9.71,9.70,
     * 9.69,9.68,9.67,9.66,9.64,9.63,9.62,9.61,9.60,9.59,9.58,9.57,
     * 9.56,9.55,9.55,9.54,9.53,9.52,9.51,9.51,9.50,9.49,9.48,9.46,
     * 9.45,9.43,9.40,9.37,9.33,9.29,9.23,9.18,9.15/

c ----- to form THERDAT common block:
c       note: CO2vf is determined in the vmrCO2 subroutine

      do 1 i = 1,LAY
       pvf(i)  = ETAG(LAY-i+1)*1000.0
       xvf(i) = alog(1000/pvf(i))
    1 continue

      do 2 i = 1, 45
      Muoo(i) = alog(Muo(i))
      Ooo(i)  = alog(Oo(i)+1.)
      N2oo(i) = alog(N2o(i))
      O2oo(i) = alog(O2o(i))
      Cpoo(i) = alog(Cpo(i))
      Goo(i)  = alog(Go(i))
    2 continue

      do 3 i = 1, LAY
      AMuvf(i) = exp(a18lin(xvf(i),xo,Muoo,1,45))
      Ovf(i)   = exp(a18lin(xvf(i),xo,Ooo,1,45)) - 1.
      CN2vf(i) = exp(a18lin(xvf(i),xo,N2oo,1,45))
      O2vf(i)  = exp(a18lin(xvf(i),xo,O2oo,1,45))
c Cpvf in (erg/g/K)
      Cpvf(i)  = exp(a18lin(xvf(i),xo,Cpoo,1,45))*10000.
      Gvf(i)   = exp(a18lin(xvf(i),xo,Goo,1,45))
    3 continue

c ----- to calculate effH
c efficiency in the O3 Hartley band. The efficiency coefficient should be
c applied to the NET (i.e., total - chemical potential) heating.

      do 10 k=1,LAY
      p = pvf(k)
      if(p.gt.1.) then
        effH(k)=1.
        else if(p.gt.1.e-2) then
          x = alog10(p)+1
          effH(k)=ChO3l(1)+ChO3l(2)*x+ChO3l(3)*x*x+ChO3l(4)*x*x*x
          else if(p.gt.1.e-4) then
            x = alog10(p)+3
            effH(k)=ChO3h(1)+ChO3h(2)*x+ChO3h(3)*x*x+ChO3h(4)*x*x*x
            else
              effH(k) = 0.6944
      end if
   10 continue

c ----- to calculate the SRC heating efficiency (for SR band efficiency = 1.)

      do 20 k = 1, LAY
      p = pvf(k)
      if(p.gt.1.e-2) then
        effSRC(k) = 0.7938
        else if(p.gt.1.e-4) then
          x = alog10(p)+3
          effSRC(k) = ChO2(1) + ChO2(2)*x + ChO2(3)*x*x + ChO2(4)*x*x*x
          else
            effSRC(k) = 0.8320
      end if
   20 continue

c ----- to define the O2 heating efficiency for 762 nm

      do 30 k = 1, LAY
      xk=xvf(k)
      if (xk.le.xO2nir(1)) then
        effO2nir(k) = 1.
      else if (xk.le. xO2nir(10)) then
        effO2nir(k) = a18lin(xk,xO2nir,effO2o,1,10)
      else
        effO2nir(k) = 0.0
      end if
   30 continue

c ----- to calculate tO2(cm-2)

      g1=O2vf(LAY)*aNa*pvf(LAY)/(AMuvf(LAY)*Gvf(LAY))
      tO2(LAY)=g1*AMuvf(LAY)/32.

      kmxm = LAY-1
      do 25 kk=1,kmxm
      k=LAY-kk
      g2=O2vf(k)*aNa*pvf(k)/(AMuvf(k)*Gvf(k))
      tO2(k)=tO2(k+1)+(g1+g2)*0.5*(xvf(k+1)-xvf(k))
      g1=g2
   25 continue

c ----- to form CHHDAT common block:
c       merging region for chemical heating and O2 (SRC&SRB) O2 heating

      X1chem = 8.3
      X2chem = 9.8
      DXchem = X2chem - X1chem

      do 4 k=1,LAY
      if (xvf(k).ge.X1chem) goto 5
    4 continue
    5 K1CHEM=k
      do 6 k=1,LAY
       if (xvf(k).le.X1chem) then
          WCHEM(k)=0.
       else if (xvf(k).le.X2chem) then
          WCHEM(k)=(xvf(K)-X1chem)/DXchem
       else
          WCHEM(k)=1.
       end if
    6 continue 

c ----- to form RADIDX common block:
c       indexes and weights to merge the upper atmosphere rad. scheme.

c * SOLAR H2O:
      X1WVS2 = 2.0
      X2WVS2 = 3.0
      X1WVS3 = 7.5
      X2WVS3 = 12.25
c * weights:
      do k=1,LAY
       if (xvf(k).le.X1WVS2) then
          WWVS2(k)=1.
       else if (xvf(k).le.X2WVS2) then
          WWVS2(k)=(X2WVS2-xvf(K))/(X2WVS2-X1WVS2)
       else
          WWVS2(k)=0.
       end if
       if (xvf(k).le.X1WVS3) then
          WWVS3(k)=1.
       else if (xvf(k).le.X2WVS3) then
          WWVS3(k)=(X2WVS3-xvf(K))/(X2WVS3-X1WVS3)
       else
          WWVS3(k)=0.
       end if
      end do

c * NIR CO2:
      X1CO2S = 7.0
      X2CO2S = 8.5
      X3CO2S = 13.75
      X4CO2S = 15.50
c * indexes:
      do k=1,LAY
      if (xvf(k).ge.X1CO2S) goto 11
      end do
   11 K1CO2S=k
      do k=K1CO2S,LAY
      if (xvf(k).ge.X3CO2S) goto 12
      end do
   12 K2CO2S=k
c * weights:
      do k=1,LAY
       if (xvf(k).le.X1CO2S) then
         WCO2S(k) =1.
         WCO2ST(k)=1.
       else if (xvf(k).le.X2CO2S) then
         WCO2S(k) =(X2CO2S-xvf(K))/(X2CO2S-X1CO2S)
         WCO2ST(k)=1.
       else if (xvf(k).le.X3CO2S) then
         WCO2S(k) =0.
         WCO2ST(k)=1.
       else if (xvf(k).le.X4CO2S) then
         WCO2S(k) =0.
         WCO2ST(k)=(X4CO2S-xvf(K))/(X4CO2S-X3CO2S)
       else
         WCO2S(k) =0.
         WCO2ST(k)=0.
       end if
      end do

c * H2O IR BANDS:
      X1WVL3 = 11.5
      X2WVL3 = 14.0
c * weights:
      do k=1,LAY
       if (xvf(k).le.X1WVL3) then
          WWVL3(k)=1.
       else if (xvf(k).le.X2WVL3) then
          WWVL3(k)=(X2WVL3-xvf(K))/(X2WVL3-X1WVL3)
       else
          WWVL3(k)=0.
       end if
      end do

c * CH4 IR BAND:
      X1CH4L = 10.5
      X2CH4L = 13.0
c * weights:
      do k=1,LAY
       if (xvf(k).le.X1CH4L) then
          WCH4L(k)=1.
       else if (xvf(k).le.X2CH4L) then
          WCH4L(k)=(X2CH4L-xvf(K))/(X2CH4L-X1CH4L)
       else
          WCH4L(k)=0.
       end if
      end do

c * O3 9.6um BAND:
      X1O3L = 10.0
      X2O3L = 15.0
c * weights:
      do k=1,LAY
       if (xvf(k).le.X1O3L) then
          WO3L(k)=1.
       else if (xvf(k).le.X2O3L) then
          WO3L(k)=(X2O3L-xvf(K))/(X2O3L-X1O3L)
       else
          WO3L(k)=0.
       end if
      end do

      return
      end
