#--------------------------------------------------------------
#Description:
#
#   This script uses the homebrew radiation code to find the
#   surface radiative cooling parameters for an atmosphere
#   containing a single greenhouse gas.  Mainly used for
#   the Venus case, but it can be used for any gas if the
#   appropriate exponential sums table is loaded.
#
#   Note: This uses miniClimt, which is in the Chapter 4
#   scripts directory. You will need to put a copy of miniClimt.py
#   in this directory, or in your modules directory, or anywhere
#   else Python can find it. 
#
#--------------------------------------------------------------

#Data on section of text which this script is associated with
Chapter = '6'
Section = '**'
Figure = '**'
#

from ClimateUtilities import *
import phys
import miniClimt as radmodel
import planets

#Path to the workbook datasets
datapath = '/Users/rtp1/Havsornen/Teaching/GeoSci232/WorkbookDatasets/'
#Path to the exponential sum tables
EsumPath = datapath + 'Chapter4Data/ExpSumTables/'

#--------------------Read in band data,set radmodel options-------------------
#Data for the full CO2 spectrum out to 10000 cm**-1. Computed
#from HITRAN using Venus surface pressure and temperature. Note
#that when we use the high-pressure table, we won't get accurate
#OLR, since the linear pressure scaling breaks down when there is
#so much line overlap.  Further, since the version of the radiation
#code used here doesn't do temperature scaling, the opacity of the
#upper atmosphere is greatly over-estimated. Even if we used
#the temperature scaling in miniClimtFancy, accurate OLR might
#not be obtained using the crude scaling over such a wide temperature
#range. 
#
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableExtendedWave.730K.92bar.Self.data')
radmodel.BackgroundGas = phys.N2 #The transparent background gas (Not used for q=1)
radmodel.GHG =phys.CO2 #The greenhouse gas
#
radmodel.pref = 92.e5 #Reference pressure at which band data is given
radmodel.Trans = radmodel.TransEsums
#radmodel.Trans = TransMalkmus   (Don't use this with solar absorption)
#radmodel.Trans = radmodel.TransOobleck (Don't use this with solar absorption

#Set the continuum, if there is one
radmodel.Continuum = radmodel.CO2Continuum
#radmodel.Continuum = radmodel.NoContinuum

#==============Now do the calculation====================================

#--------------Initializations-------------------------------------------
    



#-------Set up the pressure array------------------------------

n = 200
ps,ptop = 92.e5,100.
#p = radmodel.setpExtraGroundRes(ps,ptop,n)
p = radmodel.setpLin(ps,ptop,n)
#p = radmodel.setpLog(ps,ptop,n) #Puts extra resolution in stratosphere.
#-------------------------------------------------------------------------


#------------Set the greenhouse gas concentration------------------------
#q = (44./29.)*300.e-6 #Earthlike CO2
q = 1. #Use for Mars and Venus. When q=1, the background gas isn't used
       #In this case, remember to use the self-broadened table,
       #not the foreign-broadened table.
Rcp = radmodel.GHG.Rcp #Change if the GHG isn't dominant
#------------------------------------------------------------------------


#--------------Other parameters-------------------------------------------

radmodel.g = planets.Venus.g #Set the gravity
radmodel.LWCutoff = 5000. #Speeds up flux calculation.
                          #OK for temperatures under 1000K


Tg = 727.
T = Tg*(p/p[-1])**Rcp
#

#---------------------------------------------------------
#--------------Initializations Done-----------------------
#--------------Now do the calculation-------------------
#Add in contribution of part of spectrum outside the band table
nu1 = radmodel.bandData[0].nu1
nu2 = radmodel.bandData[-1].nu2
m = romberg(radmodel.Planck) #Evaluates definite integral

flux1 = m([.001,nu1],Tg) + m([nu2,50000.],Tg)
fluxLW,heatLW = radmodel.radCompLW(p,T,Tg,q)
print "Tg =",Tg,"ps =",ps/1.e5
print fluxLW[0],'%e'%fluxLW[-1],(fluxLW[-1])/(phys.sigma*Tg**4)
print "With out of band:",flux1
print (fluxLW[-1]+flux1)/(phys.sigma*Tg**4)