C---------------------------------------------------------------------
      SUBROUTINE BORDE_RENE(teff,gsup,pre_atm,temp_atm,rad_atm,
     @                      dmass_atm)
C----------------------------------------------------------------------
C  Teff is the effective temperature in 10^6K and  gsup the surface
c  gravity (input data).
C
C It is calculated the pressure, temperature, outer linear radius, 
c and outer mass fraction at an optical depth of tau_ross= 25.1188 
c Model atmosphere: Rene Rohrmann (it considers Hummer & Mihalas formalism
c and it includes Lyman alfpha opacity.
C Tables are read as a function of Teff and gravity from  4 files, 
c which list radius, mass, T, and P at tau=25.11. 
C  Teff goes from 10,000 K down to 2500K with a step of 100 K
C  Log10 (g) goes from 9 down to 6.5 with step 0.1 
C  (in each table, the second column corresponds to log g = 9)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      common /borde_atm/ tempeff (76),dradio (76,26),
     @          dpresion(76,26),dmassatm(76,26),dtemp(76,26)
C
      DIMENSION OP(3),ABT(3),ABG(3),OP2(3), gravedad (26)
C
      SAVE IZ
      DATA IZ / 0 /
C
      do ii=1, 26
         gravedad(ii)= 9.1d0 - 0.1d0 * dfloat(ii) 
      end do
C
      IF(IZ.EQ.0) THEN
         IZ=1
         open(unit=999,file='P.dat')
         open(unit=998,file='T.dat')
         open(unit=997,file='R.dat')
         open(unit=996,file='M.dat')
c
         do II=1, 76
            READ(999,*) tempeff(ii), (dpresion(II,JJ),JJ=1,26) ! ln P_15
            READ(998,*) tempeff(ii), (dtemp(II,JJ),JJ=1,26)    ! ln T_6
            READ(997,*) tempeff(ii), (dradio(II,JJ),JJ=1,26)   ! ln r (ext.)
            READ(996,*) tempeff(ii), (dmassatm(II,JJ),JJ=1,26) ! ln(Matm/M*)
         end do
         CLOSE(UNIT =999)
         CLOSE(UNIT =998)
         CLOSE(UNIT =997)
         CLOSE(UNIT =996)
      ENDIF
C
      tempe= 1.d6 * teff
      dlogg= dlog10(gsup)
C
      if(tempe.gt.10000.d0.or.tempe.lt.2500.) then   
      write(*,*)'en BORDE_RENE la Teff esta fuera de tabla' 
         stop
      endif
C
      if(dlogg.gt.9.0.or.dlogg.lt.6.5) then   
      write(*,*)'en BORDE_RENE la gravedad esta fuera de tabla' 
         stop
      endif
C
C  Elegimos los 3 valores de temperatura mas cercanos al punto a interpo.
C      
      DO I = 3, 76
       IF(tempeff(I).LE.tempe) THEN
                           MT  = I - 2
                           GOTO 344
                          ENDIF
      END DO
C
C  Los 3 mas cercanos al log g
C
 344  DO IJ = 3 , 26
       IF(gravedad(IJ).LE.dlogg) THEN
                             MR  = IJ - 2
                             GOTO 180
                            ENDIF
      END DO
C
180   DO II = 1 , 3
       ABT(II) = tempeff(MT-1+II)
       ABG(II) = gravedad(MR-1+II)
      END DO
C
      DO II = MT , MT+2
       I = II-MT+1
       DO JJ = MR,MR+2
          OP(JJ-MR+1) = dpresion(II,JJ)
       END DO
C
C  Se interpola en log g. La interpolacion es en logaritmos
C
       CALL SISTEMA(OP,ABG,A,B,C)
       OP2(I) = A + (B + C * dlogg ) * dlogg
      END DO
C
      CALL SISTEMA(OP2,ABT,A,B,C)
      pre_atm   = A + ( B + C * tempe ) * tempe  ! presion (ln P_15)
C
      DO II = MT , MT+2
       I = II-MT+1
       DO JJ = MR,MR+2
          OP(JJ-MR+1) = dtemp(II,JJ)
       END DO
C
C  Se interpola en log g. La interpolacion es en logaritmos
C
       CALL SISTEMA(OP,ABG,A,B,C)
       OP2(I) = A + (B + C * dlogg ) * dlogg
      END DO
C
      CALL SISTEMA(OP2,ABT,A,B,C)
      temp_atm   = A + ( B + C * tempe ) * tempe ! ln T_6
c
      DO II = MT , MT+2
       I = II-MT+1
       DO JJ = MR,MR+2
          OP(JJ-MR+1) = dradio(II,JJ)
       END DO
C
C  Se interpola en log g. La interpolacion es en logaritmos
C
       CALL SISTEMA(OP,ABG,A,B,C)
       OP2(I) = A + (B + C * dlogg ) * dlogg
      END DO
C
      CALL SISTEMA(OP2,ABT,A,B,C)
      rad_atm   = A + ( B + C * tempe ) * tempe 
      rad_atm = dexp(rad_atm) / 1.d10  ! radio en 1.d10 cm 
C
      DO II = MT , MT+2
       I = II-MT+1
       DO JJ = MR,MR+2
          OP(JJ-MR+1) = dmassatm(II,JJ)
       END DO
C
C  Se interpola en log g. La interpolacion es en logaritmos
C
       CALL SISTEMA(OP,ABG,A,B,C)
       OP2(I) = A + (B + C * dlogg ) * dlogg
      END DO
C
      CALL SISTEMA(OP2,ABT,A,B,C)
      dmass_atm   = A + ( B + C * tempe ) * tempe
C
      RETURN
      END
C---------------------------------------------------------
      SUBROUTINE SISTEMA(AC,AB,A,B,C)
C---------------------------------------------------------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION AC(3),AB(3)
      COEF1 = (AC(1) - AC(3)) / (AB(1) - AB(3))
      COEF2 = (AC(2) - AC(3)) / (AB(2) - AB(3))
      C     = (COEF1 - COEF2) / (AB(1) - AB(2))
      B     = COEF2 - C * (AB(2) + AB(3))
      A     = AC(3) - (B + C * AB(3)) * AB(3)
      RETURN
      END