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

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

If 0 eq 0 then begin

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

infile='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='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

glon=fltarr(8)
dglon=360/8
for n=0,7 do begin
 glon(n) = n*dglon
endfor
mon=10
;lat=18
fac=1
ht=39

ams=fltarr(2,12,8,37)
phs=fltarr(2,12,8,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(8)
BB=fltarr(8)
;for ss=-3,2+2*nn do begin
ss=0
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	; nn

;endfor	; ht

dens=fltarr(25,12,8,73) 
;for mon=0,11 do $
for lon=0,7 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,'97km_tidal_'+varchar(var)+'_2.5x45_field_'+strmid(101+mon,6,2)+'.dat'
printf,3,reform(dens(*,mon,*,*))
close,3
;endfor

;endfor ; var

stop
end