SUBROUTINE TAUFR1(IJ)
C     =====================
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/OPTDPT/DT(MDEPTH)
      dimension ST0(MDEPTH),SS0(MDEPTH),ab0(mdepth),
     *          tau(mdepth),taus(mdepth)
c      PARAMETER (TAUREF = 0.6666666666667)
      PARAMETER (TAUREF = 1.)
      PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
C
      FR=FREQ(IJ)
      DO ID=1,ND
         AB0(ID)=ABSO1(ID)
C        put in a floor to avoid division by zero:
         IF(AB0(ID)-SCAT1(ID).GT.1.D-100) THEN
           ST0(ID)=EMIS1(ID)/(AB0(ID)-scat1(id))
          ELSE
           ST0(ID)=0.
         ENDIF
         if(st0(id).eq.0) st0(id)=1.d-20*scat1(id)
         SS0(ID)=-SCAT1(ID)/AB0(ID)
      END DO
C
      id=1
      TAUMIN=ABSO1(1)*DEDM1
      tau(1)=taumin
C     to avoid a negative square root:
      taus(1)=sqrt(3.*ab0(id)*max(ab0(id)-scat1(id),0.d0))*DEDM1
C  
      IREF=1
      IREFs=1
      DO ID=1,ND-1
         DT(ID)=DELDMZ(ID)*(ABSOT(ID+1)+ABSOT(ID))
         tau(id+1)=tau(id)+dt(id)
C        to avoid negative square root:
         eps0=sqrt(ab0(id)*max(ab0(id)-scat1(id),0.d0))
         eps1=sqrt(ab0(id+1)*max(ab0(id+1)-scat1(id+1),0.d0))
         dts=deldm(id)*(eps0*dens1(id)+eps1*dens1(id+1))*sqrt(3.)
         taus(id+1)=taus(id)+dts
         IF(TAU(Id).LE.TAUREF.AND.TAU(Id+1).GT.TAUREF) IREF=Id
         IF(TAUs(Id).LE.TAUREF.AND.TAUs(Id+1).GT.TAUREF) IREFs=Id
      END DO
      if(iref.eq.1.and.tau(nd).le.tauref) iref=nd
      if(irefs.eq.1.and.taus(nd).le.tauref) irefs=nd
C 
      t0=1. 
      iref0=iref
      iref=irefs
      if(irefs.lt.nd) then 
         T0=LOG(TAUs(IREF+1)/TAUs(IREF))
         X0=LOG(TAUs(IREF+1)/TAUREF)/T0
         X1=LOG(TAUREF/TAUs(IREF))/T0
         DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1)
         TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1)
         abREF=EXP(LOG(ab0(IREF))*X0+LOG(ab0(IREF+1))*X1)
         scREF=EXP(LOG(scat1(IREF))*X0+LOG(scat1(IREF+1))*X1)
         STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1)
         tauef=EXP(LOG(TAU(IREF))*X0+LOG(TAU(IREF+1))*X1)
       else
         x0=1.
         x1=0.
         dmref=dm(nd)
         tref=temp(nd)
         abref=ab0(nd)
         scref=scat1(nd)
         stref=st0(nd)
         tauef=tau(nd)
      end if
      epref=(abref-scref)/abref
CX    add if statement to avoid overflow:
      IF(hk*fr/tref.lt.200.) then
         bref=1.4743e-2*(fr*1.e-15)**3/(exp(hk*fr/tref)-1.)
       ELSE
         bref=0.
      END IF
      taur=tauef*tauef
      if(tauef.gt.taur) taur=tauef
      yref=4.*ycon*tref*taur
      ALM=2.997925E18/FREQ(IJ)
      r1=rad(nfreq-kij(ij)+1,1)
      if(epref.ge.0) rb1=sqrt(epref)*bref
      if(epref.ge.0) rs1=sqrt(epref)*stref
C  
      return
      end