path='/forbesraid/processed/SABER/DAT/'

varchar=['t','u','v','w','d']

for var=0,4 do begin ; 0,4

If 0 eq 0 then begin

am=fltarr(2,8,37,161,12)
ph=fltarr(2,8,37,161,12)

infile=path+'tira_diurnal_2002_2008.nc'
fid=ncdf_open(infile)

vid=ncdf_varid(fid,'lat')
ncdf_varget,fid,vid,lat

vid=ncdf_varid(fid,'lev')
ncdf_varget,fid,vid,lev

cmp=['e3','e2','e1','s0','w1','w2']

for ss=0,5 do begin

vid=ncdf_varid(fid,'amp_'+cmp(ss)+'_'+varchar(var))
ncdf_varget,fid,vid,amp
am(0,ss,*,*,*)=amp

vid=ncdf_varid(fid,'phase_'+cmp(ss)+'_'+varchar(var))
ncdf_varget,fid,vid,pha
ph(0,ss,*,*,*)=pha

endfor
ncdf_close,fid
endif

if 0 eq 0 then begin

infile=path+'tira_semidiurnal_2002_2008.nc'

fid=ncdf_open(infile)

cmp=['e3','e2','e1','s0','w1','w2','w3','w4']

for ss=0,7 do begin

vid=ncdf_varid(fid,'amp_'+cmp(ss)+'_'+varchar(var))
ncdf_varget,fid,vid,amp
am(1,ss,*,*,*)=amp

vid=ncdf_varid(fid,'phase_'+cmp(ss)+'_'+varchar(var))
ncdf_varget,fid,vid,pha
ph(1,ss,*,*,*)=pha

endfor

ncdf_close,fid

endif

ams=fltarr(2,12,145,37)
phs=fltarr(2,12,145,37)

;for ht=0,160 do begin
for ht=39,39 do begin ; index 39 --> 97.5km
if ht le 16 then fac=exp((lev(ht)-40)/5) else fac=1
for mon=0,11 do $
for lat=0,36 do $
for nn=0,1 do begin
AA=fltarr(145)
BB=fltarr(145)
for ss=-3,2+2*nn do begin
AA=AA+am(nn,ss+3,lat,ht,mon)*cos((nn+1)*ph(nn,ss+3,lat,ht,mon)*!pi/12-(ss-nn-1)*indgen(145)*!pi/72)
BB=BB+am(nn,ss+3,lat,ht,mon)*sin((nn+1)*ph(nn,ss+3,lat,ht,mon)*!pi/12-(ss-nn-1)*indgen(145)*!pi/72)
endfor
ams(nn,mon,*,lat)=(AA^2+BB^2)^0.5*fac
if var eq 4 then ams(nn,mon,*,lat)=ams(nn,mon,*,lat)*100 
phs(nn,mon,*,lat)=atan(BB,AA)*12/!pi/(nn+1)
endfor

endfor

dens=fltarr(25,12,145,73) 
for mon=0,11 do $
for lon=0,144 do begin
thed=fltarr(25,37)
for nn=0,1 do begin
thea=ams(nn,mon,lon,*)
thep=phs(nn,mon,lon,*)  ;
for lat=0,36 do $ 
thed(*,lat)=thed(*,lat)+thea(lat)*cos((nn+1)*(indgen(25)-thep(lat))*!pi/12)
endfor
for loc=0,24 do begin
dens(loc,mon,lon,*)=interpol(reform(thed(loc,*)),indgen(37)*5-90,indgen(73)*2.5-90)
print,max(thed(loc,*)),max(dens(loc,mon,lon,*))
endfor
endfor

for mon=0,8,8 do begin ; 0,11
openw,3,'HME_tira_97km_tidal_'+varchar(var)+'_2.5x2.5_field_'+strmid(101+mon,6,2)+'.dat'
printf,3,reform(dens(*,mon,*,*))
close,3
endfor

endfor ; var

stop
end