#include "dims.h"
      subroutine vicool
      implicit none
! 
#include "params.h"
#include "radcool.h"
C
C XL(nzpm9_zpm5)  - the parameters for NLTE region (12.5 <= X <= 16.5)
      real :: amat,bmat,al
      COMMON /CO2CFG/ AMAT(43,9),BMAT(43,9)
      common /CO2CFG_PVT/ AL(zimxp,nzpm9_zpm5)
!$OMP THREADPRIVATE (/co2cfg_pvt/)
!DIR$ TASKCOMMON CO2CFG_PVT
      integer :: ig
      COMMON /PIRMGR/ IG(9)
      real :: ao3
      COMMON /PIRO3/ AO3(35,9)
!
! Local: 
! (RKO moved to radcool.h)
      real :: FU(zimxp,nzpsrf_zpm5),FO3(zimxp,nzpsrf_zpm5),H1(zimxp),
     |  H2(zimxp),H3(zimxp),TT(zimxp),Y(zimxp),ZN2(zimxp),ZO2(zimxp),
     |  ZZ(zimxp),AA1(zimxp),AA2(zimxp),D1(zimxp),D2(zimxp),
     |  ALAM(zimxp,nzpm9_zpm5)
      real :: CONST=2.63187E11  
      real :: CONSTB=9.08795E9
      real :: BOLTZ=1.38066E-16 
      real :: A10=1.5988
      integer :: len1,i,k,jj,jjs,ks,km
C
C grid levels for height integration   (p.s.h. distance = 0.25*IG)
C FU(nzpsrf_zpm5) - the LTE source functions for 15 um CO2 band at 
C   p.s.h.=0-16.5,
C FO3(45)- the same for 9.6 um O3 band at p.s.h. = 0-11, 0.25
C AL(nzpm9_zpm5) - the quantum survival probability for 
C   p.s.h. = 12.5-16.5, 0.25
C                                          (for 15 um CO2 band only)
      LEN1=zimxp
      DO K=1,nzpsrf_zpm5
        DO I=1,LEN1
          FU(I,K)=EXP(-960.217/TV(I,K))
          FO3(I,K)=EXP(-1500./TV(I,K))
        enddo
      enddo
C
C calculate the heating rates for layer below s.h.p. = 12.5
C   15 um CO2 + 9.6 um O3:
C
C     **** COOLING RATE IN BOTH O3 AND CO2 BANDS (X=2-10.5) MATRIX
C     **** APPROACH
C
      DO K=1,5
        KS = K+8
        DO I=1,LEN1
          H2(I) = (AMAT(K,1)+BMAT(K,1)*FU(I,KS))*FU(I,1)
          H3(I) = AO3(K,1)*FO3(I,1)
        enddo
        DO JJ=3,9
          JJS = KS+IG(JJ)
          DO I=1,LEN1
            H2(I) = H2(I)+(AMAT(K,JJ)+BMAT(K,JJ)*FU(I,KS))*FU(I,JJS)
            H3(I) = H3(I)+AO3(K,JJ)*FO3(I,JJS)
          enddo
        enddo
        DO I=1,LEN1
          HCO2(I,K) = H2(I)
          HO3(I,K) = H3(I)*O3(I,KS)
        enddo
      enddo
C
      DO K=6,18
        KS = K+8
        DO I=1,LEN1
          H2(I) = (AMAT(K,1)+BMAT(K,1)*FU(I,KS))*FU(I,1)
          H3(I) = AO3(K,1)*FO3(I,1)
        enddo
        DO JJ=2,9
          JJS = KS+IG(JJ)
          DO I=1,LEN1
            H2(I) = H2(I)+(AMAT(K,JJ)+BMAT(K,JJ)*FU(I,KS))*FU(I,JJS)
            H3(I) = H3(I)+AO3(K,JJ)*FO3(I,JJS)
          enddo
        enddo
        DO I=1,LEN1
          HCO2(I,K) = H2(I)
          HO3(I,K) = H3(I)*O3(I,KS)
        enddo
      enddo
C
      DO K=19,35
        KS = K+8
        DO I=1,LEN1
          H2(I) = 0.
          H3(I) = 0.
        enddo
        DO JJ=1,9
          JJS = KS+IG(JJ)
          DO I=1,LEN1
            H2(I) = H2(I)+(AMAT(K,JJ)+BMAT(K,JJ)*FU(I,KS))*FU(I,JJS)
            H3(I) = H3(I)+AO3(K,JJ)*FO3(I,JJS)
          enddo
        enddo
        DO I=1,LEN1
          HCO2(I,K) = H2(I)
          HO3(I,K) = H3(I)*O3(I,KS)
        enddo
      enddo
C
C     **** COOLING RATE IN CO2 BANDS (X=10.75-12.5, MATRIX APPROACH)
C
      DO K=36,43
        KS = K+8
        DO I=1,LEN1
          H2(I) = 0.
        enddo
        DO JJ=1,9
          JJS = KS+IG(JJ)
          DO I=1,LEN1
            H2(I) = H2(I)+(AMAT(K,JJ)+BMAT(K,JJ)*FU(I,KS))*FU(I,JJS)
          enddo
        enddo
        DO I=1,LEN1
          HCO2(I,K) = H2(I)
          HO3(I,K) = 0.
        enddo
      enddo
C
C     ****  DETERMINE LAMBDA(17)
C
      DO K=1,nzpm9_zpm5 ! see radcool.h
        DO I=1,LEN1
C
C     ****  CO2-O2 AND CO2-N2 V-T CONSTANTS
C
          TT(I) = TV(I,K+50)
          Y(I) = TT(I)**(-1./3.)
          ZN2(I) = 5.5E-17*SQRT(TT(I))+6.7E-10*EXP(-83.8*Y(I))
          ZO2(I) = 1.E-15*EXP(23.37-230.9*Y(I)+564.*Y(I)*Y(I))
C
C     ****  COLLISIONAL DEACTIVATION RATE:
C
          ZZ(I) = (SN2(I,K)*ZN2(I)+SO2(I,K)*ZO2(I)+O(I,K)*RKO)*DEN(I,K)
          ALAM(I,K) = A10/(A10+ZZ(I))
        enddo
      enddo
C
      DO I=1,LEN1
        H1(I)=HCO2(I,43)/(CO2(I,1)*(1.-ALAM(I,1))*CONSTB)
      enddo
      DO K=2,nzpm9_zpm5
        KM=K-1
        DO I=1,LEN1
          AA1(I)=1.-ALAM(I,KM)*(1.-.25*AL(I,K)-.75*AL(I,KM))
          AA2(I)=1.-ALAM(I,K)*(1.-.75*AL(I,K)-.25*AL(I,KM))
          D1(I) = -.25*(AL(I,K)+3.*AL(I,KM))
          D2(I) = .25*(3.*AL(I,K)+AL(I,KM))
          H2(I)=(AA1(I)*H1(I)-D1(I)*FU(I,KM+50)-D2(I)*FU(I,K+50))/AA2(I)
          HCO2(I,K+42)=H2(I)*CO2(I,K)*(1.-ALAM(I,K))/AM(I,K)*CONST
          HO3(I,K+42)=0.
          H1(I)=H2(I)
        enddo
      enddo
      DO I=1,LEN1
        FLUX(I) = H2(I) + FU(I,nzpsrf_zpm5)
      enddo
      RETURN
      END