;Developed by Wayne Feltz  3/15/95
;Space Science and Engineering Center (SSEC) University of Wisconsin

;  This idl program was developed to read aeri retrieval netcdf files 
;and produce atmospheric time cross sections.  These cross sections may
;be temperature, potential temperature, virtual temperature, mixing ratio, and
;relative humidity.  Total precipitable water is also calculated and stored in
;an array.  This can be ploted at the idl prompt after aerixsec has been
;terminated.  (Useful for comparing to microwave radiometer data)

;If any questions please contact Wayne Feltz, waynef@ssec.wisc.edu or
;608-265-6283

;The following modules are used with this program: 
;filter,hourform,indaxe,relhum,satmix,svpice,svpwat,totpre,virtemp,vline2,wvmr
;***************************************************************************
;Works with no ozone brf file  7/7/95 WFF
;***************************************************************************


COMMON Share1, p, time
close, 1
close, 2

set_plot, 'X'
!P.multi=[0,1,2,0,0]
!x.margin = [8,6]
!y.margin = [4,4]
font=1
window, 0, colors=256, xsize = 700, ysize = 700, retain=2
;average span

nlev=40
cin=255

;read data
as1=1
as2=2

;date=intarr(3)


for m=1,2 do begin
   if m eq 1 then filename='/home3/waynef/camexiii/andros/cdffile/980808pr.cdf'
   if m eq 2 then filename='/home3/waynef/camexiii/raman/wv0808co.007'


if m eq 1 then begin
   ;#############################################################
   ;  READ IN NETCDF AERI PHYSICAL RETRIEVALS

   fid = ncdf_open(filename)
   
   ncdf_varget, fid, 'height', alt
   ncdf_varget, fid, 'pressure', pr
   ncdf_varget, fid, 'ambientTemp', temp
   ncdf_varget, fid, 'dewpointTemp', dtemp
   ncdf_varget, fid, 'AERItime', time
   ncdf_varget, fid, 'cloudBaseHeight', cbh
   ncdf_varget, fid, 'Date', date

   
   stime=string(time)
   recno=n_elements(time)
   sdate=strcompress(string(date),/remove_all)
   recst=0
   pr=reverse(pr(0:19,*))
   temp=reverse(temp(0:19,*))
   dtemp=reverse(dtemp(0:19,*))
   
   alt=reverse(alt)
   mixr=fltarr(20,recno)
   for i = 0, recno - 1 do mixr(*,i)=wvmr(pr(*,i), $
          temp(*,i),dtemp(*,i))
          
          
;#############################################################

recst=115
recend=145

stop
minalt = 0.0
maxalt = max(alt(where(temp(0:19,0) GT 0.0)))

top=3000.
bot=0.0
ibot=19
itop=where(alt eq  top)
itop = fix(fix(itop),0)
height = 20 - where(alt eq  top)
height = fix(fix(height),0)

;matrix to contour
surfta = fltarr(recno, height)
surftha = fltarr(recno, height)
surfma = fltarr(recno, height)
surftva = fltarr(recno, height)
surfrha = fltarr(recno, height)
surfva_ns= fltarr(recno, height)

;first record
p = pr(itop:ibot,recst)
al=alt(itop:ibot)

surfta(recst, *) = temp(itop:ibot,recst)
surfma(recst,*) = mixr(itop:ibot,recst)

;Interpolate AERI temperature and mixing ratio data
for i = recst+1, recend do begin 
    surfta(i, *) = temp(itop:ibot,i) 
    surfma(i,*) = mixr(itop:ibot,i)      
endfor

result=where(surfta eq -8888.0,count)
if (count gt 1) then surfta(result)=-500.0
if (count gt 1) then surfma(result)=-500.0

for i = recst, recend do surfva_ns(i,*)=virtemp(surfta(i,*),surfma(i,*))-273.16

;smooth surf with a running average filter
for i = 0, height-1 do $  
  surfta(recst:recend,i) = filter(surfta(recst:recend, i), as1, (recend+1)-recst)
for i = 0, height-1 do $  
  surfma(recst:recend,i) = filter(surfma(recst:recend, i), as1, (recend+1)-recst)
  
;Calculate total preciptable water for later use
;Calculate theta for AERI cross section
preh2o=fltarr(recno)
pbot = (max(pr(*,0)))-1
for i = recst, recend do begin
   totpre,p, surfma(i,*), height, aeri_ini_temp,aeri_total_temp
   surftha(i,*)=surfta(i,*)*(pbot/p)^.286
   surftva(i,*)=virtemp(surfta(i,*),surfma(i,*))
   surfrha(i,*)=relhum(p,surfta(i,*),surfma(i,*))
   preh2o(i)=aeri_total_temp  
endfor 
if (count gt 1) then surfrha(result)=-500.0
  
;pressure axes (2)
altl = [0.0, 250.0, 500.0, 750.0, 1000.0, 1250.0, 1500.0, 1750.0, 2000.0, $
       2250.0, 2500.0]
nintery = (size(altl))(1)
ytname = make_array(nintery, /string)
FOR i = 0, nintery-1 DO $
  ytname(i) = string(altl(i), format = "(f6.1)") 

;desired time ticks (2)


st = 2*(fix(min(time(recst:recend))/2.0)+1)
et = 2*(fix(max(time(recst:recend))/2.0))
noet = ((et-st)/2.0)+1 	
noat = noet + 2
h = fltarr(noat)
h(0)=min(time(recst:recend))
h(1)=st
for i= 2, noat-2 do h(i)=h(i-1)+2.0
h(noat-1)=max(time(recst:recend))
ninterx = (size(h))(1)
;time tick name
xtname = make_array(ninterx, /string)
FOR i = 0, ninterx-1 DO $
  xtname(i)=string(hourform(h(i)), format="(f5.1)")
;index of the desired time tick in the raw x range [0,nprof-1]
indx = indaxe(time, h, ninterx, 1)

;***  end of netcdf file data manipulation **

endif



;***************Beginning of CART LIDAR file read and data**********
if m eq 2 then begin
   ;#############################################################
   ;  READ IN GSFC RAMAN LIDAR RETRIEVALS

;number of pressure levels
nlev = 186
;max number of profiles
dummy = ''
data = fltarr(6, nlev)
;read data
filename = pickfile(/fix,  $
       path='/home3/waynef/camexiii/raman/')
openr, 1, filename

;read total number of profiles
readf, 1, nprof
time = fltarr(nprof)

;read the first record

readf, 1, num, t1, format = '(I1,8x,f8.5)'
for j = 0, 2 do readf, 1, dummy, format='(a1)'
readf, 1, data

time=fltarr(nprof)
alt=fltarr(nlev)
mixr=fltarr(nlev,nprof)

time(0)=t1
alt(*)=data(2,*)
alt=reverse(alt)*1000
mixr(*,0)=data(3,*)

for i = 1, nprof-1 do begin
    readf, 1, num, t1, format = '(I1,8x,f8.5)'
    for j = 0, 2 do readf, 1, dummy, format='(a1)'
    readf, 1, data
    time(i)=t1
    mixr(*,i)=reverse(data(3,*),2)
endfor
close, 1    


stime=string(time)
recno=nprof
olev=n_elements(alt)
sdate=strcompress(string(date),/remove_all)
recst=0






;   mixr=reverse(mixr(0:olev-1,*))
;   alt=reverse(alt)*1000
;   alt(olev-1)=0.0
   
;#############################################################

recst=0
recend=recno-1


minalt = 0.0
maxalt = max(alt)

top=2500.
bot=0.0
ibot=olev-1
itop=32
itop = fix(fix(itop),0)
height = olev - itop
height = fix(fix(height),0)

;matrix to contour

surfma = fltarr(recno, height)

;first record

al=alt(itop:ibot)

surfma(recst,*) = mixr(itop:ibot,recst)

;Interpolate AERI temperature and mixing ratio data

for i = recst+1, recend do surfma(i,*) = mixr(itop:ibot,i)      


result=where(surfma eq -9999.0,count)

if (count gt 1) then surfma(result)=-500.0

;smooth surf with a running average filter

;for i = 0, height-1 do $  
;  surfma(recst:recend,i) = filter(surfma(recst:recend, i), as2, (recend+1)-recst)
  
;Calculate total preciptable water for later use
;Calculate theta for AERI cross section
;preh2o=fltarr(recno)
;pbot = (max(pr(*,0)))-1
;for i = recst, recend do begin
;   totpre,p, surfma(i,*), height, aeri_ini_temp,aeri_total_temp
;   preh2o(i)=aeri_total_temp  
;endfor 
;if (count gt 1) then surfrha(result)=-500.0
  
;pressure axes (2)
altl = [0.0, 250.0, 500.0, 750.0, 1000.0, 1250.0, 1500.0, 1750.0, 2000.0, $
       2250.0, 2500.0]
nintery = (size(altl))(1)
ytname = make_array(nintery, /string)
FOR i = 0, nintery-1 DO $
  ytname(i) = string(altl(i), format = "(f6.1)") 

;desired time ticks (2)


st = 2*(fix(min(time(recst:recend))/2.0)+1)
et = 2*(fix(max(time(recst:recend))/2.0))
noet = ((et-st)/2.0)+1 	
noat = noet + 2
h = fltarr(noat)
h(0)=min(time(recst:recend))
h(1)=st
for i= 2, noat-2 do h(i)=h(i-1)+2.0
h(noat-1)=max(time(recst:recend))
ninterx = (size(h))(1)
;time tick name
xtname = make_array(ninterx, /string)
FOR i = 0, ninterx-1 DO $
  xtname(i)=string(hourform(h(i)), format="(f5.1)")
;index of the desired time tick in the raw x range [0,nprof-1]
indx = indaxe(time, h, ninterx, 1)

;***  end of netcdf file data manipulation **


endif

;****************end of CART LIDAR data read ********* 
;repeat begin

; Pick variable field to be contoured
field=1      
;print,  'Temperature (1)'      
;print,  'Potential Temperature (2)'
;print,  'Mixing Ratio (3)'
;print,  'Virtual Temperature (4)'
;print,  'Relative Humidity (5)'

;repeat begin
;	read, prompt='Please pick first field to be contoured:',field1
;endrep until (field ge 1) and (field le 5)
field1=3
if (field1 eq 2) then $
	for j = recst, recend do $
		for i = height-11, height-2 do $
			if (surftha (j,i) lt surftha(j,i+1)) then $
			surftha(j,i+1)=surftha(j,i)
if m eq 1 then begin
case field1 of
	1:Begin
	   xsec1=surfta
	   type1='!17Wisconsin AERI Sounding System !C!DTEMPERATURE  (K) '
	   label='t'
	  End
	2:Begin
		xsec1=surftha
		type1='!17Wisconsin AERI Sounding System !C!DPOTENTIAL TEMPERATURE  (K) '
	   label='th'		
	  End
	3:Begin
	   xsec1=surfma
		type1='!17Wisconsin AERI Sounding System !C!DMIXING RATIO (G/KG) '
      minl1=0.0
      dc1=1.0
      ncont1=22
      label='mr'
	  End
	4:Begin
	   xsec1=surftva
		type1='!17Wisconsin AERI Sounding System !C!DVIRTUAL TEMPERATURE (K) '
		label='tv'
	  End	
	5:Begin
	   xsec1=surfrha
		type1='!17Wisconsin AERI Sounding System !C!DRelative Humidity (%) '
      minl1=10
      dc1=5
      ncont1=18
      label='rh'	
	  End	    
endcase
endif

if m eq 2 then begin
field1=3
case field1 of
	1:Begin
	   xsec1=surfta
	   type1='!17RADIOSONDE !C!DTEMPERATURE  (K) '
	   label='t'
	  End
	2:Begin
		xsec1=surftha
		type1='!17RADIOSONDE !C!DPOTENTIAL TEMPERATURE  (K) '
	   label='th'		
	  End
	3:Begin
	   xsec1=surfma
		type1='!17GSFC RAMAN LIDAR (SRL) !C!DMIXING RATIO (G/KG) '
      minl1=0.0
      dc1=1.0
      ncont1=22
	   label='mr'      
	  End
	4:Begin
	   xsec1=surftva
		type1='!17RADIOSONDE !C!DVIRTUAL TEMPERATURE (K) '
	   label='tv'		
	  End	
	5:Begin
	   xsec1=surfrha
		type1='!17RADIOSONDE !C!DRelative Humidity (%) '
      minl1=10
      dc1=5
      ncont1=18
      label='rh'	
	  End	    
endcase
endif

;contour levels

if (field1 eq 1)  then begin
;	minl1 = min(xsec1(recst:recend, zmax:zmin))
	minl1 = 282
   dc1 = 1.0
   ncont1=30
endif	

if (field1 eq 2)  then begin
;	minl1 = min(xsec1(recst:recend, zmax:zmin))
	minl1 = 292
   dc1 = 1.0
   ncont1=30
endif	

col = reverse(indgen(ncont1)*255/ncont1)
levels1 = findgen(ncont1)*dc1 + float(minl1)




;load color table
loadct, 33
tvlct,r,b,g,/get
r(cin)=255
g(cin)=255
b(cin)=255
tvlct,r,b,g


for j=0,recno-2 do begin
	dift = time(j+1)-time(j)
	print,dift
	if dift gt .3 then begin
		xsec1(j:j+1,*)=-500
	endif
endfor


read, prompt = 'Display on the screen (0) or into the Postscript file "idl.ps" (1): ', kind

;display the contour
IF (kind NE 0) and (m eq 1) THEN BEGIN
    set_plot,'ps'
    read, prompt = 'Landscape (1) or Portrait (2): ', kind
    IF kind EQ 1 THEN BEGIN
    	  device,/landscape,/color,ysize=20,xoffset=1
    ENDIF ELSE BEGIN
        device,/portrait,/color,ysize=15
    endelse    
ENDIF

case m of 
	1: begin
	     case1=xsec1
	     recno1=recno
	     level1=levels1
	     ninterx1=ninterx
	     h1=h
	     time1=time
	     recst1=recst
	     nintery1=nintery
	     xtname1=xtname
	     ytname1=ytname
	     altl1=altl
	     sdate1=sdate
	     al1=al
	     string1=type1
	   end
	2: begin
	     case2=xsec1
	     recno2=recno
	     level2=levels1
	     ninterx2=ninterx
	     h2=h
	     time2=time
	     recst2=recst
	     nintery2=nintery
	     xtname2=xtname
	     ytname2=ytname
	     altl2=altl
	     sdate2=sdate
	     al2=al
	     string2=type1
	   end
endcase	   	   	      

;if m eq 2 then begin
;  recst = 2
;  dummy=xtname(4:ninterx-1)
;  xtname=fltarr(9)
;  xtname=dummy
;  dummy2=h(4:ninterx-1)
;  h=fltarr(9)
;  h=dummy2
;  ninterx=9
;endif 


contour, xsec1(recst:recend, *), levels = levels1, /follow, /fill,xticks=ninterx-1, $
   xcharsize = 1.0, xtickv = h, xrange = [time(recst), time(recend)], $
   xstyle = 1, yticks = nintery-1, xtickname=xtname, $
   ystyle = 1, ycharsize = 1.0, ytickname = ytname, ytickv = altl, $
   yrange = [altl(0), altl(nintery-1)], charsize=1.2,$
   background = cin, color = 0, xtitle = 'Time (GMT)', $
   title = type1 + '  CAMEX III  ANDROS ISLAND, BAHAMAS  ' + sdate(0),  $
   xgridstyle = 1, xticklen = 1.0, ygridstyle = 1, yticklen = 1.0, thick = 2.0, $
   ytitle = 'Altitude (m)', font = 1, time(recst:recend), al

contour, xsec1(recst:recend,  *), levels = levels1, /follow,/overplot, $
   time(recst:recend), al, thick = 2.0, c_colors=0
   
if m eq 3 then begin
for i = 0, recno2-1 DO $
    vline2, time(i), /data, color = 10, /noerase, linestyle = 5    
endif

endfor

	                 	
tvlct,r,g,b,/get
giffile=tvrd()
write_gif, '/home3/waynef/public_html/camexiii/archive/quicklooks/'+sdate+'lr2.gif', giffile, r, g, b
  


IF kind GT 0 THEN device,/close
  
		
end