      SUBROUTINE CO2NIR(HSCO2,HR,RMU0,
     1                  XVF,CO2VF,PCD,ICD,K1CO2S,K2CO2S,WCO2S,WCO2ST,
     2                  CCO2,LAY,ILG,IL1,IL2)
C****************************************************************************
C*                                                                          *
C*                    SUBROUTINE CO2NIR                                     *
C*                                                                          *
C****************************************************************************
C
C to calculate solar heating in the Near IR CO2 bands. Method by
C Ogibalov and Fomichev [ASR, 2002] is used.
C
C                                    September, 2002, V.I. Fomichev
C              - modified for the CKD scheme, Novrmber, 2004, V.I. Fomichev
C called by MAMRAD2
C Calls A18LIN
C
C OUTPUT: HSCO2 -- solar heating in the Near IR bands in K/s at model grid
C ^^^^^^ 
C INPUT:  RMU0 -- cosine of the solar zenith angle;
C ^^^^^   CCO2 -- CO2 vmr in the troposphere
C         xvf  -- model vertical grid in pressure scale height
C                         (passed through THERDAT, calculated by THERMDAT2)
C         CO2vf-- CO2 v.m.r. (passed through THERDAT, calculated by vmrCO2,
C                             which is called from MAMRAD2)
C         CO2CL-- CO2 column amount at XR grid (0-16.5,0.25) (passed through
C                                 PIRCO2, determined by DPMCO2 subroutine,
C                                 based on vmr from CO2VP array)
C TCO2IR, QNIR -- parameterization coefficients (passed through CO2NIR,
C                                   determined in PCO2NIR subroutine)
C     PCD, ICD -- coefficients and indexes for interpolation from radiation
C                 to model vertical grid (passed through VFCLPR, calculated
C                 by PRECL2 which is called after THERMDAT2)
C K1CO2S,K2CO2S, WCO2S,WCO2ST - indexes and weights to implement the heating
C                           (passed through RADIDX, calculated in THERMDAT2)

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

      REAL HSCO2(ILG,LAY), RMU0(ILG)

      REAL  , DIMENSION(LAY)   :: XVF,CO2VF
      REAL  , DIMENSION(3,LAY) :: PCD
      REAL  , DIMENSION(LAY)   :: WCO2S,WCO2ST
      INTEGER K1CO2S, K2CO2S
      INTEGER,DIMENSION(LAY)   :: ICD

      COMMON /PIRCO2/ CO2VP(67), CO2CL(67)
      COMMON /CO2NIRP/ TCO2IR(49,10),QNIR(49,10)

c auxiliary arrays
      real HR(ILG,49), tCO2o(10), qCO2o(10)
c============================================================================
c initialization
      do k=1,LAY
      do il = IL1,IL2
      HSCO2(il,k)=0.0
      end do
      end do

      do k=1,49
      do il = IL1,IL2
      HR(il,k)=1.0
      end do
      end do

C ***** HEATING AT RADIATION VERTICAL GRID (x=2-14,0.25):

      DO 10 k=1,49

c parameterization coefficients for level "k"
      do i=1,10
      tCO2o(i)=TCO2IR(k,i)
      qCO2o(i)=QNIR(k,i)
      end do
c column CO2 amount above level "k"
      tCO2=CO2CL(k+8)*CCO2/3.600E-04

      DO 11 il=IL1, IL2

        heatCO2=0.0
        IF (RMU0(IL) .gt. 0.0) THEN

c CO2 amount along the paths
        g=log(tCO2*35./SQRT(1224.0*RMU0(IL)*RMU0(IL)+1.))

c solar heating:
        if(g.lt.tCO2o(1)) then
           heatCO2=qCO2o(1)
        else if (g.gt.tCO2o(10)) then
           heatCO2=0.0
        else
           heatCO2=a18lin(g,tCO2o,qCO2o,1,10)
        end if

        END IF      ! RMU0

        HR(il,k)=heatCO2

   11 CONTINUE
   10 CONTINUE

C ***** HEATING AT MODEL GRID

      DO 20 k=K1CO2S,K2CO2S-1
      kk = ICD(k)
        do 21 il = IL1,IL2
        HSCO2(il,k) = CO2vf(k)*(PCD(1,k)*HR(il,kk)+
     +                PCD(2,k)*HR(il,kk+1)+PCD(3,k)*HR(il,kk+2))
   21   continue
   20 CONTINUE
      DO 22 k=K2CO2S,LAY 
          do 23 il = IL1,IL2
          HSCO2(il,k) = WCO2ST(k)*CO2vf(k)*HR(il,48)
   23     continue
   22 CONTINUE

      return
      end
