PROGRAM paper2 IMPLICIT NONE !!!!!!INPUT VARIABLES!!!!!!!!!!!!!!!! !Model parameters REAL*8 :: Pmin,Pmax !minimal and maximal pressure of the model. INTEGER, PARAMETER :: N=100 ! Number of levels in the model CHARACTER*150 :: OutputFile ! Name of the output file CHARACTER*150 :: OutputFileCSV ! Name of the output file !Planet parameters REAL*8 :: Tint !Temperature for the internal flux REAL*8 :: Teq0 !Temperature of the incoming flux REAL*8 :: mu !Cosine of the stellar angle mu=1/sqrt(3) for averaged profiles REAL*8 :: f !Parameter for averaged profiles f=0.5 for dayside average and f=0.25 for planet-average REAL*8 :: grav !Gravity of the planet (m/s**2) !Options CHARACTER*10 :: ROSS ! Can be either 'USER' when Kappa is defined by the user or 'AUTO' to use the opacities of Freedman et al. 2008 as fitted by Valencia et al. 2013 CHARACTER*10 :: COEFF ! Can be either 'USER' when the input coefficients for Gv1,Gv2,Betav,Beta,Gp or 'AUTO' to use the fit provided by Parmentier et al. 2014 CHARACTER*10 :: COMP ! Can be 'SOLAR' for a solar-composition atmosphere or 'NOTIO' for an atmosphere without TiO/VO. CHARACTER*10 :: STAR ! Can only be 'SUN' for now. CHARACTER*10 :: ALBEDO !If 'USER' uses the value of Ab, if 'AUTO', uses the fit of Parmentier at al. 2014. CHARACTER*10 :: CONV !If YES calculates a radiative/convective profile, if NO, calculates only the radiative solution !Atmosphere parameters REAL*8 :: Gv1 !Gamma_v1=Kv1/Kr is the ratio of the first visible opacity to the Rosseland mean opacity see Parmentier et al. 2014 REAL*8 :: Gv2 !Gamma_v2=Kv2/Kr is the ratio of the second visible opacity to the Rosseland mean opacity see Parmentier et al. 2014 REAL*8 :: Gv3 !Gamma_v2=Kv3/Kr is the ratio of the third visible opacity to the Rosseland mean opacity see Parmentier et al. 2014 REAL*8 :: Gp !Gamma_P=Kp/Kr ratio of the Planck to the Rosseland mean opacity see Parmentier et al. 2014 REAL*8 :: Beta ! Relative spectral width of the first thermal band REAL*8 :: Ab ! Bond albedo of the planet !!!!!!INPUTS/OUTPUTS!!!!!!!!!! REAL*8, DIMENSION(N) ::Kappa ! Rosseland mean opacity at each atmospheric level, either calculated or set by the user !!!!!!!OUTPUTS!!!!!!!!!!!!! REAL*8 :: Teff ! Effective temperature of the model REAL*8, DIMENSION(N) :: TAU ! Rosseland optical depth REAL*8, DIMENSION(N) :: P ! Pressure REAL*8, DIMENSION(N) :: T ! Temperature REAL*8, DIMENSION(N) :: gradad ! Adiabatic gradient REAL*8, DIMENSION(N) :: gradrad ! Radiative gradients REAL*8 :: PRC ! Pressure at the radiative/convective boudary !!!!!!!INTERNAL VARIABLES INTEGER :: i,K CHARACTER*4 :: xtemp OutputFile='PTprofile.dat' OutputFileCSV='PTprofile.csv' !!!!!!Model parameters Pmin=1e1 !Minimal pressure (Pa) Pmax=1e9 ! Maximal pressure (Pa) !!!!!!Planet parameters Tint=100!Internal temperature (K) Teq0=1000!!Equilibrium temperature for zero albedo (K) f=0.25 !f=0.25 for a planet-average profile, f=0.5 for a dayside average profile mu=1/sqrt(3.0) ! Cosine of the irradiation angle. mu=1/sqrt(3) is the mean mu to use for a dayside or a planet average. ! Teff0=1000 !Effective temperature for zero albedo grav=25! gravity (m/s**2) !!!!!Options ROSS="AUTO" !This can be AUTO or USER. If USER, the Rosseland mean opacity must be provided below for each level. If AUTO, the Rosseland opacities are calculated from the fit of Valencia et al. 2013 COEFF="AUTO" !This can be AUTO or USER whether you want to use the fit of the coefficients provided in Parmentier et al. 2013 or set your own coefficients. STAR="SUN" COMP="SOLAR" !This can be SOLAR for a solar-composition atmosphere or NOTIO for an atmosphere without TiO/VO ALBEDO="AUTO" !This can be AUTO or USER. If USER, the Bond albedo must be provided below. If AUTO, the Bond Albedo is calculated from the fit of Parmentier et al. 2014 CONV="YES" ! If "YES", a convective zone is calculated at the bottom of the model using the convective gradient of Parmentier et al. 2014. If "NO" only the radiative solution is provided !!!!!Input variables, used only if COEFF="USER" if(COEFF.eq."USER") then Gv1=10 Gv2=0.1 Gv3=0.01 ! 0.5 is two bands of equal width for the absorption of the stellar flux Gp=10! Gp=10 Planck opacities are ten times higher than Rosseland ones. beta=0.5! The two thermal bands have the same width endif !!!!!Input variables, used only if ALBEDO="USER" if(ALBEDO.eq."USER") then Ab=0.3!!Bond albedo set by the user if Albedo="USER". endif !!!!!Input variables, used only if ROSS="USER" if(ROSS.eq."USER") then DO i=1,N Kappa(i)=1e-1 ! kg/m^2 Here I set the Rosseland mean opacity, if ROSS="CONSTANT". Kappa(i) can be any functional form of the pressure P(i) ENDDO endif !!!!!!!!!!!!!!END OF INPUTS!!!!!!!!!!!!!!!!!!!!! !Set the pressure levels DO i=1,N P(i)=10**(log10(Pmin)+(log10(Pmax)-log10(Pmin))* & &1.0*(i-1)/(N-1))!The pressure varies from Pmin to Pmax in logspace with N points! ENDDO !Calculate the Temperature call tprofile(Tint,Teq0,mu,f,grav,Gv1,Gv2,Gv3,Gp,Beta,Kappa, & &Ab,Teff,ROSS,COEFF,COMP,STAR,CONV,ALBEDO,TAU,P,T,N,PRC,gradad, & &gradrad) !Output formats 1400 format (f10.2,1X,f10.2,1X,f10.2,1X,f10.2,1X,f10.2,1X,f6.2,1X, & &e9.3,1X,e9.3,1X,f5.4,1X,e9.3,1X,f5.4,1X,f5.4,1X & &,A10,A10,A10,A10,A10,A10' Teff,Tint,Teq0,mu,f,grav,Gv1,Gv2,Gv3,Gp& &,Beta,Ab,ROSS,COEFF,COMP,STAR,CONV,ALBEDO') 1500 format( & & 'c',t1,'(1)N',t6,'(2)P(Pa)',t17,'(3)T(K)' & & ,t27,'(4)Tau', t37,'(5)Kr(m^2/kg)', t50, '(6)Gradad' & & , t60, '(7)Gradrad', t72,'(8)PRC (Pa)') 9900 format (1X,I3,",",E9.3,",",F10.3,",",E9.3,",",E9.3,",",F8.5,"," & &,F8.5,",",E9.3,",",E9.3,",",E9.3,",",F5.3,",",E9.3,",",F5.3,",", & &F6.3,",",F6.1,",",F6.1,",",F6.1,",",F5.3,",",F5.3,",",F8.3,",", & &A10,",",A10,",",A10,",",A10,",",A10,",",A10) 9800 format (1X,I3,1X,E9.3,1X,F10.3,1X,E9.3,1X,E9.3,1X,E9.3,1X, & &E9.3,1X,E9.3,1X,E9.3,1X,E9.3,1X,E9.3,1X,E9.3) write(xtemp,'(I4.4)') Int((1-Ab)**0.25*(4*f)**0.25*Teq0) ! OutputFile='./profilesNoTiO/PTprofile-Tmu-'//xtemp//'.dat' !Open and write output file open(unit=40,form='formatted',status='UNKNOWN',file=OutputFileCSV& &) write(40,'(A)')'N,P(Pa),T(K),Tau,Kappa(m^2/kg),gradad,gradrad,PRC& &,Gv1,Gv2,Gv3,Gp,Beta,Ab,Teff(K),Tint(K),Teq0(K),mu,f,grav(m/s^2), & &ROSS,COEFF,COMP,STAR,CONV,ALBEDO' DO i=1,N write(40,9900)i,P(i),T(i),Tau(i),Kappa(i),gradad(i),gradrad(i),& &PRC,Gv1,Gv2,Gv3,Gp,Beta,Ab,Teff,Tint,Teq0,mu,f,grav,ROSS,COEFF, & &COMP,STAR,CONV,ALBEDO ENDDO close(40) open(unit=40,form='formatted',status='UNKNOWN',file=OutputFile) write(40,1400) Teff,Tint,Teq0,mu,f,grav,Gv1,Gv2,Gv3,Gp,Beta,Ab,& &ROSS,COEFF,COMP,STAR,CONV,ALBEDO write(40,*),' ' write(40,*),' ' write(40,1500) DO i=1,N write(40,9800) i,P(i),T(i),Tau(i),Kappa(i),gradad(i),gradrad(i)& &,PRC ENDDO close(40) ! ENDDO write(*,*) 'Done, PTprofile is available in files: "', & &TRIM(OutputFile),'" or "',TRIM(OutputFileCSV),'"' END PROGRAM paper2