;pro batch_rtinput,args
;------------------------------------------------------------------------------
;Title:		batch_rtinput.pro
;Purpose:	Run the rtinput program which creates input files for
;			the Streamer and BugsRad RT codes from the retrieved
;			microphysical and sounding data.
;
;Note:		Requires 1-minute radar microphysics data
;
;Needed:	glas_idltemplate.dat, interpolated_positionyear.dat,
;			albedo_shebayear_crrel.dat, fluxes_year.dat,mcclatcheyozone.dat
;			zenith.pro
;
;Author:	Matthew Shupe
;Date:		1/22/02
;Modified:	1/22/02
;-------------------------------------------------------------------------------

;TO DO:
; ***Put in default soundings if not enough exist.
; ***Put in defaults or no data flags if ancillary data does not exist.

args='1998042[789]*'

;Set up bulk parameters for running the program
minave=10.0 ;Averaging time in minutes (don't use prime numbers!!)
microdir='/data2/sheba/microphysics'
radiodir='/data1/sheba/radiosonde/rawrws'
ancildir='/data1/sheba/heatrate'
rtdir='/data1/sheba/heatrate'

;Set up common blocks
common com_sonde,rtemp,rpress,rhum,ralt,rtime
common com_micro,lc,ic,lr,ir,msk,time,height,icr,irr,icd,ird,ics,irs,lcr,lrr,lc1,lr1,lc2,lr2,ixr,ixd,ixs
common com_state,alb,lat,lon,lwds,lwus,swds,swus,o3mix,o3pr

;Get a list of files and create a list of data arguments according to "args"
cd,microdir,current=orig_dir
farg='shbmcrs1microC1*'+args+'*.cdf'
rfiles=findfile(farg,count=numf)
fargs=strmid(rfiles,19,8)

;Read in all of the ancillary data that will be needed.
cd,ancildir
restore,'glas_idltemplate.dat'
restore,'interpolated_positionyear.dat'
ptime=time
restore,'albedo_shebayear_crrel.dat'
restore,'fluxes_year.dat'
ftime=time
restore,'mcclatcheyozone.dat'

for f=0,numf-1 do begin
	print,fargs[f]
	;Get all the necessary radiosonde file names
	cd,radiodir
	rarg='shbglastxt*'+fargs[f]+'.*'
	rfile=findfile(rarg,count=rnum)
	exhr=fltarr(rnum)*0.0
	if f gt 0 then begin
		rarg='shbglastxt*'+fargs[f-1]+'.*'
		rfile0=findfile(rarg,count=rnum0)
	endif else rnum0=0
	if f lt numf-1 then begin
		rarg='shbglastxt*'+fargs[f+1]+'.*'
		rfile2=findfile(rarg,count=rnum2)
	endif else rnum2=0
	if rnum0 gt 0 then begin
		if rnum gt 0 then begin
			rfile=[rfile0[rnum0-1],rfile] 
			exhr=[-24.0,fltarr(rnum)*0.0]
		endif else begin
			rfile=[rfile0[rnum0-1]]
			exhr=[-24.0]
		endelse
		rnum=rnum+1
	endif
	if rnum2 gt 0 then begin
		if rnum gt 0 then begin
		 	rfile=[rfile,rfile2[0]]
			exhr=[exhr,24.0] 
		endif else begin 
			rfile=[rfile2[0]]
			exhr=[24.0]
		endelse
		rnum=rnum+1
	endif
	if rnum eq 0 then begin
		print,'No sounding data for ',fargs[f]
	endif

	;read in micro data, select the pertinent beam(s)
	;	make average of microphysics surrounding the sounding launch time.
	cd,microdir
	fid=ncdf_open(rfiles[f])
	ncdf_varget,fid,ncdf_varid(fid,'time'),time
	ncdf_varget,fid,ncdf_varid(fid,'mask'),mask
	ncdf_varget,fid,ncdf_varid(fid,'height'),height
	ncdf_varget,fid,ncdf_varid(fid,'idmn'),idmn
	ncdf_varget,fid,ncdf_varid(fid,'idmn_s'),idmn_s
	ncdf_varget,fid,ncdf_varid(fid,'idmn_d'),idmn_d
	ncdf_varget,fid,ncdf_varid(fid,'iwc'),iwc
	ncdf_varget,fid,ncdf_varid(fid,'iwc_s'),iwc_s
	ncdf_varget,fid,ncdf_varid(fid,'iwc_d'),iwc_d
	ncdf_varget,fid,ncdf_varid(fid,'iextcoeff_d'),iextcoeff_d
	ncdf_varget,fid,ncdf_varid(fid,'iextcoeff'),iextcoeff
	ncdf_varget,fid,ncdf_varid(fid,'iextcoeff_s'),iextcoeff_s
	ncdf_varget,fid,ncdf_varid(fid,'lre'),lre
	ncdf_varget,fid,ncdf_varid(fid,'lre_s1'),lre_s1
	ncdf_varget,fid,ncdf_varid(fid,'lre_s2'),lre_s2
	ncdf_varget,fid,ncdf_varid(fid,'lwc'),lwc
	ncdf_varget,fid,ncdf_varid(fid,'lwc_s1'),lwc_s1
	ncdf_varget,fid,ncdf_varid(fid,'lwc_s2'),lwc_s2
	ncdf_close,fid

	;Average the radar data to a "minave" minute grid
	numbeam=24.0*60.0/minave
	times=(minave/2.0/60.0)+findgen(numbeam)*minave/60.0
	nhts=n_elements(height)
	ir=fltarr(numbeam,nhts)*0.0
	ic=ir & lc=ir & lr=ir & msk=ir & icr=ic & irr=ir & icd=ir & ird=ir & irs=ir & ics=ir & lcr=ir & lrr=ir & lc1=ir & lc2=ir & lr1=ir & lr2=ir & ixd=ir & ixr=ir & ixs=ir
	for i=0,numbeam-1 do begin
		iwh=where(time ge times[i]-minave/2.0/60.0 and time lt times[i]+minave/2.0/60.0,num)
		bidmn=idmn[iwh,*]
		biwc=iwc[iwh,*]
		bidmn_s=idmn_s[iwh,*]
		biwc_s=iwc_s[iwh,*]
		bidmn_d=idmn_d[iwh,*]
		biwc_d=iwc_d[iwh,*]
		blre=lre[iwh,*]
		blwc=lwc[iwh,*]
		blre_s1=lre_s1[iwh,*]
		blwc_s1=lwc_s1[iwh,*]
		blre_s2=lre_s2[iwh,*]
		blwc_s2=lwc_s2[iwh,*]
		biextcoeff_d=iextcoeff_d[iwh,*]
		biextcoeff=iextcoeff[iwh,*]
		biextcoeff_s=iextcoeff_s[iwh,*]
		msk[i,*]=mask[fix(num/2.0),*]
		for h=0,nhts-1 do begin
			tmp=[bidmn[*,h],bidmn_s[*,h],bidmn_d[*,h]]
			iwh=where(tmp gt 0 and tmp lt 1000,num)
			if num ge minave/2.0 then ir[i,h]=mean(tmp[iwh])*1.52
			tmp=[biwc[*,h],biwc_s[*,h],biwc_d[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then ic[i,h]=mean(tmp[iwh])
			tmp=[bidmn[*,h]]
			iwh=where(tmp gt 0 and tmp lt 1000,num)
			if num ge minave/2.0 then irr[i,h]=mean(tmp[iwh])*1.52
			tmp=[biwc[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then icr[i,h]=mean(tmp[iwh])
			tmp=[bidmn_d[*,h]]
			iwh=where(tmp gt 0 and tmp lt 1000,num)
			if num ge minave/2.0 then ird[i,h]=mean(tmp[iwh])*1.52
			tmp=[biwc_d[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then icd[i,h]=mean(tmp[iwh])
			tmp=[bidmn_s[*,h]]
			iwh=where(tmp gt 0 and tmp lt 1000,num)
			if num ge minave/2.0 then irs[i,h]=mean(tmp[iwh])*1.52
			tmp=[biwc_s[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then ics[i,h]=mean(tmp[iwh])
			tmp=[blre[*,h],blre_s1[*,h],blre_s2[*,h]]
			iwh=where(tmp gt 0 and tmp lt 50,num)
			if num ge minave/2.0 then lr[i,h]=mean(tmp[iwh])
			tmp=[blwc[*,h],blwc_s1[*,h],blwc_s2[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then lc[i,h]=mean(tmp[iwh])
			tmp=[blre[*,h]]
			iwh=where(tmp gt 0 and tmp lt 50,num)
			if num ge minave/2.0 then lrr[i,h]=mean(tmp[iwh])
			tmp=[blwc[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then lcr[i,h]=mean(tmp[iwh])
			tmp=[blre_s1[*,h]]
			iwh=where(tmp gt 0 and tmp lt 50,num)
			if num ge minave/2.0 then lr1[i,h]=mean(tmp[iwh])
			tmp=[blwc_s1[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then lc1[i,h]=mean(tmp[iwh])
			tmp=[blre_s2[*,h]]
			iwh=where(tmp gt 0 and tmp lt 50,num)
			if num ge minave/2.0 then lr2[i,h]=mean(tmp[iwh])
			tmp=[blwc_s2[*,h]]
			iwh=where(tmp gt 0 and tmp lt 0.5,num)
			if num ge minave/2.0 then lc2[i,h]=mean(tmp[iwh])
			tmp=[biextcoeff[*,h]]
			iwh=where(tmp gt 0 and tmp lt 10,num)
			if num ge minave/2.0 then ixr[i,h]=mean(tmp[iwh])
			tmp=[biextcoeff_d[*,h]]
			iwh=where(tmp gt 0 and tmp lt 10,num)
			if num ge minave/2.0 then ixd[i,h]=mean(tmp[iwh])
			tmp=[biextcoeff_s[*,h]]
			iwh=where(tmp gt 0 and tmp lt 10,num)
			if num ge minave/2.0 then ixs[i,h]=mean(tmp[iwh])
		endfor
	endfor


	;read in radiosonde data, make sure data are in good shape.
	; TO DO:
	; ****put in more catches to fix sounding problems
	; ****make sure alt is good.
	; ****remove alt jumps of >1km
	cd,radiodir
	rtemp=fltarr(rnum,10000)*0-999.0
	rpress=rtemp & rhum=rtemp & ralt=rtemp
	rtime=fltarr(rnum)*0-999.0
	for i=0,rnum-1 do begin
		ts=strmid(rfile[i],25,6)
		rs=read_ascii(rfile[i],template=glas_idltemplate)
		iwh=where(rs.alt eq max(rs.alt))
		iwh=iwh[0]
		if iwh eq 0 then begin
			rtime[i]=-999
		endif else begin
			cut=0
			jwh=where(rs.alt[0:iwh] gt 20000 and rs.rh[0:iwh] gt 4,num)
			for j=0,num-1 do begin
				if jwh[num-1-j] eq n_elements(rs.alt[0:iwh])-1-j then cut=cut+1
			endfor
			rtime[i]=float(strmid(ts,0,2))+float(strmid(ts,2,2))/60.0+float(strmid(ts,4,2))/3600.+exhr[i]
			rpress[i,0:iwh-cut]=rs.press[0:iwh-cut]
			rtemp[i,0:iwh-cut]=rs.temp[0:iwh-cut]
			rhum[i,0:iwh-cut]=rs.rh[0:iwh-cut]
			ralt[i,0:iwh-cut]=rs.alt[0:iwh-cut]/1000.
		endelse
	endfor
	iwh=where(rtime ne -999,num)
	if num ne rnum then begin
		rtime=rtime[iwh]
		rtemp=rtemp[iwh,*]
		rpress=rpress[iwh,*]
		rhum=rhum[iwh,*]
		ralt=ralt[iwh,*]
	endif

	;Read in position data interpolate to the radar times.
	yr=float(strmid(fargs[f],0,4))
	mo=float(strmid(fargs[f],4,2))
	da=float(strmid(fargs[f],6,2))
	jd=julday(mo,da,yr)-julday(1,1,yr)+1
	if jd lt 280 then pptime=ptime-365.0 else pptime=ptime
	iwh=where(fix(pptime) eq jd,np)
	if np gt 0 then begin
		lat=interpol(north[iwh],(pptime[iwh]-jd)*24.0,times)
		lon=interpol(west[iwh],(pptime[iwh]-jd)*24.0,times)
	endif else print,'No Lat/Lon data for ',fargs[f]

	;Read in albedo data an interpolate to the radar times
	;  ***if albedo data does not exist then fix it at.......********
	iwh=where(fix(alb_c_time) eq jd,np)
	if np gt 0 then begin
		altm=(alb_c_time[iwh]-jd)*24.0
		al=alb_c[iwh]
		alb=interpol(al,altm,times)
		iwh=where(finite(alb) eq 1)
		alb=interpol(alb[iwh],times[iwh],times)
	endif else begin
		print,'No Alb dat for ',fargs[f]
		alb=times*0.0+0.8
	endelse

	;Read in the broadband flux data, put in array that is the same
	;	length as the radar times but only fill in the appropriate values.
	if jd lt 280 then fftime=ftime-365.0 else fftime=ftime
	lwds=fltarr(n_elements(times))*0.0-999.0
	lwus=lwds
	swds=lwds
	swus=lwds
	iwh=where(fix(ftime) eq jd,num)
	if num gt 0 then begin
		tm=(ftime[iwh]-jd)*24.0+0.5
		for i=0,num-1 do begin
			wh=closest(tm[i],times,0)
			lwds[wh]=lwd[iwh[i]]
			lwus[wh]=lwu[iwh[i]]
			swds[wh]=swd[iwh[i]]
			swus[wh]=swu[iwh[i]]
		endfor
	endif

	;Make McClatchey standard ozone profile by averaging summer and winter profiles according to date
	;June is considered Summer and December is Winter
	month=fix(strmid(fargs[f],4,2))
	ratiow=[5,4,3,2,1,0,1,2,3,4,5,6]/6.0
	ratios=[1,2,3,4,5,6,5,4,3,2,1,0]/6.0
	o3ht=ozoneheight_summer		 ;Km
	o3pr=ozonepress_summer		  ;mb
	o3mix=ozonemixing_winter*ratiow[month-1]+ozonemixing_summer*ratios[month-1]  ;kg/kg

	;call rtinput routine with all the input data
	;	TO DO:  ****make the variables into a common block so the call isn't so long.
	time=times
	cd,rtdir
	rtinput,fargs[f]

endfor
cd,orig_dir
end