;Id = NCDF_OPEN('/aim/d/waccm/smedyear/03/wa319_2x_smedyear.cam2.h1.1951-03-21-00000.nc', /NOWRITE)

;Id = NCDF_OPEN('/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03_MeanTZ3.nc', /NOWRITE)
IdWAM = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_mar.nc', /NOWRITE)
;IdWAM = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_mar.nc', /NOWRITE)
;IdWAM = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_mar.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/MTM/WAX3548GWSMnQUP50_2x.cam2.h1.1954_MarchMeanTZ3.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/MTM/wa319_2x_sminyear.cam2.h1.1951-06_MeanTZ3.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03-11-00000.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03-27-00000.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/MTM/WAX3548GWSMnQUP50_2x.cam2.h1.1954-03-11-00000.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/wa319_TIPHYS/MTM/WAX3548GWSMnQUP50_2x.cam2.h1.1954-03-27-00000.nc', /NOWRITE)

;; Coordinates

ncdf_varget, IdWAM, 'lon', lonWAM         ; longitude
ncdf_varget, IdWAM, 'lat', latWAM         ; latitude
ncdf_varget, IdWAM, 'lev', zplevWAM       ; log pressure, -ln(5.e-7/p) (p in hPa)
ncdf_varget, IdWAM, 'time', time_minWAM   ; time in minites

nxWAM = n_elements(lonWAM)
nyWAM = n_elements(latWAM)
nzWAM = n_elements(zplevWAM)
ntimeWAM = n_elements(time_minWAM)

;; Field variables

ncdf_varget, IdWAM, 'temp', tnWAM      ;neutral temperature
;ncdf_varget, IdWAM, 'u', unWAM      ;neutral zonal wind
;ncdf_varget, IdWAM, 'v', vnWAM      ;neutral zonal wind

NCDF_CLOSE, IdWAM

;ntime2=ntime/2+1
;nshift=(nx+1)/2-1
;frq=findgen(ntime2)
;xwvn=-findgen(nx)+nshift
;
;zonalfreq0,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;
;  The following are for the reconstruction and don't apply to amplitude plots made here of uspc only ucomp
;
;nfrq = -1
;nwv = 1 ;migrating diurnal wavenumber 1
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompd
;
;nfrq = -2
;nwv = 2 ;migrating semi-diurnal wavenumber 2
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompsd
;
;nfrq = -3
;nwv = 3 ;migrating semi-diurnal wavenumber 3
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcomptd
;
;nfrq = -4
;nwv = 4 ;migrating semi-diurnal wavenumber 4
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompqd
;
;tcompd_sd_td_qd = tcompd + tcompsd + tcomptd + tcompqd

ntime2WAM=ntimeWAM/2+1
nshiftWAM=(nxWAM+1)/2-1
frqWAM=findgen(ntime2WAM)
xwvnWAM=-findgen(nxWAM)+nshiftWAM
;
;  The following are for the reconstruction and don't apply to amplitude plots made here of uspc only ucomp
;
nfrq = -1
nwv = 1 ;migrating diurnal wavenumber 1

zonalfreq1,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcompdWAM
;zonalfreq1,unWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcompdWAM
;zonalfreq1,vnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcompdWAM

;nfrq = -2
;nwv = 2 ;migrating semi-diurnal wavenumber 2
;
;zonalfreq1,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcompsdWAM
;
;nfrq = -3
;nwv = 3 ;migrating semi-diurnal wavenumber 3
;
;zonalfreq1,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcomptdWAM
;
;nfrq = -4
;nwv = 4 ;migrating semi-diurnal wavenumber 4
;
;zonalfreq1,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,nwv,nfrq,tcompqdWAM
;
;tcompd_sd_td_qdWAM = tcompdWAM + tcompsdWAM + tcomptdWAM + tcompqdWAM

;separate procedure to handle the plotting.

;save, lat, lon, zplev, tcompd_sd_td_qd, file='tcomp_d_sd_td_sd_MarchMean_GWHR.sav'

;plotting, un, lat, zplev, uspc, ucomp
;plotting, tn, lat, zplev, uspc, tcomp

nwvn = 12

uwvtotalWAM = fltarr(nxWAM,nyWAM,nzWAM,ntimeWAM)
;uwvtotal = fltarr(nx,ny,nz,ntime)
;uwvtotal_4to12 = fltarr(nx,ny,nz,ntime)
;uwvtotal_2 = fltarr(nx,ny,nz,ntime)
;uwvtotal_4 = fltarr(nx,ny,nz,ntime)
;uwvtotal_5 = fltarr(nx,ny,nz,ntime)
;uwvtotal_6 = fltarr(nx,ny,nz,ntime)
;uwvtotal_456 = fltarr(nx,ny,nz,ntime)
;uwvtotal_7to12 = fltarr(nx,ny,nz,ntime)
;uwv1UT = fltarr(nx,ny,nz,nwvn+1)

;FOR it = 0, ntime-1 DO BEGIN
;FOR it = 0, ntimeWAM-1 DO BEGIN

;  field_in = fltarr(nx,ny,nz)
;  field_in = REFORM(tn[*,*,*,it])
;  
;  zonalwvnall,field_in,nx,ny,nz,nwvn,uspc,uphs,uwv

;  field_inWAM = fltarr(nxWAM,nyWAM,nzWAM)
;  field_inWAM = REFORM(tnWAM[*,*,*,it])
; 
;  zonalwvnall,field_inWAM,nxWAM,nyWAM,nzWAM,nwvn,uspcWAM,uphsWAM,uwvWAM

;  FOR iwvn=1,nwvn DO uwvtotalWAM[*,*,*,it] = uwvtotalWAM[*,*,*,it] + uwvWAM[*,*,*,iwvn]
;  FOR iwvn=1,nwvn DO uwvtotal[*,*,*,it] = uwvtotal[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=4,nwvn DO uwvtotal_4to12[*,*,*,it] = uwvtotal_4to12[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=4,4 DO uwvtotal_4[*,*,*,it] = uwv[*,*,*,iwvn]
;  FOR iwvn=2,2 DO uwvtotal_2[*,*,*,it] = uwv[*,*,*,iwvn]
;  FOR iwvn=5,5 DO uwvtotal_5[*,*,*,it] = uwvtotal_5[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=6,6 DO uwvtotal_6[*,*,*,it] = uwvtotal_6[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=4,6 DO uwvtotal_456[*,*,*,it] = uwvtotal_456[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=7,nwvn DO uwvtotal_7to12[*,*,*,it] = uwvtotal_7to12[*,*,*,it] + uwv[*,*,*,iwvn]

;  IF it EQ 0 THEN uwv1UT = uwv[*,*,*,*]
;  IF it EQ 3 THEN uwv1UT = uwv[*,*,*,*]

;ENDFOR

;plotting_waccm, tn, lat, zplev, uspc, tcompd
plotting_wam, tnWAM, latWAM, zplevWAM, uspcWAM, tcompdWAM
;plotting_wam_u, unWAM, latWAM, zplevWAM, uspcWAM, tcompdWAM
;plotting_wam_v, vnWAM, latWAM, zplevWAM, uspcWAM, tcompdWAM

;plot_wam_t_prof, tnWAM, zplevWAM, uspcWAM
;plot_wam_u_prof, unWAM, zplevWAM, uspcWAM
;plot_wam_v_prof, vnWAM, zplevWAM, uspcWAM

;plotting_recon_wam, tnWAM, latWAM, lonWAM, uspcWAM, tcompd_sd_td_qdWAM, uwvtotalWAM

;plotting_recon, tn, lat, zplev, uspc, tcompd_sd_td_qd, uwvtotal

;zonalwvn,tn,nx,ny,nz,uspc4,uphs4,uwv4

;plotting_wvn, tn, uwv1UT, lon

;plot_utave_lonave_t, uwvtotal_7to12, ntime, nx, ny, nz,lat


;plot_utave_lonave_z_t, uphs, 5, nx, ny, nz, zgm

END