!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