#------------------------------------------------------------
#Reads temperature and humidity sounding data from an
#ascii file and computes the infrared fluxes using ClimT.
#
#**UNDER CONSTRUCTION
#  (Not used in the text or problems but a handy utility to
#   have someday)
#ToDo:
#      **Set it up to use CEPEX or IGRA soundings using readTabele
#        and interpolation

#Bug: qsat is computed incorrectly: The denominator should have
#the actual pressure, not the surface pressure.  The code with the
#bug in it was used to make the first figure in the RealClimate
#article "A Busy Week for Water Vapor."  The changes are not
#extreme, however, and since we only described the calculation
#as based on an idealized moisture profile, the main point doesn't
#change. Upper levels were somewhat drier than they should have been.
#The OLR changes from 311 W/m**2 to 297 W/m**2
#
#
#Query: How does the NCAR radiation code (and ClimT) make
#use of the surface temperature? The downwelling flux
#seems to be sensitive to surface temperature, whereas it
#shouldn't be.

import climt_lite as climt
from ClimateUtilities import *
import phys


#---Instantiate the radiation model---
r = climt.radiation()

#Set global constants
ps = 1000.
rh = .5 #Relative humidity aloft
rhbdd = .8 #Relative humidity in b.l.
#
#Read in data and interpolate to the right number of
#levels, if necessary. Remember, ClimT uses g/kg as the
#units for water vapor.
#
#(Temporary hard-coded array of pressure/temperature sounding)
#**ToDo: Replace this with read of sounding data, including humidity
pin = [ 1000.,   925.,   850.,   700.,   600.,   500.,   400.,   300.,   250.,
              200.,   150.,   100.,    70.,    50.,    30.,    20.,    10.]
Tin = [ 299.850006,  299.929993,  295.299988,  280.799988,  271.530029,  263.480011,
             252.030014,  235.780014,  225.650009,  214.980011,  212.350006,
             213.630005,  216.400024,  218.930023,  222.950012,  228.650009,
             232.750015]
pin.reverse() #Use the ordering ccmrad wants
Tin.reverse() #Use the ordering ccmrad wants
pin = numpy.array(pin)
Tin = numpy.array(Tin)

#
#---Set up pressure array----
#(was typecode = 'd'). Default float is a double
p = ( numpy.arange(r.nlev, dtype= 'float' )+ 0.1 ) * ps/r.nlev


#Interpolate temperature
pI = interp(pin,Tin)
T = [pI(pp) for pp in p]



#Create an idealized water vapor profile
#ToDo: Read in actual humidity data from sounding
T = numpy.array(T)
#Old version: qsat = [622.*phys.satvpg(TT)/1.e5 for TT in T] #q in g/Kg !
qsat = [phys.satvpg(TT) for TT in T] 
qsat = numpy.array(qsat)
qsat = 622.*qsat/(100.*p) #q in g/Kg !
q = numpy.where(p < 200.,1.e-9,rh*qsat)
q = numpy.where(p>=850.,rhbdd*qsat,q)



#Set parameters
co2 = 4.*280.
ch4 = 1.e-30
Ts = T[-1]


#Do the radiation calculation
r(p=p,ps=ps,T=T,Ts=Ts,q=q,co2=co2,ch4=ch4)
#Put out the results
print "OLR:",-r.lwflx[0]
print "Surface Back Radiation:", r.lwflx[-1]+ phys.sigma*Ts**4