SUBROUTINE fst88(heatra,grnflx,topflx, qh2o,press,p,delp,delp2,temp,t,  &
     cldfac,nclds,ktop,kbtm,camt, co21,co2nbl,co2sp1,co2sp2,  &
     var1,var2,var3,var4,cntval, toto3,tphio3,totphi,totvo2,  &
     emx1,emx2,empl)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !subroutine fst88 is the main computation module of the
  !long-wave radiation code. in it all "emissivity" calculations,
  !including calls to table lookup subroutines. also,after calling
  !subroutine "spa88", final combined heating rates and ground
  !flux are obtained.
  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

  INTEGER :: i,k,kp,klen,kk,icnt,kclds,kmin,kmax,j1,j3

  REAL,INTENT(out),DIMENSION(imax,l)          :: heatra
  REAL,INTENT(out),DIMENSION(imax)            :: grnflx,topflx
  REAL,INTENT(in out),DIMENSION(imax)         :: emx2
  REAL,INTENT(in out),DIMENSION(imax,l)       :: var1
  REAL,INTENT(in out),DIMENSION(imax,lp1,lp1) :: co21,cldfac
  REAL,INTENT(in out),DIMENSION(imax,lp1)     :: qh2o,camt,co2sp1,co2sp2
  REAL,INTENT(in),DIMENSION(imax,llp1)        :: empl
  REAL,INTENT(in),DIMENSION(imax,l)           :: delp,delp2,co2nbl,var2,var3
  REAL,INTENT(in),DIMENSION(imax,l)           :: var4
  REAL,INTENT(in),DIMENSION(imax,lp1)         :: press,p,temp,t,cntval
  REAL,INTENT(in),DIMENSION(imax,lp1)         :: toto3,tphio3,totphi,totvo2
  REAL,INTENT(in),DIMENSION(imax)             :: emx1

  INTEGER,INTENT(in),DIMENSION(imax)         :: nclds
  INTEGER,INTENT(in),DIMENSION(imax,lp1)    :: ktop(imax,lp1),kbtm(imax,lp1)

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


  REAL,DIMENSION(imax) :: gxcts,flx1e1,delptc,ptop,pbot,ftop,fbot 
  REAL,DIMENSION(imax,2) :: emspec,fxosp,dtsp
  REAL,DIMENSION(imax,l) :: excts,ctso3,cts,e1cts2,e1ctw2,rlog,to3spc
  REAL,DIMENSION(imax,lp1) :: avephi,emiss,emissb,co2sp,to3sp,oss,css,ss1,ss2
  REAL,DIMENSION(imax,lp1) :: tc,dtc,csour,avvo2,heatem,over1d,to31d,cont1d
  REAL,DIMENSION(imax,lp1) :: avmo3,avpho3,vtmp3,dsorc,fac1,delpr1,delpr2
  REAL,DIMENSION(imax,lp1) :: emisdg,contdg,to3dg,flxnet,ixo,vsum1,flxthk,z1
  REAL,DIMENSION(imax,lp1) :: e1cts1,e1ctw1,fxo,dt,fxoe2,dte2,flx,totevv,cnttau
  REAL,DIMENSION(imax,lp1) :: e1flx
  REAL,DIMENSION(imax,llp1) :: c,c2,alp,csub,csub2,emd,tpl
  REAL,DIMENSION(imax,lp1,nbly) :: sorc


  EQUIVALENCE (alp,c,csub),(csub2,c2)
  EQUIVALENCE (fac1,dsorc,over1d,delpr2,flxnet)
  EQUIVALENCE (delpr1,heatem)
  EQUIVALENCE (ixo,avvo2,flxthk,to3dg)
  EQUIVALENCE (z1,avmo3,contdg)
  EQUIVALENCE (emisdg,vsum1,avpho3)
  EQUIVALENCE (emd(1,1),e1cts1(1,1)),(emd(1,lp2),e1cts2(1,1))
  EQUIVALENCE (tpl(1,1),e1ctw1(1,1)),(tpl(1,lp2),e1ctw2(1,1))

  !first section is table lookup for source function and
  !derivative (b and db/dt).also,the nlte co2 source function
  !is obtained in calcs. below, decrementing the index by 9
  !accounts for the tables beginning at t=100k.
  !at t=100k.

  DO  i=1,imax*lp1
     !temp. indices for e1,source
     vtmp3(i,1)=AINT(temp(i,1)*hp1)
     fxo(i,1)=vtmp3(i,1)-9.
     dt(i,1)=temp(i,1)-ten*vtmp3(i,1)
     !integer index for source (used immediately)
     ixo(i,1)=fxo(i,1)
  END DO
  DO  i=1,imax*l
     !temp. indices for e2 (kp=1 layer not used in flux calculations)
     vtmp3(i,1)=AINT(t(i,2)*hp1)
     fxoe2(i,1)=vtmp3(i,1)-9.
     dte2(i,1)=t(i,2)-ten*vtmp3(i,1)
  END DO
  !special case to handle kp=lp1 layer and special e2 calcs.
  DO  i=1,imax
     fxoe2(i,lp1)=fxo(i,l)
     dte2(i,lp1)=dt(i,l)
     fxosp(i,1)=fxoe2(i,lm1)
     fxosp(i,2)=fxo(i,lm1)
     dtsp(i,1)=dte2(i,lm1)
     dtsp(i,2)=dt(i,lm1)
  END DO

  !source function for combined band 1

  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),1)
        dsorc(i,k)=dsrce(ixo(i,k),1)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,1)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 2
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),2)
        dsorc(i,k)=dsrce(ixo(i,k),2)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,2)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 3
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),3)
        dsorc(i,k)=dsrce(ixo(i,k),3)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,3)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 4
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),4)
        dsorc(i,k)=dsrce(ixo(i,k),4)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,4)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 5
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),5)
        dsorc(i,k)=dsrce(ixo(i,k),5)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,5)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 6
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),6)
        dsorc(i,k)=dsrce(ixo(i,k),6)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,6)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 7
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),7)
        dsorc(i,k)=dsrce(ixo(i,k),7)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,7)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for combined band 8
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),8)
        dsorc(i,k)=dsrce(ixo(i,k),8)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,8)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 9 (560-670 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),9)
        dsorc(i,k)=dsrce(ixo(i,k),9)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,9)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 10 (670-800 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),10)
        dsorc(i,k)=dsrce(ixo(i,k),10)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,10)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 11 (800-900 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),11)
        dsorc(i,k)=dsrce(ixo(i,k),11)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,11)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 12 (900-990 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),12)
        dsorc(i,k)=dsrce(ixo(i,k),12)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,12)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 13 (990-1070 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),13)
        dsorc(i,k)=dsrce(ixo(i,k),13)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,13)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO
  !source function for band 14 (1070-1200 cm-1)
  DO  i=1,imax
     DO  k=1,lp1
        vtmp3(i,k)=source(ixo(i,k),14)
        dsorc(i,k)=dsrce(ixo(i,k),14)
     END DO
  END DO
  DO  i=1,imax*lp1
     sorc(i,1,14)=vtmp3(i,1)+dt(i,1)*dsorc(i,1)
  END DO

  !obtain special source functions for the 15 um band (csour)
  !and the window region (ss1)
  DO  i=1,imax*lp1
     ss1(i,1)=sorc(i,1,11)+sorc(i,1,12)+sorc(i,1,14)
  END DO
  DO  i=1,imax*lp1
     csour(i,1)=sorc(i,1,9)+sorc(i,1,10)
  END DO

  !compute temp**4 (tc) and vertical temperature differences
  !(oss,css,ss2,dtc). all these will be used later in flux computations.
  DO  i=1,imax*lp1
     tc(i,1)=(temp(i,1)*temp(i,1))**2
  END DO
  DO  i=1,imax*l
     oss(i,2)=sorc(i,2,13)-sorc(i,1,13)
     css(i,2)=csour(i,2)-csour(i,1)
     dtc(i,2)=tc(i,2)-tc(i,1)
     ss2(i,2)=ss1(i,2)-ss1(i,1)
  END DO

  !the followimg is a drastic rewrite of the radiation code to
  !(largely) eliminate three-dimensional arrays. the code works
  !on the following principles:
  !  let k = fixed flux level, kp = varying flux level
  !  then flux(k)=sum over kp : (deltab(kp)*tau(kp,k))
  !       over all kp's, from 1 to lp1.
  !          we can break down the calculations for all k's as follows:
  !          for all k's k=1 to lp1:
  !              flux(k)=sum over kp : (deltab(kp)*tau(kp,k))  (1)
  !                      over all kp's, from k+1 to lp1
  !          and
  !              for kp from k+1 to lp1:
  !                 flux(kp) = deltab(k)*tau(k,kp)              (2)
  !          now if tau(k,kp)=tau(kp,k) (symmetrical arrays)
  !          we can compute a 1-dimensional array tau1d(kp) from
  !          k+1 to lp1, each time k is incremented.
  !          equations (1) and (2) then become:
  !             tau1d(kp) = (values for tau(kp,k) at the particular k)
  !             flux(k) = sum over kp : (deltab(kp)*tau1d(kp))   (3)
  !             flux(kp) = deltab(k)*tau1d(kp)                   (4)
  !         the terms for tau (k,k) and other special terms (for
  !         nearby layers) must, of course, be handled separately, and
  !         with care.
  !      compute "upper triangle" transmission functions for
  !      the 9.6 um band (to3sp) and the 15 um band (over1d). also,
  !      the
  !      stage 1...compute o3 ,over transmission fctns and avephi
  !---do k=1 calculation (from flux layer kk to the top) separately
  !   as vectorization is improved,and ozone cts transmissivity
  !   may be extracted here.
  DO  i=1,imax*l
     avephi(i,1)=totphi(i,2)
  END DO
  !in order to properly evaluate emiss integrated over the (lp1)
  !layer, a special evaluation of emiss is done. this requires
  !a special computation of avephi, and it is stored in the
  !(otherwise vacant) lp1'th position
  DO  i=1,imax
     avephi(i,lp1)=avephi(i,lm1)+emx1(i)
  END DO
  !compute fluxes for k=1
  CALL e1e290(e1cts1,e1cts2,e1flx,e1ctw1,e1ctw2,emiss,  &
       fxo,dt,fxoe2,dte2,avephi,temp,t)
  DO  i=1,imax*l
     fac1(i,1)=bo3rnd(2)*tphio3(i,2)/toto3(i,2)
     to3spc(i,1)=haf*(fac1(i,1)*  &
          (SQRT(one+(four*ao3rnd(2)*toto3(i,2))/fac1(i,1))-one))
     !for k=1, to3sp is used instead of to31d (they are equal in this
     !case); to3sp is passed to spa90, while to31d is a work-array.
     to3sp(i,1)=EXP(hm1ez*(to3spc(i,1)+sko3r*totvo2(i,2)))
     over1d(i,1)=EXP(hm1ez*(SQRT(ab15wd*totphi(i,2))+ skc1r*totvo2(i,2)))
     !because all continuum transmissivities are obtained from the
     !2-d quantity cnttau (and its reciprocal totevv) we store both
     !of these here. for k=1, cont1d equals cnttau
     cnttau(i,1)=EXP(hm1ez*totvo2(i,2))
     totevv(i,1)=1./cnttau(i,1)
  END DO
  DO  k=1,l
     DO  i=1,imax
        co2sp(i,k+1)=over1d(i,k)*co21(i,1,k+1)
     END DO
  END DO
  DO  i=1,imax*l
     co21(i,2,1)=co21(i,2,1)*over1d(i,1)
  END DO
  !rlog is the nbl amount for the 15 um band calculation
  DO  i=1,imax
     rlog(i,1)=over1d(i,1)*co2nbl(i,1)
  END DO
  !the terms when kp=1 for all k are the photon exchange with
  !the top of the atmosphere, and are obtained differently than
  !the other calculations
  DO  k=2,lp1
     DO  i=1,imax
        flx(i,k)= (tc(i,1)*e1flx(i,k) +ss1(i,1)*cnttau(i,k-1)  &
             +sorc(i,1,13)*to3sp(i,k-1) +csour(i,1)*co2sp(i,k))  &
             *cldfac(i,1,k)
     END DO
  END DO
  DO  i=1,imax
     flx(i,1)= tc(i,1)*e1flx(i,1)+ss1(i,1)+sorc(i,1,13) +csour(i,1)
  END DO
  !the kp terms for k=1...
  DO  kp=2,lp1
     DO  i=1,imax
        flx(i,1)=flx(i,1)+(oss(i,kp)*to3sp(i,kp-1) +ss2(i,kp)*cnttau(i,kp-1)  &
             +css(i,kp)*co21(i,kp,1) +dtc(i,kp)*emiss(i,kp-1))*cldfac(i,kp,1)
     END DO
  END DO

  !subroutine spa88 is called to obtain exact cts for water
  !co2 and o3, and approximate cts co2 and o3 calculations.
  CALL spa88(excts,ctso3,gxcts,sorc,csour, cldfac,temp,press,var1,var2,  &
       p,delp,delp2,totvo2,to3sp,to3spc, co2sp1,co2sp2,co2sp)

  !this section computes the emissivity cts heating rates for 2
  !emissivity bands: the 0-160,1200-2200 cm-1 band and the 800-
  !990,1070-1200 cm-1 band. the remaining cts comtributions are
  !contained in ctso3, computed in spa88.
  DO  i=1,imax
     vtmp3(i,1)=1.
  END DO
  DO  i=1,imax*l
     vtmp3(i,2)=cnttau(i,1)*cldfac(i,2,1)
  END DO
  DO  i=1,imax*l
     cts(i,1)=radcon*delp(i,1)*(tc(i,1)*  &
          (e1ctw2(i,1)*cldfac(i,2,1)-e1ctw1(i,1)*cldfac(i,1,1)) +  &
          ss1(i,1)*(vtmp3(i,2)-vtmp3(i,1)))
  END DO
  DO  i=1,imax*l
     vtmp3(i,1)=tc(i,1)*(cldfac(i,1,1)*(e1cts1(i,1)-e1ctw1(i,1)) -  &
          cldfac(i,2,1)*(e1cts2(i,1)-e1ctw2(i,1)))
  END DO
  DO  i=1,imax
     flx1e1(i)=tc(i,lp1)*cldfac(i,lp1,1)* (e1cts1(i,lp1)-e1ctw1(i,lp1))
  END DO
  DO  k=1,l
     DO  i=1,imax
        flx1e1(i)=flx1e1(i)+vtmp3(i,k)
     END DO
  END DO

  !now repeat flux calculations for the k=2..lm1  cases.
  !calculations for flux level l and lp1 are done separately, as all
  !emissivity and co2 calculations are special cases or nearby layers.
  DO  k=2,lm1
     klen=k
     DO  kk=1,lp1-k
        DO  i=1,imax
           avephi(i,kk+k-1)=totphi(i,kk+k)-totphi(i,k)
        END DO
     END DO
     DO  i=1,imax
        avephi(i,lp1)=avephi(i,lm1)+emx1(i)
     END DO
     !compute emissivity fluxes (e2) for this case. note that
     !we have omitted the nearby later case (emiss(i,k,k)) as well
     !as all cases with k=l or lp1. but these cases have always
     !been handled as special cases, so we may as well compute
     !their fluxes separastely.
     CALL e290(emissb,emiss,avephi,klen,fxoe2,dte2)
     DO  kk=1,lp1-k
        DO  i=1,imax
           avmo3(i,kk+k-1)=toto3(i,kk+k)-toto3(i,k)
           avpho3(i,kk+k-1)=tphio3(i,kk+k)-tphio3(i,k)
           avvo2(i,kk+k-1)=totvo2(i,kk+k)-totvo2(i,k)
           cont1d(i,kk+k-1)=cnttau(i,kk+k-1)*totevv(i,k-1)
        END DO
     END DO
     DO  i=1,imax*(lp1-k)
        fac1(i,k)=bo3rnd(2)*avpho3(i,k)/avmo3(i,k)
        vtmp3(i,k)=haf*(fac1(i,k)*  &
             (SQRT(one+(four*ao3rnd(2)*avmo3(i,k))/fac1(i,k))-one))
        to31d(i,k)=EXP(hm1ez*(vtmp3(i,k)+sko3r*avvo2(i,k)))
        over1d(i,k)=EXP(hm1ez*(SQRT(ab15wd*avephi(i,k))+ skc1r*avvo2(i,k)))
        co21(i,k+1,k)=over1d(i,k)*co21(i,k+1,k)
     END DO
     DO  kp=k+1,lp1
        DO  i=1,imax
           co21(i,k,kp)=over1d(i,kp-1)*co21(i,k,kp)
        END DO
     END DO
     !rlog is the nbl amount for the 15 um band calculation
     DO  i=1,imax
        rlog(i,k)=over1d(i,k)*co2nbl(i,k)
     END DO
     !the kp terms for arbirrary k..
     DO  kp=k+1,lp1
        DO  i=1,imax
           flx(i,k)=flx(i,k)+(oss(i,kp)*to31d(i,kp-1) +ss2(i,kp)*cont1d(i,kp-1)  &
                +css(i,kp)*co21(i,kp,k) +dtc(i,kp)*emiss(i,kp-1))*cldfac(i,kp,k)
        END DO
     END DO
     DO  kp=k+1,lp1
        DO  i=1,imax
           flx(i,kp)=flx(i,kp)+(oss(i,k)*to31d(i,kp-1) + &
                ss2(i,k)*cont1d(i,kp-1) + &
                css(i,k)*co21(i,k,kp) + &
                dtc(i,k)*emissb(i,kp-1))*cldfac(i,k,kp)
        END DO
     END DO
  END DO
321 CONTINUE

  !now do k=l case. since the kp loop is length 1, many simplifi-
  !cations occur. also, the co2 quantities (as well as the emiss
  !quantities) are computed in the nbl sedction; therefore, we want
  !only over,to3 and cont1d (over(i,l),to31d(i,l) and cont1d(i,l)
  !according to the notation. thus no call is made to the e290
  !subroutine.
  !         the third section calculates boundary layer and nearby layer
  !     corrections to the transmission functions obtained above. methods
  !     are given in ref. (4).
  !          the following ratios are used in various nbl calculations:
  !   the remaining calculations are for :
  !                        1) the (k,k) terms, k=2,lm1;
  !                        2) the (l,l) term
  !                        3) the (l,lp1) term
  !                        4) the (lp1,l) term
  !                        5) the (lp1,lp1) term.
  !     each is uniquely handled; different flux terms are computed
  !     differently
  !          fourth section obtains water transmission functions
  !     used in q(approx) calculations and also makes nbl corrections:
  !     1) emiss (i,j) is the transmission function matrix obtained
  !     by calling subroutine e1e288;
  !     2) "nearby layer" corrections (emiss(i,i)) are obtained
  !     using subroutine e3v88;
  !     3) special values at the surface (emiss(l,lp1),emiss(lp1,l),
  !     emiss(lp1,lp1)) are calculated.
  !obtain arguments for e1e288 and e3v88:
  DO  i=1,imax
     tpl(i,1)=temp(i,l)
     tpl(i,lp1)=haf*(t(i,lp1)+temp(i,l))
     tpl(i,llp1)=haf*(t(i,l)+temp(i,l))
  END DO
  DO  k=2,l
     DO  i=1,imax
        tpl(i,k)=t(i,k)
        tpl(i,k+l)=t(i,k)
     END DO
  END DO

  !e2 functions are required in the nbl calculations for 2 cases,
  !denoted (in old code) as (l,lp1) and (lp1,lp1)
  DO  i=1,imax
     avephi(i,1)=var2(i,l)
     avephi(i,2)=var2(i,l)+empl(i,l)
  END DO
  CALL e2spec(emiss,avephi,fxosp,dtsp)

  !call e3v88 for nbl h2o transmissivities
  CALL e3v88(emd,tpl,empl)

  !compute nearby layer and special-case transmissivities for emiss
  !using methods for h2o given in ref. (4)
  DO  k=2,l
     DO  i=1,imax
        emisdg(i,k)=emd(i,k+l)+emd(i,k)
     END DO
  END DO

  !note that emx1/2 (pressure scaled paths) are now computed in lwr88
  DO  i=1,imax
     emspec(i,1)=(emd(i,1)*empl(i,1)-emd(i,lp1)*empl(i,lp1))/  &
          emx1(i) + quartr*(emiss(i,1)+emiss(i,2))
     emisdg(i,lp1)=two*emd(i,lp1)
     emspec(i,2)=two*(emd(i,1)*empl(i,1)-emd(i,llp1)*empl(i,llp1))/ emx2(i)
  END DO
  DO  i=1,imax
     fac1(i,l)=bo3rnd(2)*var4(i,l)/var3(i,l)
     vtmp3(i,l)=haf*(fac1(i,l)*  &
          (SQRT(one+(four*ao3rnd(2)*var3(i,l))/fac1(i,l))-one))
     to31d(i,l)=EXP(hm1ez*(vtmp3(i,l)+sko3r*cntval(i,l)))
     over1d(i,l)=EXP(hm1ez*(SQRT(ab15wd*var2(i,l))+ skc1r*cntval(i,l)))
     cont1d(i,l)=cnttau(i,l)*totevv(i,lm1)
     rlog(i,l)=over1d(i,l)*co2nbl(i,l)
  END DO
  DO  i=1,imax*l
     rlog(i,1)=LOG(rlog(i,1))
  END DO
  DO  i=1,imax*lm1
     delpr1(i,2)=delp(i,2)*(press(i,2)-p(i,2))
     alp(i,lp1)=-SQRT(delpr1(i,2))*rlog(i,2)
  END DO
  DO  i=1,imax*l
     delpr2(i,2)=delp(i,1)*(p(i,2)-press(i,1))
     alp(i,1)=-SQRT(delpr2(i,2))*rlog(i,1)
  END DO
  DO  i=1,imax
     alp(i,ll)=-rlog(i,l)
     alp(i,llp1)=-rlog(i,l)*SQRT(delp(i,l)*(p(i,lp1)-press(i,lm1)))
  END DO
  !the first computation is for the 15 um band,with the
  !for the combined h2o and co2 transmission function.
  !perform nbl computations for the 15 um band
  !the statement function sf in prev. versions is now explicitly
  !evaluated.
  DO  i=1,imax*llp1
     c(i,1)=alp(i,1)*(hmp66667+alp(i,1)*(quartr+alp(i,1)*hm6666m2))
  END DO
  DO  i=1,imax
     co21(i,lp1,lp1)=one+c(i,l)
     co21(i,lp1,l)=one+(delp2(i,l)*c(i,ll)-(press(i,l)-p(i,l))*  &
          c(i,llm1))/(p(i,lp1)-press(i,l))
     co21(i,l,lp1)=one+((p(i,lp1)-press(i,lm1))*c(i,llp1)-  &
          (p(i,lp1)-press(i,l))*c(i,l))/(press(i,l)-press(i,lm1))
  END DO
  DO  k=2,l
     DO  i=1,imax
        co21(i,k,k)=one+haf*(c(i,lm1+k)+c(i,k-1))
     END DO
  END DO

  !compute nearby-layer transmissivities for the o3 band and for the
  !one-band continuum band (to3 and emiss2). the sf2 function is
  !used. the method is the same as described for co2 in ref (4).
  DO  i=1,imax*lm1
     csub(i,2)=cntval(i,2)*delpr1(i,2)
     csub(i,lp1)=cntval(i,1)*delpr2(i,2)
  END DO
  !the sf2 function in prev. versions is now explicitly evaluated
  DO  i=1,imax*llm2
     csub2(i,2)=sko3r*csub(i,2)
     c(i,2)=csub(i,2)*(hmp5+csub(i,2)*(hp166666-csub(i,2)*h41666m2))
     c2(i,2)=csub2(i,2)*(hmp5+csub2(i,2)* (hp166666-csub2(i,2)*h41666m2))
  END DO
  DO  i=1,imax
     contdg(i,lp1)=1.+c(i,llm1)
     to3dg(i,lp1)=1.+c2(i,llm1)
  END DO
  DO  k=2,l
     DO  i=1,imax
        contdg(i,k)=one+haf*(c(i,k)+c(i,lm1+k))
        to3dg(i,k)=one+haf*(c2(i,k)+c2(i,lm1+k))
     END DO
  END DO
  !now obtain fluxes

  !for the diagonal terms...
  DO  k=2,lp1
     DO  i=1,imax
        flx(i,k)=flx(i,k)+(dtc(i,k)*emisdg(i,k) +ss2(i,k)*contdg(i,k)  &
             +oss(i,k)*to3dg(i,k) +css(i,k)*co21(i,k,k))*cldfac(i,k,k)
     END DO
  END DO
  !for the two off-diagonal terms...
  DO  i=1,imax
     flx(i,l)=flx(i,l)+(css(i,lp1)*co21(i,lp1,l) +dtc(i,lp1)*emspec(i,2)  &
          +oss(i,lp1)*to31d(i,l) +ss2(i,lp1)*cont1d(i,l))*cldfac(i,lp1,l)
     flx(i,lp1)=flx(i,lp1)+(css(i,l)*co21(i,l,lp1) +oss(i,l)*to31d(i,l)  &
          +ss2(i,l)*cont1d(i,l) +dtc(i,l)*emspec(i,1))*cldfac(i,l,lp1)
  END DO

  !final section obtains emissivity heating rates,
  !total heating rates and the flux at the ground
  !calculate the emissivity heating rates
  DO  i=1,imax*l
     heatem(i,1)=radcon*(flx(i,2)-flx(i,1))*delp(i,1)
  END DO
  !calculate the total heating rates
  DO  i=1,imax*l
     heatra(i,1)=heatem(i,1)-cts(i,1)-ctso3(i,1)+excts(i,1)
  END DO
  !calculate the flux at each flux level using the flux at the
  !top (flx1e1+gxcts) and the integral of the heating rates (vsum1)
  DO  i=1,imax*l
     vsum1(i,1)=heatra(i,1)*delp2(i,1)*radcon1
  END DO
  DO  i=1,imax
     topflx(i)=flx1e1(i)+gxcts(i)
     flxnet(i,1)=topflx(i)
  END DO
  !only the surface value of flux (grnflx) is needed unless
  !the thick cloud section is invoked.
  DO  k=2,lp1
     DO  i=1,imax
        flxnet(i,k)=flxnet(i,k-1)+vsum1(i,k-1)
     END DO
  END DO
  DO  i=1,imax
     grnflx(i)=flxnet(i,lp1)
  END DO

  !this is the thick cloud section.optionally,if thick cloud
  !fluxes are to be "convectively adjusted",ie,df/dp is constant,
  !for cloudy part of grid point, the following code is executed.
  !first,count the number of clouds along the lat. row. skip the
  !entire thick cloud computation of there are no clouds.
  icnt=0
  DO  i=1,imax
     icnt=icnt+nclds(i)
  END DO
  IF (icnt == 0) go to 6999
  !find the maximum number of clouds in the latitude row
  kclds=nclds(1)
  DO  i=2,imax
     kclds=MAX(nclds(i),kclds)
  END DO

  !obtain the pressures and fluxes of the top and bottom of
  !the nc'th cloud (it is assumed that all ktop and kbtm's have
  !been defined!).
  DO  kk=1,kclds
     kmin=lp1
     kmax=0
     DO  i=1,imax
        j1=ktop(i,kk+1)
        j3=kbtm(i,kk+1)
        IF (j3 > j1) THEN
           ptop(i)=p(i,j1)
           pbot(i)=p(i,j3+1)
           ftop(i)=flxnet(i,j1)
           fbot(i)=flxnet(i,j3+1)
           !obtain the "flux derivative" df/dp (delptc)
           delptc(i)=(ftop(i)-fbot(i))/(ptop(i)-pbot(i))
           kmin=MIN(kmin,j1)
           kmax=MAX(kmax,j3)
        END IF
     END DO
     kmin=kmin+1
     !calculate the tot. flux chg. from the top of the cloud, for
     !all levels.
     DO  k=kmin,kmax
        DO  i=1,imax
           IF(ktop(i,kk+1) < k .AND. k <= kbtm(i,kk+1)) THEN
              z1(i,k)=(p(i,k)-ptop(i))*delptc(i)+ftop(i)
              !original flxnet(i,k)=flxnet(i,k)*(one-camt(i,kk+1)) +
              !original1 z1(i,k)*camt(i,kk+1)
              flxnet(i,k)=z1(i,k)
           END IF
        END DO
     END DO
  END DO
  !using this flux chg. in the cloudy part of the grid box, obtain
  !the new fluxes, weighting the clear and cloudy fluxes:again, only
  !the fluxes in thick-cloud levels will eventually be used.
  !     do 6051 k=1,lp1
  !     do 6051 i=1,imax
  !     flxnet(i,k)=flxnet(i,k)*(one-camt(i,nc)) +
  !    1            z1(i,k)*camt(i,nc)
  !051  continue
  !merge flxthk into flxnet for appropriate levels.
  !     do 1401 k=1,lp1
  !     do 1401 i=1,imax
  !     if (k.gt.itop(i) .and. k.le.ibot(i)
  !    1  .and.  (nc-1).le.nclds(i))  then
  !          flxnet(i,k)=flxthk(i,k)
  !     endif
  !401  continue
  !******end of cloud loop*****

6001 CONTINUE
6999 CONTINUE

  !the final step is to recompute the heating rates based on the
  !revised fluxes:
  DO  i=1,imax*l
     heatra(i,1)=radcon*(flxnet(i,2)-flxnet(i,1))*delp(i,1)
  END DO
  !the thick cloud section ends here.
END SUBROUTINE fst88
