PROGRAM     paper1
      IMPLICIT NONE
      REAL*8 :: Tint
      REAL*8 :: Tirr
      REAL*8 :: mu
      REAL*8 :: grav      
      REAL*8 :: Gv1
      REAL*8 :: Gv2
      REAL*8 :: Betav
      REAL*8 :: Gp
      REAL*8 :: Beta
      REAL*8 :: Pmin,Pmax 
      REAL*8 :: bid
      INTEGER, PARAMETER :: N=100 ! Number of levels in the atmosphere
      INTEGER :: i
      REAL*8, DIMENSION(N) :: TAU
      INTEGER*4 :: leng
      CHARACTER*150 :: OutputFile
      REAL*8, DIMENSION(N) :: P
      REAL*8, DIMENSION(N) :: T
      REAL*8, DIMENSION(N) ::Kappa

      ! Default values
      Tint=100
      Tirr=2000
      mu=1/sqrt(3.0) ! Mean mu for a planet or dayside average
      grav=25 !SI units
      Gv1=1 !Kv1=Kr/10
      Gv2=1 ! Kv2=Kr/10
      Betav=0.5 ! beta=0.5 is two bands of equal width for the absorption of the stellar flux
      Gp=100 ! Gp=10 Planck opacities are ten times higher than Rosseland ones.
      beta=0.5 ! The two thermal bands have the same width
      OutputFile="ptprofile.csv"

!Set the pressure levels
      Pmin=1e-2 !Minimal pressure
      Pmax=1e10 ! Maximal pressure
      DO i=1,N
       P(i)=10**(log10(Pmin)+(log10(Pmax)-log10(Pmin))*                 &
     &1.0*(i-1)/(N-1))!The pressure varies from taumin to taumax in logspace with N points
      ENDDO

      DO i=1,N
      Kappa(i)=1e-1 ! kg/m^2 Rosseland mean opacity at each level i. Any function of the pressure P(i) can be used here.
      ENDDO

      !Calculate the Temperature 
      call tprofile(Tint,Tirr,mu,grav,Gv1,Gv2,Betav,Gp,Beta,Kappa,      &
     &TAU,P,T,N)
9800   format(ES9.3,12(",",ES9.3))
      !Output
       open(unit=40,form='formatted',status='UNKNOWN',file=OutputFile)
       write(40,'(A)')'P(Pa),T(K),Tau,OpaRosseland(m^2/kg),beta,Gp,Betav&
     &,Gv1,Gv2,grav(m/s^2),mu,Tirr(K),Tint(K)'
        DO i=1,N
         write(40,9800)P(i),T(i),Tau(i),Kappa(i),beta,Gp,Betav,Gv1,Gv2  &
     &,grav,mu,Tirr,Tint
        ENDDO
       close(40)
      write(*,*)'Done, output in file ',Outputfile
      END PROGRAM paper1