SUBROUTINE lwr88(heatra,grnflx,topflx,press,temp,rh2o,qo3,cldfac,camt,&
     nclds,ktop,kbtm)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !subroutine lwr88 computes temperature-corrected co2 transmission
  !functions and also computes the pressure grid and layer optical
  !paths.
  !
  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
  INTEGER, PARAMETER :: ko2=12,ko21=ko2+1,ko2m=ko2-1

  INTEGER  :: i,k,kp

  REAL,INTENT(in),DIMENSION(imax,lp1)        :: press,temp
  REAL,INTENT(in),DIMENSION(imax,l)          :: rh2o,qo3
  REAL,INTENT(in out),DIMENSION(imax,l)      :: heatra
  REAL,INTENT(in out), DIMENSION(imax)       :: grnflx,topflx
  REAL,INTENT(in out),DIMENSION(imax,lp1)    :: camt
  REAL,INTENT(in out),DIMENSION(imax,lp1,lp1) :: cldfac
  INTEGER,INTENT(in out),DIMENSION(imax)     :: nclds
  INTEGER,INTENT(in out),DIMENSION(imax,lp1) :: ktop,kbtm

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

  REAL,DIMENSION(imax) :: emx1,emx2,vsum1,vsum2,a1,a2
  REAL,DIMENSION(imax,l) :: qh2o,delp,delp2,co2nbl,var1,var2,var3,var4,co2mr
  REAL,DIMENSION(imax,l) :: co2md,co2m2d,vv,vsum4
  REAL,DIMENSION(imax,lp1) :: t,p,co2sp1,co2sp2,cntval,toto3,tphio3,totphi
  REAL,DIMENSION(imax,lp1) :: totvo2,co2r,dift,co2r1,dco2d1,d2cd21,d2cd22,co2r2
  REAL,DIMENSION(imax,lp1) :: dco2d2,tdav,tstdav,vsum3,dco2dt,d2cdt2,texpsl
  REAL,DIMENSION(imax,lp1) :: tlsqu
  REAL,DIMENSION(imax,llp1) :: empl
  REAL,DIMENSION(imax,lp1,lp1) :: co21

  EQUIVALENCE (vsum3,tlsqu,texpsl)
  EQUIVALENCE (vv,vsum4)

  !compute flux pressures (p) and differences (delp2,delp)
  !compute flux level temperatures (t) and continuum temperature
  !corrections (texpsl)
  DO  k=2,l
     DO  i=1,imax
        p(i,k)=haf*(press(i,k-1)+press(i,k))
        t(i,k)=haf*(temp(i,k-1)+temp(i,k))
     END DO
  END DO
  DO  i=1,imax
     p(i,1)=zero
     p(i,lp1)=press(i,lp1)
     t(i,1)=temp(i,1)
     t(i,lp1)=temp(i,lp1)
  END DO
  DO  k=1,l
     DO  i=1,imax
        delp2(i,k)=p(i,k+1)-p(i,k)
        delp(i,k)=one/delp2(i,k)
     END DO
  END DO

  !compute argument for cont.temp.coeff.
  !(this is 1800.(1./temp-1./296.))
  DO  k=1,lp1
     DO  i=1,imax
        texpsl(i,k)=h18e3/temp(i,k)-h6p08108
        !then take exponential
        texpsl(i,k)=EXP(texpsl(i,k))
     END DO
  END DO

  !compute optical paths for h2o and o3, using the diffusivity
  !approximation for the angular integration (1.66). obtain the
  !unweighted values(var1,var3) and the weighted values(var2,var4).
  !the quantities h3m4(.0003) and h3m3(.003) appearing in the var2 and
  !var4 expressions are the approximate voigt corrections for h2o and
  !o3,respectively.
  DO  k=1,l
     DO  i=1,imax
        qh2o(i,k)=rh2o(i,k)*diffctr
        !vv is the layer-mean pressure (in atm),which is not the same as
        !the level pressure (press)
        vv(i,k)=haf*(p(i,k+1)+p(i,k))*p0inv
        var1(i,k)=delp2(i,k)*qh2o(i,k)*ginv
        var3(i,k)=delp2(i,k)*qo3(i,k)*diffctr*ginv
        var2(i,k)=var1(i,k)*(vv(i,k)+h3m4)
        var4(i,k)=var3(i,k)*(vv(i,k)+h3m3)
        !compute optical path for the h2o continuum, using roberts coeffs.
        !(betinw),and temp. correction (texpsl). the diffusivity factor
        !(which cancels out in this expression) is assumed to be 1.66. the
        !use of the diffusivity factor has been shown to be a significant
        !source of error in the continuum calcs.,but the time penalty of
        !an angular integration is severe.
        cntval(i,k)=texpsl(i,k)*rh2o(i,k)*var2(i,k)*betinw/ (rh2o(i,k)+rath2omw)
     END DO
  END DO

  !compute summed optical paths for h2o,o3 and continuum
  DO  i=1,imax
     totphi(i,1)=zero
     toto3(i,1)=zero
     tphio3(i,1)=zero
     totvo2(i,1)=zero
  END DO
  DO  k=2,lp1
     DO  i=1,imax
        totphi(i,k)=totphi(i,k-1)+var2(i,k-1)
        toto3(i,k)=toto3(i,k-1)+var3(i,k-1)
        tphio3(i,k)=tphio3(i,k-1)+var4(i,k-1)
        totvo2(i,k)=totvo2(i,k-1)+cntval(i,k-1)
     END DO
  END DO

  !emx1 is the additional pressure-scaled mass from press(l) to
  !p(l). it is used in nearby layer and emiss calculations.
  !emx2 is the additional pressure-scaled mass from press(l) to
  !p(lp1). it is used in calculations between flux levels l and lp1.
  DO  i=1,imax
     emx1(i)=qh2o(i,l)*press(i,l)*(press(i,l)-p(i,l))*gp0inv
     emx2(i)=qh2o(i,l)*press(i,l)*(p(i,lp1)-press(i,l))*gp0inv
  END DO

  !empl is the pressure scaled mass from p(k) to press(k) (index 2-lp1)
  !or to press(k+1) (index lp2-ll)
  DO  i=1,imax*l
     empl(i,2)=qh2o(i,1)*p(i,2)*(p(i,2)-press(i,1))*gp0inv
  END DO
  DO  i=1,imax*lm1
     empl(i,lp2)=qh2o(i,2)*p(i,2)*(press(i,2)-p(i,2))*gp0inv
  END DO
  DO  i=1,imax
     empl(i,1)=var2(i,l)
     empl(i,llp1)=empl(i,ll)
  END DO

  !compute weighted temperature (tdav) and pressure (tstdav) integrals
  !for use in obtaining temp. difference bet. sounding and std.
  !temp. sounding (dift)
  DO  i=1,imax
     tstdav(i,1)=zero
     tdav(i,1)=zero
  END DO
  DO  k=1,lp1
     DO  i=1,imax
        vsum3(i,k)=temp(i,k)-stemp(k)
     END DO
  END DO
  DO  k=1,l
     DO  i=1,imax
        vsum2(i)=gtemp(k)*delp2(i,k)
        vsum1(i)=vsum2(i)*vsum3(i,k)
        tstdav(i,k+1)=tstdav(i,k)+vsum2(i)
        tdav(i,k+1)=tdav(i,k)+vsum1(i)
     END DO
  END DO

  !evaluate coefficients for co2 pressure interpolation (a1,a2)
  DO  i=1,imax
     a1(i)=(press(i,lp1)-p0xzp8)/p0xzp2
     a2(i)=(p0-press(i,lp1))/p0xzp2
  END DO

  !perform co2 pressure interpolation on all inputted transmission
  !functions and temp. derivatives
  !successively computing co2r,dco2dt and d2cdt2 is done to save
  !storage (at a slight loss in computation time)
  DO  k=1,lp1
     DO  i=1,imax
        co2r1(i,k)=a1(i)*co231(k)+a2(i)*co238(k)
        d2cd21(i,k)=h1m3*(a1(i)*c2d31(k)+a2(i)*c2d38(k))
        dco2d1(i,k)=h1m2*(a1(i)*cdt31(k)+a2(i)*cdt38(k))
        co2r2(i,k)=a1(i)*co271(k)+a2(i)*co278(k)
        d2cd22(i,k)=h1m3*(a1(i)*c2d71(k)+a2(i)*c2d78(k))
        dco2d2(i,k)=h1m2*(a1(i)*cdt71(k)+a2(i)*cdt78(k))
     END DO
  END DO
  DO  k=1,l
     DO  i=1,imax
        co2mr(i,k)=a1(i)*co2m51(k)+a2(i)*co2m58(k)
        co2md(i,k)=h1m2*(a1(i)*cdtm51(k)+a2(i)*cdtm58(k))
        co2m2d(i,k)=h1m3*(a1(i)*c2dm51(k)+a2(i)*c2dm58(k))
     END DO
  END DO

  !compute co2 temperature interpolations for all bands,using dift
  !the case where k=1 is handled first. we are now replacing
  !3-dimensional arrays by 2-d arrays, to save space. thus this
  !calculation is for (i,kp,1)
  DO  kp=2,lp1
     DO  i=1,imax
        dift(i,kp)=tdav(i,kp)/tstdav(i,kp)
     END DO
  END DO
  DO  i=1,imax
     co21(i,1,1)=1.0
     co2sp1(i,1)=1.0
     co2sp2(i,1)=1.0
  END DO
  DO  kp=2,lp1
     DO  i=1,imax
        !calculations for kp>1 for k=1
        co2r(i,kp)=a1(i)*co251(kp,1)+a2(i)*co258(kp,1)
        dco2dt(i,kp)=h1m2*(a1(i)*cdt51(kp,1)+a2(i)*cdt58(kp,1))
        d2cdt2(i,kp)=h1m3*(a1(i)*c2d51(kp,1)+a2(i)*c2d58(kp,1))
        co21(i,kp,1)=co2r(i,kp)+dift(i,kp)*(dco2dt(i,kp)+  &
             haf*dift(i,kp)*d2cdt2(i,kp))
        !calculations for (effectively) kp=1,k>kp. these use the
        !same value of dift due to symmetry
        co2r(i,kp)=a1(i)*co251(1,kp)+a2(i)*co258(1,kp)
        dco2dt(i,kp)=h1m2*(a1(i)*cdt51(1,kp)+a2(i)*cdt58(1,kp))
        d2cdt2(i,kp)=h1m3*(a1(i)*c2d51(1,kp)+a2(i)*c2d58(1,kp))
        co21(i,1,kp)=co2r(i,kp)+dift(i,kp)*(dco2dt(i,kp)+  &
             haf*dift(i,kp)*d2cdt2(i,kp))
     END DO
  END DO

  !the transmission functions used in spa88 may be computed now.
  !(in the 250 loop,dift really should be (i,1,k), but dift is
  !invariant with respect to k,kp,and so (i,1,k)=(i,k,1))
  DO  k=2,lp1
     DO  i=1,imax
        co2sp1(i,k)=co2r1(i,k)+dift(i,k)*(dco2d1(i,k)+haf*dift(i,k)* d2cd21(i,k))
        co2sp2(i,k)=co2r2(i,k)+dift(i,k)*(dco2d2(i,k)+haf*dift(i,k)* d2cd22(i,k))
     END DO
  END DO

  !next the case when k=2...l
  DO  k=2,l
     DO  kp=k+1,lp1
        DO  i=1,imax
           dift(i,kp)=(tdav(i,kp)-tdav(i,k))/ (tstdav(i,kp)-tstdav(i,k))
           co2r(i,kp)=a1(i)*co251(kp,k)+a2(i)*co258(kp,k)
           dco2dt(i,kp)=h1m2*(a1(i)*cdt51(kp,k)+a2(i)*cdt58(kp,k))
           d2cdt2(i,kp)=h1m3*(a1(i)*c2d51(kp,k)+a2(i)*c2d58(kp,k))
           co21(i,kp,k)=co2r(i,kp)+dift(i,kp)*(dco2dt(i,kp)+  &
                haf*dift(i,kp)*d2cdt2(i,kp))
           co2r(i,kp)=a1(i)*co251(k,kp)+a2(i)*co258(k,kp)
           dco2dt(i,kp)=h1m2*(a1(i)*cdt51(k,kp)+a2(i)*cdt58(k,kp))
           d2cdt2(i,kp)=h1m3*(a1(i)*c2d51(k,kp)+a2(i)*c2d58(k,kp))
           co21(i,k,kp)=co2r(i,kp)+dift(i,kp)*(dco2dt(i,kp)+  &
                haf*dift(i,kp)*d2cdt2(i,kp))
        END DO
     END DO
  END DO

  !finally the case when k=kp,k=2..lp1
  DO  k=2,lp1
     DO  i=1,imax
        dift(i,k)=haf*(vsum3(i,k)+vsum3(i,k-1))
        co2r(i,k)=a1(i)*co251(k,k)+a2(i)*co258(k,k)
        dco2dt(i,k)=h1m2*(a1(i)*cdt51(k,k)+a2(i)*cdt58(k,k))
        d2cdt2(i,k)=h1m3*(a1(i)*c2d51(k,k)+a2(i)*c2d58(k,k))
        co21(i,k,k)=co2r(i,k)+dift(i,k)*(dco2dt(i,k)+ haf*dift(i,k)*d2cdt2(i,k))
     END DO
  END DO

  !we aren't doing nbl tfs on the 100 cm-1 bands .
  DO  k=1,l
     DO  i=1,imax
        co2nbl(i,k)=co2mr(i,k)+vsum3(i,k)*(co2md(i,k)+haf* vsum3(i,k)*co2m2d(i,k))
     END DO
  END DO

  !compute temp. coefficient based on t(k) (see ref.2)
  DO  k=1,lp1
     DO  i=1,imax
        IF (t(i,k) <= h25e2) THEN
           tlsqu(i,k)=b0+(t(i,k)-h25e2)* (b1+(t(i,k)-h25e2)*  &
                (b2+b3*(t(i,k)-h25e2)))
        ELSE
           tlsqu(i,k)=b0
        END IF
     END DO
  END DO

  !apply to all co2 tfs
  DO  k=1,lp1
     DO  kp=1,lp1
        DO  i=1,imax
           co21(i,kp,k)=co21(i,kp,k)*(one-tlsqu(i,kp))+tlsqu(i,kp)
        END DO
     END DO
  END DO
  DO  k=1,lp1
     DO  i=1,imax
        co2sp1(i,k)=co2sp1(i,k)*(one-tlsqu(i,1))+tlsqu(i,1)
        co2sp2(i,k)=co2sp2(i,k)*(one-tlsqu(i,1))+tlsqu(i,1)
     END DO
  END DO
  DO  k=1,l
     DO  i=1,imax
        co2nbl(i,k)=co2nbl(i,k)*(one-tlsqu(i,k))+tlsqu(i,k)
     END DO
  END DO
  CALL 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)
END SUBROUTINE lwr88
