      SUBROUTINE PCOOLG2(H,T,CO2,O2,SN2,O,AM,ZCO2O,FLUX,
     &                   SU,lambda,H2,ILG,IL1,IL2)

C****************************************************************************
C*                                                                          *
C*                       SUBROUTINE  PCOOLG2                                *
C*                                                                          *
C****************************************************************************
C
C  Modified from PCOOLG to be used with CKD scheme:
C    (1) 9.6 um O3 band is taken out (results from the CKD scheme are used)
C                                   July, 2004. V. I. Fomichev
C
C  to calculate the heating rate in both, 15 um CO2 and 9.6 um O3 bands.
C  NLTE conditions are taken into account for CO2 band. Irregular vertical
C  grid, to account for heat exchange in the "LTE-layer" (x=2-12.5), is used.
C  The subroutine can be used for an arbitrary CO2 volume mixing ratio
C  with the parameterization coefficients calculated by the CCO2GR.
C
C                                    November, 1996.  V.I. Fomichev

C Called by MAMRAD2
C Calls nothing

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

C input:
C T(longitude,67)  - temperature (K) at x(=0-16.5, 0.25) levels
C CO2,O2,N2,O(17) - volume mixing ratios (vmr's) for corresponding gases
C                   at X = 12.5 - 16.5, with a step of 0.25
C AM(17) - air molecular weight at X = 12.5 - 16.5, 0.25
C ZCO2O    - collisional deactivation rate constant (in cm3/sec)
C ILG = IL1, IL2 - longitude indexes
C parameterization coeficients for CO2 (calculated by CCO2GR) come from
C the CO2CFG and PIRO3 common blocks.
C the following data are defined in the BLOCK DATA DPMCO2:
C irregular vertical grid (ig(9)) comes from the PIRMGR common block
C pressure scale height arrays (0(2)-16.5,0.25) - PIRGRD common block
C some constants for cooling rate calculation - PIRCONS common block

C output:
C  H(ILG,59) - heating rates in erg/g/sec at X = 2-16.5, 0.25
C  FLUX(ILG) - upward flux at X=16.5 to calculate the heating rate
C              above this level.
C
C Note: 1) Actually, any arbitrary vertical grid could be utilized above
C          from x=13.75. To do so, the coefficients for the reccurence
C          formula must be determined at a proper grid (using CCO2GR) and all
C          input information should be given at this grid too.
C       2) As of now, ZCO2O is recommended to vary in a range of (1.5-6)e-12
C          with a mean value of 3e-12 cm3/sec, without T-dependence
C       3) To get the heating rate values in K/sec, the heating rates
C          calculated here should be divided by Cp - specific heat at constant
C          pressure. Where Cp = R/AM * {7/2*SUM[Rv,i - for 2-atom gases] +
C          5/2*SUM[Rv,i - for 1-atom gases]}, R - gas constant, Rv - vmr's.

      REAL H(ILG,59),FLUX(ILG),T(ILG,67),
     &                         CO2(17),O2(17),SN2(17),O(17),AM(17)

      common /CO2CFG/ amat(43,9),bmat(43,9),al(17)
      common /PIRMGR/ ig(9)
      common /PIRO3/  ao3(35,9)
      common /PIRGRD/ x(67),xh(59)
      common /PIRCONS/ const, constb, boltz, a10, aku

C INTERNAL arrays:

C SU(ILG,67)  - exponential part of the Planck function for CO2 band
C               at x(67) grid
C lambda(ILG,17) - quantum survival probability at the X grid for levels
C                  above x=12.5 (indexes 51-67)

      real SU(ILG,67), lambda(ILG,17)

c auxiliary arrays
C
      real H2(ILG)
C

c to determine SU(ILG,67)

      do 3 i = 1,67
      do 4 il = IL1, IL2
        SU(il,i)=exp(aku/T(il,i))
    4 continue
    3 continue

c to determine lambda

      do 5 i=1,17
      do 6 il = IL1, IL2

c ----- CO2-O2 and CO2-N2 V-T constants:
        tt=T(il,i+50)
        y=tt**(-1./3.)
        zn2=5.5e-17*sqrt(tt)+6.7e-10*exp(-83.8*y)
        zo2=1.e-15*exp(23.37-230.9*y+564.*y*y)

c ----- air number density:
        d=1.e6*exp(-x(i+50))/(tt*boltz)

c ----- collisional deactivation rate:
        z=(SN2(i)*zn2+O2(i)*zo2+O(i)*ZCO2O)*d

        lambda(il,i) = a10/(a10+z)
    6 continue
    5 continue

c cooling rate in the 15 um CO2 band (x=2-12.5, matrix approach is used)

      do 7 i=1,5
      is = i+8
      do 8 il = IL1, IL2
        H(il,i)=(amat(i,1)+bmat(i,1)*SU(il,is))*SU(il,1)
    8 continue
      do 9 j=3,9
      js=is+ig(j)
      do 10 il = IL1, IL2
        H(il,i)=H(il,i)+(amat(i,j)+bmat(i,j)*SU(il,is))*SU(il,js)
   10 continue
    9 continue
    7 continue

      do 12 i=6,18
      is = i+8
      do 13 il = IL1, IL2
        H(il,i)=(amat(i,1)+bmat(i,1)*SU(il,is))*SU(il,1)
   13 continue
      do 14 j=2,9
      js=is+ig(j)
      do 15 il = IL1, IL2
        H(il,i)=H(il,i)+(amat(i,j)+bmat(i,j)*SU(il,is))*SU(il,js)
   15 continue
   14 continue
   12 continue

      do 17 i=19,43
      is = i+8
      do 18 il = IL1, IL2
        H(il,i)=0.
   18 continue
      do 19 j=1,9
      js=is+ig(j)
      do 20 il = IL1, IL2
        H(il,i)=H(il,i)+(amat(i,j)+bmat(i,j)*SU(il,is))*SU(il,js)
   20 continue
   19 continue
   17 continue

c calculate the heating rate for x=12.75-16.5 (the reccurence formula is used)

c --- to form the boundary condition at x=12.5
      do 26 il = IL1, IL2
      H2(il)=H(il,43)/(CO2(1)*(1.-lambda(il,1))*constb)
   26 continue

c --- the reccurence formula
      do 27 I=2,17
      im=i-1
      d1=-.25*(AL(i)+3.*AL(im))
      d2=.25*(3.*AL(i)+AL(im))
      do 28 il = IL1, IL2
        hh=((1.-lambda(il,im)*(1.+d1))*H2(il)-d1*SU(il,im+50)-
     *                     d2*SU(il,i+50))/(1.-lambda(il,i)*(1.-d2))
        H2(il)=hh
        H(il,i+42)=hh*CO2(i)*(1.-lambda(il,i))/AM(i)*const
   28 continue
   27 continue

c to determine FLUX
c cooling rate above x=16.5 is suggested to be calculated by equation
c           H(i) = const/AM(i)*CO2(i)*(1.-lambda(i))*(FLUX-SU(i))
c no any parameterization coefficients are needed and an arbitrary hight grid
c can be used above x=16.5 level

      do 29 il = IL1, IL2
      FLUX(il) = H2(il) + SU(il,67)
   29 continue

      return
      end
