      SUBROUTINE MAMRAD4(HRS, HRL, 
     1                   HRSH2O, HRSO3, HRSCO2, HRSO2,
     2                   HRLH2O, HRLO3, HRLCO2, HRLCH4, HRLN2O, 
     3                   T, RMU,
     4                   xvf, pvf, o2vf, cn2vf, ovf, amuvf, 
     5                   cpvf, gvf, co2vf,
     6                   effh, effsrc, to2, effo2nir,
     7                   x1chem, x2chem, dxchem, wchem, k1chem,
     8                   wwvs2, wwvs3, wco2s, wco2st,
     9                   wwvl3, wch4l, wo3l, k1co2s, k2co2s,
     A                   pdr, pcd, idr, icd,
     B                   CO2_PPM, LAY,ILG,IL1,IL2)

C****************************************************************************
C*                                                                          *
C*                       SUBROUTINE MAMRAD4 
C*                                                                          *
C****************************************************************************
C
C to merge/adapt the original CKD scheme with/to the MA schemes/conditions
C
C                                   Apr, 2011 - M. Lazare (adapted for gcm16).                                     
C                                   July, 2004, V. I. Fomichev

C Called by physics driver.
C Calls CO2NIR, PCOOLG2

C INPUT
C ^^^^^   (all the CKD arrays are ordered from the top down):
C
C HRSH2O(ILG,LAY,3) - solar heating in the H2O bands.
C                     1: 8400-14500 (+contribution from O3 in this band);
C                        set to 0 above 100 hPa
C                     2: 4200-8400; 3: 2500-4200. 
C HRSO3(ILG,LAY,9) - solar heating in the O3 bands.
C                    1: 42000-50000; 2: 37400-42000; 3: 35700-37400;
C                    4: 34000-35700; 5: 32185-34000; 6: 30300-32185;
C                    7: 25000-30300; 8: 20000-25000; 9: 14500-20000.
C                      (+contribution from H2O below 100hPa in this band)
C HRSCO2(ILG,LAY) - solar heating in the CO2 near-IR bands. Contributions
C                   from intervals (4200-8400)+(2500-4200).
C HRSO2(ILG,LAY,2) - solar heating due to O2.
C                    1: 42000-50000; 2: (14500-20000)+(8400-14500).
C
C HRLH2O(ILG,LAY,9) - IR heating in H2O bands(cm-1). 
C                     1: 2200-2500 (+contribution from N2O in this band); 
C                     2: 1900-2200 (+N2O); 3: 1400-1900; 4: 1100-1400 (+CFCs);
C                     5: 980-1100 (+CO2+CFCs); 6: 800-980 (+CFCs)
C                     7: 540-800 (N2O+O3+H2O are included in HRLCO2(x,x,2));
C                     8: 340-540; 9: 0-340.
C HRLO3(ILG,LAY) - IR heating in the 9.6um O3 band (interval: 980-1100)
C HRLCO2(ILG,LAY,2)- IR heating in CO2 bands.
C                    1: 2200-2500 (4.3um band);
C                    2: 540-980 (15um band, plus N2O+O3+H2O in 540-800).
C HRLCH4(ILG,LAY) - IR heating in 7-9um CH4 band (interval 1100-1400).
C HRLN2O(ILG,LAY) - IR heating in 7-9um N2O band (interval 1100-1400).
C
C Input arrays formed in THERMDAT2, VMRCO2, PRECL2, and DPMCO2 (surface --> top):
C
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)
C  /VFCLPR/:
C     PDR(3,67)- coefficients for interpolation from X(N) to XR(67) grid
C     IDR(67)- the indexes of levels (for X(N)) for second order
C              interpolation from X(N) to XR(67) grid
C     PCD(3,N),ICD(N) - coefficients and indexes of levels for
C                       interpolation from XC(59) to X(N)
C
C T(ILG,LAY) - temperature (K) for mid-layer levels.
C CO2_PPM - CO2 vmr in the troposphere
C RMU(ILG) - cosine of the zenith angle
C LAY - number of levels, ILG - number of longitudes (from IL1 to IL2)
C
C OUTPUT
C ^^^^^^ (should be ordered similar to the input, i.e. top -> surface)
C HRS, HRL(ILG,LAY) - total solar and IR heating from CKD and the MA.

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

      REAL HRS(ILG,LAY), HRL(ILG,LAY)
      REAL HRSH2O(ILG,LAY,3), HRSO3(ILG,LAY,9),
     &     HRSCO2(ILG,LAY), HRSO2(ILG,LAY,2)
      REAL HRLH2O(ILG,LAY,9), HRLO3(ILG,LAY), HRLCO2(ILG,LAY,2),
     &     HRLCH4(ILG,LAY), HRLN2O(ILG,LAY)
c
      real  , dimension(3,lay) :: pcd
      real  , dimension(3,67)  :: pdr
      real  , dimension(lay)   :: xvf,pvf,o2vf,cn2vf,ovf,amuvf,
     1                            cpvf,gvf,co2vf,effh,effsrc,to2,
     2                            effo2nir
      real  , dimension(lay)   :: wwvs2,wwvs3,wco2s,wco2st,wwvl3,
     1                            wch4l,wo3l,wchem
      integer,dimension(lay)   :: icd
      integer,dimension(67)    :: idr 
      real x1chem, x2chem, dxchem
      integer k1chem, k1co2s, k2co2s
c
      REAL T(ILG,LAY), RMU(ILG)
      INTEGER LAY, ILG, IL1,IL2
 
C INTERNAL ARRAYS:

C Chemical potential energy generated by O3 and O2 (only Hz) photolysis
C                 (ordered top --> surface)

      REAL CHPOT(ILG,LAY)

C ARRAYS USED in the matrix parameterization for the the NIR CO2 bands

      REAL HSCO2(ILG,LAY), HSPCO2(ILG,49)

C ARRAYS USED in the matrix parameterization for the 15um CO2 bands

      REAL TPCO2(ILG,67), HPCO2(ILG,59), FLPCO2(ILG), SPCO2(ILG,67),
     &     ALPCO2(ILG,17), H2PCO2(ILG),
     &     CO2R(17),O2R(17),SN2R(17),OR(17),AMR(17) 

      common /PIRGRD/ xf(67), xhf(59)
C=====================================================================
C **************** SOLAR HEATING ********************************

C ---------------------------------------------------------------
C I. O3 HEATING IN 9 SPECTRAL INTERVALS (14500-50000) AND O2 IN 42000-50000
C

C Calculate the chemical potential energy and net heating (i.e. take out
C chemical potential energy and introduce efficiency in the Hartley O3 band).
C
C Part of the O3 and O2 heating converted to the chemical potential energy
C in each spectral interval with quantum energy of hv:
C      O3: po3=DO3/hv, DO3=1.01eV=1.666E-19 J 
C      02: p02=DO2/hv, DO2=5.12eV=8.203E-19 J 
C po3: 42000-50000: 0.9043
C po3: 42000-50000: 0.1837; 37400-42000: 0.2120; 35700-37400: 0.2296;
C      34000-35700: 0.2408; 32185-34000: 0.2536; 30300-32185: 0.2687;
C      25000-30300: 0.3061; 20000-25000: 0.3774; 14500-20000: 0.4989.

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
c O2 Hz band:
        HRS(I,K)   = HRSO2(I,K,1)*0.0957
        CHPOT(I,K) = HRSO2(I,K,1)*0.9043 
c O3 Ha band:
        O3HT = HRSO3(I,K,1)+HRSO3(I,K,2)+HRSO3(I,K,3)+
     &         HRSO3(I,K,4)+HRSO3(I,K,5)
        CHEM = HRSO3(I,K,1)*0.1837 + HRSO3(I,K,2)*0.2120 +
     &         HRSO3(I,K,3)*0.2296 + HRSO3(I,K,4)*0.2408 +
     &         HRSO3(I,K,5)*0.2536
        CHPOT(I,K) = CHPOT(I,K) + CHEM
        HRS(I,K) = HRS(I,K) + (O3HT - CHEM)*effH(KK)
c other (Hu&Sh) O3 bands:
        O3HT = HRSO3(I,K,6)+HRSO3(I,K,7)+HRSO3(I,K,8)+HRSO3(I,K,9)
        CHEM = HRSO3(I,K,6)*0.2687 + HRSO3(I,K,7)*0.3061 +
     &         HRSO3(I,K,8)*0.3774 + HRSO3(I,K,9)*0.4989
        CHPOT(I,K) = CHPOT(I,K) + CHEM
        HRS(I,K) = HRS(I,K) + (O3HT - CHEM)
        END DO
      END DO

c return chemical potential energy:
c Note: when chemical heating is implementing, introduce efficiency for
c       O+O2+M -> O3+M reaction (contribution of the O3 hot bands)
      DO K = 1, LAY
        DO I = IL1, IL2
        HRS(I,K) = HRS(I,K) + CHPOT(I,K)
        END DO
      END DO

C ---------------------------------------------------------------
C II. O2 HEATING IN SRB AND SRC (only include if the model lid above X1chem)
C

c Solar heating due to absorption in Schumann-Runge c continuum (SRC) and
c Schumann-Runge bands by the method of Strobel (JRG, c vol 83, p 6225, 1978)
c with taking into account the heating efficiency for SRC from
c Mlynczack&Solomon (JGR, vol 98, p 10517, 1993) (efficiency for SRB =1)
c Only NET (i.e. without chemical potential) heating is included.
c The SRC heating efficiency and O2 column amount is calculated in THERMDAT2
c and passed through common block THERSH (arrays: effSRC and tO2)

      IF(xvf(LAY).GT.X1chem) THEN

      DO 20 K = K1CHEM, LAY
c loop is done from the surface up
      KK = LAY+1-K
        DO 10 I = IL1, IL2
        heatO2 = 0.0
        if (RMU(I) .gt. 0.0) THEN
c  * O2 amount: downward looking paths:
            g=tO2(K)*35./SQRT(1224.0*RMU(I)*RMU(I)+1.)
c  * SRC:
            SRC1=2.716E6*exp(-1.E-17*g)
            SRC2=(5.902E23*exp(-2.9E-19*g)-3.312E23*exp(-1.7E-18*g)
     1           -2.590E23*exp(-1.15E-17*g))/g
            SRC = (SRC1+SRC2)*effSRC(K)
c  * SRB:
            SRB=1.463E5
            if(g.gt.1.e18) SRB=1./(1.113E-24*g+5.712E-15*sqrt(g))
c  * O2 heating (K/s):
            heatO2=O2vf(K)*(SRC+SRB)/AMuvf(K)/Cpvf(K)
c  * to include smoothly between X1chem and X2chem
            heatO2=WCHEM(K)*heatO2
        endif
c  * net (chemical potential energy is taken out) O2 heating is added
c  * to HRS(LAY,ILG) array.

        HRS(I,KK) = HRS(I,KK) + heatO2

   10   CONTINUE
   20 CONTINUE

      END IF

C ---------------------------------------------------------------
C III. H2O HEATING IN 3 SPECTRAL INTERVALS (2500-14500)
C

c to account for NLTE (based on the source function calculation) the
c following approximation used:
c 1. 8400-14500 (0.69-1.19um)- zero heating above 100 hPa is already forced
c 2. 4200-8400 (1.19-2.38um)- heating smoothly reduced to 0 between x=2-3
c 3. 2500-4200 (2.38-4.00um)- heating smoothly reduced to 0 between x=7.5-12.25
c   (the dominating/strongest bands in 4200-8400: 011, 101 -> 000, the first
c    one is under NLTE above ~5 km, the second one - above the ground, so we
c    include these bands only in the troposphere)
c   (the dominating/strongest band in 2500-4200 is 001-000, 2.7um. It's under
c    NLTE above ~40-50 km, the source function departs from the LTE by ~10 km
c    below than that of the 2.7um CO2 band. Calculations for the NIR CO2, 
c    1.05-4.3um (2.7um also dominates), show that LTE is valid up to x~8.1 and
c    at x=12.5 NLTE heating is ~0.1 of the LTE value. So, we introduced a
c    similar decrease but started lower and, so, reduced the LTE heating in
c    2500-4200 interval 0 between x=7.5-12.25).

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        H2OS=HRSH2O(I,KK,1)+HRSH2O(I,KK,2)*WWVS2(K)+
     &       HRSH2O(I,KK,3)*WWVS3(K)
        HRS(I,KK) = HRS(I,KK) + H2OS
        END DO
      END DO

C ---------------------------------------------------------------
C IV. O2 HEATING IN 2 SPECTRAL INTERVALS (8400-20000)
C

c to accout for NLTE the efficiency from Mlynczak&Marshall is used
c (GRL, 1996, vol. 23, p. 657, Fig. 1.)

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        O2HT=HRSO2(I,KK,2)*effO2nir(K)
        HRS(I,KK) = HRS(I,KK) + O2HT
        END DO
      END DO
 
C ---------------------------------------------------------------
C V. CO2 HEATING IN THE NEAR-IR BANDS (3 SPECTRAL INTERVALS: 2200-8400)
C         (note, interval 2200-2500 belongs to the longwave part)
C

c parameterization by Ogibalov and Fomichev (Adv.Space Res, 2003, p. 759)
c is used in the NLTE region above x=8.5. Between x=7.0-8.5 the heating
c rates are merged with the CKD results.

      CALL CO2NIR(HSCO2,HSPCO2,RMU,
     1            XVF,CO2VF,PCD,ICD,K1CO2S,K2CO2S,WCO2S,WCO2ST,
     2            CO2_PPM,LAY,ILG,IL1,IL2)

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        CO2HT=WCO2S(K)*HRSCO2(I,KK) + (1.-WCO2S(K))*HSCO2(I,K)
        HRS(I,KK) = HRS(I,KK) + CO2HT
        END DO
      END DO

C **************** LONGWAVE COOLING ********************************

C ---------------------------------------------------------------
C I. H2O+N2O in 2 INTERVALS: 2200-2500 (4-4.54um); 1900-2200 (4.54-5.26um)
C

c Main effect in the MA: (weak) heating due to absorption of solar radiation.
c H2O band 010-001 + 010-100 in 1900-2200 cm-1.
c Source function: 010 - NLTE above ~75 km; 001 - NLTE above ~40-50 km;
c  100 - NLTE as low as 3-5 km.
c N2O band 001-000. Source function: NLTE above ~55km.
c Weights WWVS3 (2500-4200 H2O) is used: reduced to 0 between x=7.5-12.25).

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        HRL(I,KK) = WWVS3(K)*(HRLH2O(I,KK,1)+HRLH2O(I,KK,2))
     &            + WCO2S(K)*HRLCO2(I,KK,1)
        END DO
      END DO

C ---------------------------------------------------------------
C II. H2O in 1400-1900 cm-1 (5.26-7.14 um)

c Main effect in the MA: (weak) cooling near the stratopause; (weak) heating
c  below 30 km and above 55 km (partly due to absorption of solar radiation).
c Dominating band 010-000 (6.3um). Source function: NLTE above ~75 km.
c Heating rate (Lopez-Puertas and Taylor, 2001): NLTE small below ~85 km and
c  very strong above 100 km, NLTE heating above 100 km is ~0.
c Reduce heating rates to 0 between x=11.5-14.

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        HRL(I,KK) = HRL(I,KK) + WWVL3(K)*HRLH2O(I,KK,3)
        END DO
      END DO

C ---------------------------------------------------------------
C III. H2O+CFS, CH4, N2O in 1100-1400 cm-1 (7.14-9.09um)

c Main effect in the MA: (weak) cooling near the stratopause; (weak) heating
c  below 25 km and above 60 km. CH4 dominates (especially above x=10).

c H2O: continuum + wings of 6.3um. (no NLTE study for continuum).
c CFS: no NLTE study.
c N2O bands 020-000; 100-000 (dominating). Source function: NLTE above ~75km 
c  for both day and night. 
c CH4 bands 0001-0000 (dominates); 0101-0100; 0002-0001. Source function 
c  for 0001: NLTE above 65-70 km.
c Reduce heating rate to 0 between
c   x=11.5-14 for H2O+CFS and N2O  (similar to 1400-1900);
c   x=10.5-13.0 for CH4 (by 1.0 lower than for 1400-1900).

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        H2OL = WWVL3(K)*(HRLH2O(I,KK,4)+HRLN2O(I,KK))+
     1         WCH4L(K)*HRLCH4(I,KK)
        HRL(I,KK) = HRL(I,KK) + H2OL
        END DO
      END DO


C ---------------------------------------------------------------
C IV. H2O+CO2+CFS, O3 in 980-1100 cm-1 (9.09-10.2um)

c Main effect in the MA: cooling near the stratopause; (weak) heating
c  below ~30 km and above 60 km. O3 dominates (the only contributor above x=8)
c H2O+CO2+CFS: no NLTE study.
c O3 bands 001-000 (9.6um, dominates); 100-000. Source functions for
c  fundamental bands under NLTE above 70-75 km without diurnal variation (note
c  strong diurnal variation for hot bands, but they are less important).
c Heating rate (Fomichev et al., 1994, Lopez-Puertas and Taylor, 2001): NLTE
c  above ~75 km.
c Reduce heating rates to 0 between
c  x=11.5-14 for H2O+CO2+CFS (similar to 1400-1900);
c  x=10-15 for O3. Based on F.et al., 1994 results without chemical pumping:
c      O3 heating =0 at ~x=15; ratio NLTE/LTE ~0.3 at x=13.5 - consistent
c      with L.-P.&T. (2001); starts to deviate from LTE at ~x=10). 

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        H2OL = WWVL3(K)*HRLH2O(I,KK,5)+WO3L(K)*HRLO3(I,KK)
        HRL(I,KK) = HRL(I,KK) + H2OL
        END DO
      END DO

C ---------------------------------------------------------------
C V. H2O+CFS in 800-980 cm-1 (10.2-12.5um) (CO2 contribution included later)

c Main effect in the MA: zero above ~x=1. no NLTE study.
c Reduce to 0 between x=11.5-14 (similar to 1400-1900).

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        HRL(I,KK) = HRL(I,KK) + WWVL3(K)*HRLH2O(I,KK,6)
        END DO
      END DO

C ---------------------------------------------------------------
C VI. H2O+N2O+O3 in 540-800cm-1 (12.5-18.52um) (CO2 included later)

c Main effect in the MA:(weak) cooling at ~x=3-7, zero above.
c H2O: no NLTE study
c N2O band 010-000. Source function: NLTE above ~75km (Gusev&Shved,JASTP,1997)
c O3 band 010-000 (14.1um). Source function: NLTE above ~75 km. 
c Heating rate for the 14.1um O3 band is studied by Manuilova et al. (2001).
c  O3 band contributes mainly below ~50 km (likely dominates H2O+N2O in this
c  spectral interval.
c Reduce heating rates to 0 between x=11.5-14 (similar to 1400-1900).

      DO K = 1, LAY
      KK = LAY+1-K
        DO I = IL1, IL2
        HRL(I,KK) = HRL(I,KK) + WWVL3(K)*HRLH2O(I,KK,7)
        END DO
      END DO

C ---------------------------------------------------------------
C VII. H2O in 2 INTERVALS: 340-540 (18.52-29.4um); 0-340 cm-1 (> 29.4um)

c Main effect in the MA: cooling, 0-340 dominates (especially above x=10)
c Rotational band. LTE (e.g., see Gordiets&Markov,1982).

      DO K = 1, LAY
        DO I = IL1, IL2
        HRL(I,K) = HRL(I,K) + HRLH2O(I,K,8)+HRLH2O(I,K,9)
        END DO
      END DO

C ---------------------------------------------------------------
C VIII. CO2 in 2 INTERVALS: 540-980 (10.2-18.52um)

c Main effect in the MA: dominates in the MA cooling
c parameterization by Fomichev et al. (1998) is used in the NLTE region
c above x=10.75. Between x=9.25-10.75 the NLTE heating rates are merged with
c values provided by the CKD scheme.

c *** Prepare input data for the NLTE parameterization ***

c * CO2-O deactivation rate constant (cm3/s):
      ZCO2O = 3.e-12

c * T interpolation to the parameterization grid (x=0-12.25)
      DO I = 1,50
        IF(xf(I).LE.xvf(LAY)) THEN
          K = IDR(I)
          KK= LAY+1-K
          DO IL = IL1,IL2
            TPCO2(IL,I) = PDR(1,I)*T(IL,KK)+PDR(2,I)*T(IL,KK-1)+
     &                                      PDR(3,I)*T(IL,KK-2)
          END DO
        ELSE
          DO IL = IL1,IL2
            TPCO2(IL,I) = T(IL,1)
          END DO
        END IF
      END DO

c * T, O2, O, N2, CO2 and AM interpolation to the param. grid (x=12.5-16.5)
      DO I = 51, 67
        IF(xf(I).LE.xvf(LAY)) THEN
          K = IDR(I)
          KK= LAY+1-K
          CO2R(I-50)=PDR(1,I)*CO2vf(K)+PDR(2,I)*CO2vf(K+1)+
     :                                  PDR(3,I)*CO2vf(K+2)
          O2R(I-50)=PDR(1,I)*O2vf(K)+PDR(2,I)*O2vf(K+1)+
     :                                  PDR(3,I)*O2vf(K+2)
          SN2R(I-50)=PDR(1,I)*CN2vf(K)+PDR(2,I)*CN2vf(K+1)+
     :                                  PDR(3,I)*CN2vf(K+2)
          OR(I-50)=PDR(1,I)*Ovf(K)+PDR(2,I)*Ovf(K+1)+
     :                                  PDR(3,I)*Ovf(K+2)
          AMR(I-50)=PDR(1,I)*AMuvf(K)+PDR(2,I)*AMuvf(K+1)+
     :                                  PDR(3,I)*AMuvf(K+2)
          DO IL = IL1,IL2
            TPCO2(IL,I) = PDR(1,I)*T(IL,KK)+PDR(2,I)*T(IL,KK-1)+
     :                                      PDR(3,I)*T(IL,KK-2)
          END DO
        ELSE
          CO2R(I-50)=CO2vf(LAY)
          O2R(I-50)=O2vf(LAY)
          SN2R(I-50)=CN2vf(LAY)
          OR(I-50)=Ovf(LAY)
          AMR(I-50)=AMuvf(LAY)
          DO IL = IL1,IL2
            TPCO2(IL,I) = T(IL,1)
          END DO
        END IF
      END DO

c * cooling rates at the parameterization grid
      CALL PCOOLG2(HPCO2,TPCO2,CO2R,O2R,SN2R,OR,AMR,ZCO2O,FLPCO2,
     1             SPCO2,ALPCO2,H2PCO2,ILG,IL1,IL2)

c * cooling rates at the model grid
      DO K = 1, LAY
         KK = LAY+1-K
        IF(xvf(K).LE.9.25) THEN
          DO IL = IL1,IL2
           HRL(IL,KK) = HRL(IL,KK) + HRLCO2(IL,KK,2)
          END DO
        ELSE IF(xvf(K).LE.10.75) THEN
          I = ICD(K)
          WCO2L = (10.75-xvf(K))/1.5
          DO IL = IL1,IL2
           CO2CLP = (PCD(1,K)*HPCO2(IL,I)+PCD(2,K)*HPCO2(IL,I+1)+
     &               PCD(3,K)*HPCO2(IL,I+2))/Cpvf(K)
           HRL(IL,KK) = HRL(IL,KK) +
     &                  WCO2L*HRLCO2(IL,KK,2)+(1.-WCO2L)*CO2CLP
          END DO
        ELSE IF(xvf(K).LE.16.5) THEN
          I = ICD(K)
          DO IL = IL1,IL2
           CO2CLP = (PCD(1,K)*HPCO2(IL,I)+PCD(2,K)*HPCO2(IL,I+1)+
     &               PCD(3,K)*HPCO2(IL,I+2))/Cpvf(K)
           HRL(IL,KK) = HRL(IL,KK) + CO2CLP
          END DO
        ELSE
c cooling rate above x=16.5
          DO IL = IL1,IL2
c quantum survival probability (alambda):
c --- V-T constants for O2 and N2
            tt=T(IL,KK)
            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(-xvf(K))/(tt*1.38066e-16)
c --- collisional deactivation rate:
            z=(CN2vf(K)*zn2+O2vf(K)*zo2+Ovf(K)*ZCO2O)*d
            alambda = a10/(a10+z)
c --- source function
            s = exp(aku/tt)
c --- cooling rate
            CO2CLP = const/AMuvf(K)*CO2vf(K)*
     &                            (1.-alambda)*(FLPCO2(IL)-s)/Cpvf(K)
            HRL(IL,KK) = HRL(IL,KK) + CO2CLP
          END DO
        END IF
      END DO

      RETURN
      END
