!This subroutine computes the PT profile for an irradiated plane-parallel atmosphere using the analytical solution of Parmentier & Guillot 2013.
! Units in the code are SI. (P in Pa, Kappa in kg/m^2, T in K, grav in m²/s etc...)
! This code is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (see http://creativecommons.org/licenses/by-nc-sa/3.0/ for more details).

      SUBROUTINE tprofile(Tint,Tirr,mu,grav,Gv1,Gv2,Betav,Gp,Beta,      &
     &Kappa,TAU,P,T,N)
      IMPLICIT NONE
      REAL*8, INTENT(IN) :: Tint !Temperature for the internal flux
      REAL*8, INTENT(IN) :: Tirr !Irradiation temperature
      REAL*8, INTENT(IN) :: mu !Stellar angle
      REAL*8, INTENT(IN) :: grav      !Gravity of the planet
      REAL*8, INTENT(INOUT) :: Gv1 !Gamma_v1=Kv1/Kr see Parmentier2013c
      REAL*8, INTENT(INOUT) :: Gv2 !Gamma_v2=Kv2/Kr see Parmentier2013c
      REAL*8, INTENT(INOUT) :: Betav ! Width of the first  visible band over the second one
      REAL*8, INTENT(INOUT) :: Gp  !Gamma_P=Kp/Kr
      REAL*8, INTENT(INOUT) :: Beta ! Width of the first thermal band over the second one
      INTEGER, INTENT(IN) :: N ! Number of levels in the atmosphere
      REAL*8, DIMENSION(N), INTENT(IN) :: P !Pressure
      REAL*8, DIMENSION(N), INTENT(INOUT) :: Kappa !Rosseland mean opacity used in the model. if OPA="CSTE", Kappa(i)=Kappa(1), if OPA="VALENCIA" Kappa is re-calculated in function of P and T.
      REAL*8, DIMENSION(N), INTENT(OUT) :: TAU ! Optical depth
      REAL*8, DIMENSION(N), INTENT(OUT) :: T ! Temperature

! Now all the internal variables
	REAL*8 :: aGp,bGp,cGp,aB,bB,cB,aGv1a,bGv1a,aGv1b,bGv1b,aGv1c,bGv1c&
     &,aGv1d,bGv1d,aGv2a,bGv2a,aGv2b,bGv2b,aGv2c,bGv2c,aGv2d,bGv2d,aBva &
     &,aBvb,aBvc,aBvd,bBva,bBvb,bbvc,Bbvd
	REAL*8 :: At1,At2,AA1v1,AA1v2,AA2v1,AA2v2
	REAL*8 :: a0,a1,a2v1,a2v2,a3v1,a3v2,b0,b1v1,b1v2,b2v1,b2v2,b3v1,  &
     &b3v2
	REAL*8 :: A,B,Cv1,Cv2,Dv1,Dv2,Ev1,Ev2
	REAL*8 :: R,Gv1s,Gv2s,TauLim,G1,G2 
      INTEGER :: i,j

!Error message if the input parameters are out of bounds
      if (Tint<0) then
       Print*, "Error : Tint must be positive"
       STOP
      endif

      if (Tirr<0) then 
       Print*, "Error : Tirr must be positive"
       STOP
      endif

      if (mu>1) then
       PRINT*, "Error : mu=cos(theta) must be smaller than one"
       STOP
      elseif (mu<0) then
       PRINT*, "Error : mu=cos(theta) must be positive"
       STOP
      endif

      if (grav<0) then
       Print*, "Error, gravity must be positive"
       STOP
      endif

! If the coefficients do not come from the fit, checking the values of the coefficients
       if (Gv1<0) then
        Print*, "Error : Gv1 must be positive"
        STOP
       endif

       if (Gv2<0) then
        Print*, "Error : Gv2 must be positive"
        STOP
       endif

       if (Betav<0) then
        Print*, "Error : Betav must be positive"
        STOP
       elseif (Betav>1) then 
        Print*, "Error : Betav must be lower than 1"
        STOP
       endif

       if (Gp<1) then
        Print*, "Error : Gp must be bigger than 1"
        STOP
       endif

       if (Beta<0) then
        Print*, "Error : Beta must be positive"
        STOP
       elseif (Beta>1) then 
        Print*, "Error : Beta must be lower than 1"
        STOP
       endif 

      !We calculate different useful quantities
      Gv1s=Gv1/mu
      Gv2s=Gv2/mu 
      R=1+(Gp-1)/(2*Beta*(1-Beta))+sqrt(((Gp-1)/(2*Beta*(1-Beta)))**2+  &
     &(Gp-1)/(2*Beta*(1-Beta)))
      G1=Beta+R-Beta*R
      G2=G1/R
      TauLim=1/(G1*G2)*sqrt(Gp/3.0)

      !Now we calculate the coefficients 
      At1=G1**2*log(1+1/(TauLim*G1))
      At2=G2**2*log(1+1/(TauLim*G2))
      AA1v1=G1**2*log(1+Gv1s/G1)
      AA2v1=G2**2*log(1+Gv1s/G2)
      AA1v2=G1**2*log(1+Gv2s/G1)
      AA2v2=G2**2*log(1+Gv2s/G2)
      a0=1/G1+1/G2

      a1=-1.0/(3*TauLim**2)*(Gp/(1-Gp)*(G1+G2-2)/(G1+G2)+(G1+G2)*TauLim-&
     &(At1+At2)*TauLim**2)!a1->infinity for Gp->1 but the product a1*b0 -> 0.

      a2v1=TauLim**2/(Gp*Gv1s**2)*((3*G1**2-Gv1s**2)*(3*G2**2-Gv1s**2)* &
     &(G1+G2)-3*Gv1s*(6*G1**2*G2**2-Gv1s**2*(G1**2+G2**2)))/(1-Gv1s**2* &
     &TauLim**2)

      a2v2=TauLim**2/(Gp*Gv2s**2)*((3*G1**2-Gv2s**2)*(3*G2**2-Gv2s**2)* &
     &(G1+G2)-3*Gv2s*(6*G1**2*G2**2-Gv2s**2*(G1**2+G2**2)))/(1-Gv2s**2* &
     &TauLim**2)

      a3v1=-TauLim**2*(3*G1**2-Gv1s**2)*(3*G2**2-Gv1s**2)*(AA2v1+AA1v1) &
     &/(Gp*Gv1s**3*(1-Gv1s**2*TauLim**2))

      a3v2=-TauLim**2*(3*G1**2-Gv2s**2)*(3*G2**2-Gv2s**2)*(AA2v2+AA1v2) &
     &/(Gp*Gv2s**3*(1-Gv2s**2*TauLim**2))

      b0=1.0/(G1*G2/(G1-G2)*(At1-At2)/3-(G1*G2)**2/sqrt(3*Gp)-(G1*G2)**3&
     &/((1-G1)*(1-G2)*(G1+G2)))

      b1v1=G1*G2*(3*G1**2-Gv1s**2)*(3*G2**2-Gv1s**2)*TauLim**2/         &
     &(Gp*Gv1s**2*(Gv1s**2*TauLim**2-1))

      b1v2=G1*G2*(3*G1**2-Gv2s**2)*(3*G2**2-Gv2s**2)*TauLim**2/         &
     &(Gp*Gv2s**2*(Gv2s**2*TauLim**2-1))

      b2v1=3*(G1+G2)*Gv1s**3/((3*G1**2-Gv1s**2)*(3*G2**2-Gv1s**2))
      b2v2=3*(G1+G2)*Gv2s**3/((3*G1**2-Gv2s**2)*(3*G2**2-Gv2s**2))
      b3v1=(AA2v1-AA1v1)/(Gv1s*(G1-G2))
      b3v2=(AA2v2-AA1v2)/(Gv2s*(G1-G2))

      A=1.0/3*(a0+a1*b0)
      B=-1.0/3*(G1*G2)**2/Gp*b0
      Cv1=-1.0/3*(b0*b1v1*(1+b2v1+b3v1)*a1+a2v1+a3v1)
      Cv2=-1.0/3*(b0*b1v2*(1+b2v2+b3v2)*a1+a2v2+a3v2)
      Dv1=1.0/3*(G1*G2)**2/Gp*b0*b1v1*(1+b2v1+b3v1)
      Dv2=1.0/3*(G1*G2)**2/Gp*b0*b1v2*(1+b2v2+b3v2)
      Ev1=(3-(Gv1s/G1)**2)*(3-(Gv1s/G2)**2)/(9*Gv1s*((Gv1s*TauLim)**2-1))
      Ev2=(3-(Gv2s/G1)**2)*(3-(Gv2s/G2)**2)/(9*Gv2s*((Gv2s*TauLim)**2-1))

!Now I calculate the temperature in function of Tau
       TAU(1)=Kappa(1)/grav*P(1)
      T(1)=(3.0*Tint**4/4*(TAU(1)+A+B*exp(-TAU(1)/TauLim))+3.0*Betav    &
     &*Tirr**4/4*mu*(Cv1+Dv1*exp(-TAU(1)/TauLim)+Ev1*exp(-Gv1s*TAU(1))) &
     &+3.0*(1-Betav)*Tirr**4/4*mu*(Cv2+Dv2*exp(-TAU(1)/TauLim)+Ev2*     &
     &exp(-Gv2s*TAU(1))))**0.25
       DO i=2,N
        TAU(i)=Tau(i-1)+Kappa(i)/grav*(P(i)-P(i-1))
      T(i)=(3.0*Tint**4/4*(TAU(i)+A+B*exp(-TAU(i)/TauLim))+3.0*Betav    &
     &*Tirr**4/4*mu*(Cv1+Dv1*exp(-TAU(i)/TauLim)+Ev1*exp(-Gv1s*TAU(i))) &
     &+3.0*(1-Betav)*Tirr**4/4*mu*(Cv2+Dv2*exp(-TAU(i)/TauLim)+Ev2*     &
     &exp(-Gv2s*TAU(i))))**0.25
       ENDDO


      END SUBROUTINE tprofile