SUBROUTINE spa88(excts,ctso3,gxcts,sorc,csour,cldfac,temp,press,var1,var2,  &
     p,delp,delp2,totvo2,to3sp,to3spc,co2sp1,co2sp2,co2sp)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !subroutine spa88 computes exact cts heating rates and fluxes and
  !corresponding cts emissivity quantities for h2o,co2 and o3.
  IMPLICIT NONE

  INCLUDE "parmeta.f90"

  INTEGER,PARAMETER :: l=lm,imax=2*im-1,ncol=imax,nblw=163,nblx=47,nbly=15
  INTEGER,PARAMETER :: nblm=nbly-1,lp1=l+1,lp2=l+2,lp3=l+3,lm1=l-1,lm2=l-2
  INTEGER,PARAMETER :: lm3=l-3,ll=2*l,llp1=ll+1,llp2=ll+2,llp3=ll+3,llm1=ll-1
  INTEGER,PARAMETER :: llm2=ll-2,llm3=ll-3,lp1m=lp1*lp1,lp1m1=lp1m-1
  INTEGER,PARAMETER :: lp1v=lp1*(1+2*l/2),lp121=lp1*nbly,ll3p=3*l+2,nb=12
  INTEGER,PARAMETER :: inlte=3,inltep=inlte+1,nnlte=56,lp1i=imax*lp1
  INTEGER,PARAMETER :: llp1i=imax*llp1,ll3pi=imax*ll3p,nb1=nb-1,ko2=12
  INTEGER,PARAMETER :: ko21=ko2+1,ko2m=ko2-1

  REAL,INTENT(in),DIMENSION(imax,lp1,nbly) :: sorc
  REAL,INTENT(in),DIMENSION(imax,lp1,lp1) :: cldfac
  REAL,INTENT(in),DIMENSION(imax,lp1)     :: temp,to3sp,co2sp
  REAL,INTENT(in),DIMENSION(imax,l)       :: var1,var2,delp
  REAL,INTENT(in out),DIMENSION(imax,lp1) :: csour,press,co2sp1,totvo2,co2sp2,p
  REAL,INTENT(in out),DIMENSION(imax,l)   :: delp2,to3spc
  REAL,INTENT(out),DIMENSION(imax,l)      :: excts,ctso3
  REAL,INTENT(out),DIMENSION(imax)        :: gxcts

  INCLUDE "hcon.f90"
  INCLUDE "rnddta.f90" 

  REAL,DIMENSION(imax,l) :: phitmp,psitmp,tt,fac1,fac2,topm,topphi,f,ff,ag,agg
  REAL,DIMENSION(imax,l) :: x,y
  REAL,DIMENSION(imax,lp1) :: ctmp,ctmp3,ctmp2

  EQUIVALENCE (f,ag,phitmp)
  EQUIVALENCE (ff,agg,psitmp)

  INTEGER :: k,i
  !compute temperature quantities for use in program
  DO  k=1,l
     DO  i=1,imax
        x(i,k)=temp(i,k)-h25e2
        y(i,k)=x(i,k)*x(i,k)
     END DO
  END DO

  !initialize ctmp(i,1),ctmp2(i,1),ctmp3(i,1) to unity; these are
  !transmission fctns at the top.
  DO  i=1,imax
     ctmp(i,1)=one
     ctmp2(i,1)=1.
     ctmp3(i,1)=1.
  END DO

  !begin loop on frequency bands (1)***
  !calculation for band 1 (combined band 1)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(1)*x(i,1)+bpcm(1)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(1)*x(i,1)+btpcm(1)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(1)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(1)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*fac1(i,1)/SQRT(1.+fac2(i,1)))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=sorc(i,1,1)*(ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,1)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,1)-sorc(i,l,1)))
  END DO

  !calculation for band 2 (combined band 2)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(2)*x(i,1)+bpcm(2)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(2)*x(i,1)+btpcm(2)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(2)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(2)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*fac1(i,1)/SQRT(1.+fac2(i,1)))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,2)* (ctmp(i,2)-ctmp(i,1))
  END DO
  !---gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,2)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,2)-sorc(i,l,2)))
  END DO

  !calculation for band 3 (combined band 3)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(3)*x(i,1)+bpcm(3)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(3)*x(i,1)+btpcm(3)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(3)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(3)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*fac1(i,1)/SQRT(1.+fac2(i,1)))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,3)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,3)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,3)-sorc(i,l,3)))
  END DO

  !calculation for band 4 (combined band 4)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(4)*x(i,1)+bpcm(4)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(4)*x(i,1)+btpcm(4)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(4)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(4)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*fac1(i,1)/SQRT(1.+fac2(i,1)))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,4)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,4)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,4)-sorc(i,l,4)))
  END DO

  !calculation for band 5 (combined band 5)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(5)*x(i,1)+bpcm(5)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(5)*x(i,1)+btpcm(5)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(5)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(5)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(5)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,5)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,5)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,5)-sorc(i,l,5)))
  END DO

  !calculation for band 6 (combined band 6)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(6)*x(i,1)+bpcm(6)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(6)*x(i,1)+btpcm(6)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(6)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(6)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(6)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,6)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,6)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,6)-sorc(i,l,6)))
  END DO

  !calculation for band 7 (combined band 7)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(7)*x(i,1)+bpcm(7)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(7)*x(i,1)+btpcm(7)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(7)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(7)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(7)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,7)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,7)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,7)-sorc(i,l,7)))
  END DO

  !calculation for band 8 (combined band 8)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(8)*x(i,1)+bpcm(8)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(8)*x(i,1)+btpcm(8)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(8)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(8)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(8)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,8)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,8)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,8)-sorc(i,l,8)))
  END DO

  !calculation for band 9 ( 560-670 cm-1; includes co2)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(9)*x(i,1)+bpcm(9)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(9)*x(i,1)+btpcm(9)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(9)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(9)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(9)*totvo2(i,2)*sko2d))*co2sp1(i,2)
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,9)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,9)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,9)-sorc(i,l,9)))
  END DO

  !calculation for band 10 (670-800 cm-1; includes co2)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(10)*x(i,1)+bpcm(10)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(10)*x(i,1)+btpcm(10)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(10)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(10)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(10)*totvo2(i,2)*sko2d))*co2sp2(i,2)
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,10)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,10)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,10)-sorc(i,l,10)))
  END DO

  !calculation for band 11 (800-900 cm-1)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(11)*x(i,1)+bpcm(11)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(11)*x(i,1)+btpcm(11)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(11)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(11)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(11)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,11)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,11)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,11)-sorc(i,l,11)))
  END DO

  !calculation for band 12 (900-990 cm-1)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(12)*x(i,1)+bpcm(12)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(12)*x(i,1)+btpcm(12)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(12)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(12)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(12)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,12)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,12)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,12)-sorc(i,l,12)))
  END DO

  !calculation for band 13 (990-1070 cm-1; includes o3))
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(13)*x(i,1)+bpcm(13)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(13)*x(i,1)+btpcm(13)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(13)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(13)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(13)*totvo2(i,2)*sko2d +to3spc(i,1)))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,13)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,13)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,13)-sorc(i,l,13)))
  END DO

  !calculation for band 14 (1070-1200 cm-1)
  !obtain temperature correction (capphi,cappsi),then multiply
  !by optical path (var1,var2) to compute temperature-corrected
  !optical path and mean pressure for a layer (phitmp,psitmp)
  DO  i=1,imax*l
     f(i,1)=h44194m2*(apcm(14)*x(i,1)+bpcm(14)*y(i,1))
     ff(i,1)=h44194m2*(atpcm(14)*x(i,1)+btpcm(14)*y(i,1))
     ag(i,1)=(h1p41819+f(i,1))*f(i,1)+one
     agg(i,1)=(h1p41819+ff(i,1))*ff(i,1)+one
     phitmp(i,1)=var1(i,1)*(((( ag(i,1)*ag(i,1))**2)**2)**2)
     psitmp(i,1)=var2(i,1)*(((( agg(i,1)*agg(i,1))**2)**2)**2)
  END DO

  !obtain optical path,mean pressure from the top to the pressure
  !p(k) (topm,topphi)
  DO  i=1,imax
     topm(i,1)=phitmp(i,1)
     topphi(i,1)=psitmp(i,1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        topm(i,k)=topm(i,k-1)+phitmp(i,k)
        topphi(i,k)=topphi(i,k-1)+psitmp(i,k)
     END DO
  END DO

  !tt is the cloud-free cts transmission function
  DO  i=1,imax*l
     fac1(i,1)=acomb(14)*topm(i,1)
     fac2(i,1)=fac1(i,1)*topm(i,1)/(bcomb(14)*topphi(i,1))
     tt(i,1)=EXP(hm1ez*(fac1(i,1)/SQRT(one+fac2(i,1))+  &
          betacm(14)*totvo2(i,2)*sko2d))
     ctmp(i,2)=tt(i,1)*cldfac(i,2,1)
  END DO

  !excts is the cts cooling rate accumulated over frequency bands
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)+sorc(i,1,14)* (ctmp(i,2)-ctmp(i,1))
  END DO

  !gxcts is the exact cts top flux accumulated over frequency bands
  DO  i=1,imax
     gxcts(i)=gxcts(i)+cldfac(i,lp1,1)*(tt(i,l)*sorc(i,l,14)+  &
          (haf*delp(i,l)*(tt(i,lm1)*(p(i,lp1)-press(i,l)) +  &
          tt(i,l)*(p(i,lp1)+press(i,l)-two*p(i,l)))) * (sorc(i,lp1,14)-sorc(i,l,14)))
  END DO

  !obtain cts flux at the top by integration of heating rates and
  !using cts flux at the bottom (current value of gxcts). note
  !that the pressure quantities and conversion factors have not
  !been included either in excts or in gxcts. these cancel out, thus
  !reducing computations!
  DO  k=1,l
     DO  i=1,imax
        gxcts(i)=gxcts(i)-excts(i,k)
     END DO
  END DO

  !now scale the cooling rate (excts) by including the pressure
  !factor (delp) and the conversion factor (radcon)
  DO  i=1,imax*l
     excts(i,1)=excts(i,1)*radcon*delp(i,1)
  END DO
  !this is the end of the exact cts computations; at this point
  !excts has its appropriate value.

  !compute approximate cts heating rates for 15um and 9.6 um bands
  !(ctso3)
  DO  i=1,imax*l
     ctmp2(i,2)=co2sp(i,2)*cldfac(i,2,1)
     ctmp3(i,2)=to3sp(i,1)*cldfac(i,2,1)
  END DO
  DO  i=1,imax*l
     ctso3(i,1)=radcon*delp(i,1)* (csour(i,1)*(ctmp2(i,2)-ctmp2(i,1)) +  &
          sorc(i,1,13)*(ctmp3(i,2)-ctmp3(i,1)))
  END DO
END SUBROUTINE spa88
