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

Id = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-03_MeanTZ3.nc', /NOWRITE)
IdUV = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-03_MeanUVZ3.nc', /NOWRITE)
;Id = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-06_MeanTZ3.nc', /NOWRITE)
;IdUV = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-06_MeanUVZ3.nc', /NOWRITE)
;Id = NCDF_OPEN('/net/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03_MeanTZ3.nc', /NOWRITE)
;IdUV = NCDF_OPEN('/net/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03_MeanUV.nc', /NOWRITE)
;Id = NCDF_OPEN('/net/aim/d/waccm/sminyear/12/wa319_2x_sminyear.cam2.h1.1951-12_MeanTZ3.nc', /NOWRITE)
;IdUV = NCDF_OPEN('/net/aim/d/waccm/sminyear/12/wa319_2x_sminyear.cam2.h1.1951-12_MeanUVZ3.nc', /NOWRITE)
;Id = NCDF_OPEN('/net/aim/d/waccm/sminyear/06/wa319_2x_sminyear.cam2.h1.1951-06-Mean.nc', /NOWRITE)
;IdUV = NCDF_OPEN('/net/aim/d/waccm/sminyear/06/wa319_2x_sminyear.cam2.h1.1951-06-Mean.nc', /NOWRITE)
IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_mar.nc', /NOWRITE)
IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_mar.nc', /NOWRITE)
IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_mar.nc', /NOWRITE)
;IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_dec.nc', /NOWRITE)
;IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_dec.nc', /NOWRITE)
;IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_dec.nc', /NOWRITE)
;IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_jun.nc', /NOWRITE)
;IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_jun.nc', /NOWRITE)
;IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_jun.nc', /NOWRITE)
;Id = NCDF_OPEN('/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03_MeanUV.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, Id, 'lon', lon         ; longitude
ncdf_varget, Id, 'lat', lat         ; latitude
ncdf_varget, Id, 'lev', zplev       ; log pressure, -ln(5.e-7/p) (p in hPa)
ncdf_varget, Id, 'time', time_min   ; time in minutes

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 minutes

;; Field variables

ncdf_varget, Id, 'Z3', zg     ;geopotential height
zgm = zg * (1. + zg/6370000.0) /1000.

nx = n_elements(lon)
ny = n_elements(lat)
nz = n_elements(zplev)
ntime = n_elements(time_min)

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

; Dimensions lon,lat,lev,time

ncdf_varget, Id, 'T', tn      ;neutral temperature

ncdf_varget, IdUV, 'U', un      ;neutral zonal wind

ncdf_varget, IdUV, 'V', vn     ;neutral meridional wind

;ncdf_varget, Id, 'OMEGA', om     ;neutral vertical wind

; Dimensions lon,lat,lev,time

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

NCDF_CLOSE, Id
NCDF_CLOSE, IdUV

NCDF_CLOSE, IdWAM
NCDF_CLOSE, IdUWAM
NCDF_CLOSE, IdVWAM

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
;zonalfreq0,un,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;zonalfreq0,vn,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
;zonalfreq1,un,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompd
;zonalfreq1,vn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompd

ntime2WAM=ntimeWAM/2+1
nshiftWAM=(nxWAM+1)/2-1
frqWAM=findgen(ntime2WAM)
xwvnWAM=-findgen(nxWAM)+nshiftWAM
;
;  The following are for the 2d reconstruction and don't apply to amplitude plots made here of 
;  uspc only ucomp
;
;nfrq = -1
;nwv = 1 ;migrating diurnal wavenumber 1
;
;  The following is for calculation and plotting the t u, and v amplitude and phase line plots 
;  of WACCM and WAM for diurnal, semi-diurnal, and ter-diurnal migrating tides.
;
;xxx
;zonalfreq0,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;zonalfreq0,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,uphsWAM
;plot_waccm_wam_t_prof_phase, zplev, zplevWAM, uphs, uphsWAM
;plot_waccm_wam_t_prof, zplev, zplevWAM, uspc, uspcWAM
;zonalfreq0,un,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;zonalfreq0,unWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,uphsWAM
;plot_waccm_wam_u_prof_phase, zplev, zplevWAM, uphs, uphsWAM
;plot_waccm_wam_u_prof, zplev, zplevWAM, uspc, uspcWAM
;zonalfreq0,vn,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;zonalfreq0,vnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,uphsWAM
;plot_waccm_wam_v_prof_phase, zplev, zplevWAM, uphs, uphsWAM
;plot_waccm_wam_v_prof, zplev, zplevWAM, uspc, uspcWAM

;
;plot_waccm_wam_tuv_wvn,tn,tnWAM, un,unWAM,vn,vnWAM

;zonalfreq0,tn,ntime,nx,ny,nz,nshift,ntime2,tspc,tphs
;zonalfreq0,un,ntime,nx,ny,nz,nshift,ntime2,uspc,uphs
;zonalfreq0,vn,ntime,nx,ny,nz,nshift,ntime2,vspc,vphs

;plot_waccm_tuv_prof, zplev, tspc, uspc, vspc
;plot_waccm_tuv_wn5frq1to12_prof, zplev, tspc, uspc, vspc

;zonalfreq0,tnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,tspcWAM,tphsWAM
;zonalfreq0,unWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,uspcWAM,uphsWAM
;zonalfreq0,vnWAM,ntimeWAM,nxWAM,nyWAM,nzWAM,nshiftWAM,ntime2WAM,vspcWAM,vphsWAM

;plot_wam_tuv_prof, zplevWAM, tspcWAM, uspcWAM, vspcWAM
;plot_wam_tuv_wn5frq1to12_prof, zplevWAM, tspcWAM, uspcWAM, vspcWAM
;plot_waccm_wam_tuvamp_wn5frq5_latlev, lat, zplev, tspc, uspc, vspc, latWAM, zplevWAM, tspcWAM, uspcWAM, vspcWAM
;
;  Need to convert from phase in radians from zonalfreq0 to phase in degrees
;
;tphs = tphs / !dtor
;uphs = uphs / !dtor
;vphs = vphs / !dtor
;tphsWAM = tphsWAM / !dtor
;uphsWAM = uphsWAM / !dtor
;vphsWAM = vphsWAM / !dtor

;plot_waccm_wam_tuv_wnphs1to6wm_prof, zplev, tphs, uphs, vphs, zplevWAM, tphsWAM, uphsWAM, vphsWAM ;zonalfreq0
;plot_waccm_wam_tuvphs_wn5frq5_latlev, lat, zplev, tphs, uphs, vphs, latWAM, zplevWAM, tphsWAM, uphsWAM, vphsWAM

;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

;plot_waccm_wam_t_prof_eq, zplev, zplevWAM, tn, tnWAM, un, unWAM, vn, vnWAM

;plot_waccm_wam_t_prof_phase, zplev, zplevWAM, uphs, uphsWAM
;plot_waccm_wam_u_prof_phase, zplev, zplevWAM, uphs, uphsWAM
;plot_waccm_wam_v_prof_phase, zplev, zplevWAM, uphs, uphsWAM

;plot_waccm_wam_t_prof, zplev, zplevWAM, uspc, uspcWAM
;plot_waccm_wam_u_prof, zplev, zplevWAM, uspc, uspcWAM
;plot_waccm_wam_v_prof, zplev, zplevWAM, uspc, uspcWAM

;plotting_waccm, tn, lat, zplev, uspc, tcompd

;plotting_waccm_u, un, lat, zplev, uspc, tcompd

;plotting_waccm_v, vn, lat, zplev, uspc, tcompd
 
;plot_waccm_t_prof, tn, zplev, uspc
;plot_waccm_u_prof, un, zplev, uspc
;plot_waccm_v_prof, un, zplev, uspc
 
;plot_recon_waccm,tcompsd,lat,lon,zplev,time_min,nx,ny,nz,ntime

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

;plot_recon_lon, tcompd, lon, lev

;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

nwvn = 12

twv1ZAllT = FLTARR(nx,ny,ntime,nwvn+1)
uwv1ZAllT = FLTARR(nx,ny,ntime,nwvn+1)
vwv1ZAllT = FLTARR(nx,ny,ntime,nwvn+1)

twv1TAllZ = FLTARR(nx,ny,nz,nwvn+1)
uwv1TAllZ = FLTARR(nx,ny,nz,nwvn+1)
vwv1TAllZ = FLTARR(nx,ny,nz,nwvn+1)

tphs1TAllZ = FLTARR(nx,ny,nz)
uphs1TAllZ = FLTARR(nx,ny,nz)
vphs1TAllZ = FLTARR(nx,ny,nz)

ntw = 0 ; Do time 00Z for WAM
nt = 16 ; Do time 16Z for WACCM

FOR it = 0, ntime-1 DO BEGIN

  field_in = fltarr(nx,ny,nz)
  field_in = REFORM(tn[*,*,*,it])
  
  zonalwvnall,field_in,nx,ny,nz,nwvn,uspc,tphs,twv

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

  field_in = REFORM(vn[*,*,*,it])
  
  zonalwvnall,field_in,nx,ny,nz,nwvn,uspc,vphs,vwv

;  Grab the pressure level near 285km.  WACCM is top down
;
;  Find altitude index for this time
;    
  FOR iLn=0,nX-1 DO BEGIN

    FOR iLT = 0, nY-1 DO BEGIN

      zWLnLtT = zgm[iLn,iLt,*,iT]
      zInd = WHERE(zWLnLtT GE 275.0 AND zWLnLtT LE 295.0)

      IF zInd[0] EQ -1 THEN BEGIN
  
        MESSAGE, 'Cannot find altitude index in WACCM'
     
      ENDIF ELSE BEGIN
  
        iZ = zInd[0]
     
      ENDELSE

;print, 'iZ is ', iZ

      FOR iwvn=0,nwvn DO twv1ZAllT[iLn,iLt,it,iwvn] = twv[iLn,iLt,iZ,iwvn]
      FOR iwvn=0,nwvn DO uwv1ZAllT[iLn,iLt,it,iwvn] = uwv[iLn,iLt,iZ,iwvn]
      FOR iwvn=0,nwvn DO vwv1ZAllT[iLn,iLt,it,iwvn] = vwv[iLn,iLt,iZ,iwvn]

    ENDFOR

  ENDFOR

;  Grab the specifed UT 

  IF it eq nt THEN BEGIN 
    twv1TAllZ = twv
    uwv1TAllZ = uwv
    vwv1TAllZ = vwv
    
    tphs1TAllZ = tphs
    uphs1TAllZ = uphs
    vphs1TAllZ = vphs
  ENDIF

ENDFOR

twv1ZAllTWAM = FLTARR(nxWAM,nyWAM,ntimeWAM,nwvn+1)
uwv1ZAllTWAM = FLTARR(nxWAM,nyWAM,ntimeWAM,nwvn+1)
vwv1ZAllTWAM = FLTARR(nxWAM,nyWAM,ntimeWAM,nwvn+1)

twv1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM,nwvn+1)
uwv1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM,nwvn+1)
vwv1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM,nwvn+1)

tphs1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM)
uphs1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM)
vphs1TAllZWAM = FLTARR(nxWAM,nyWAM,nzWAM)

FOR it = 0, ntimeWAM-1 DO BEGIN

  field_in = fltarr(nxWAM,nyWAM,nzWAM)
  field_in = REFORM(tnWAM[*,*,*,it])
  
  zonalwvnall,field_in,nxWAM,nyWAM,nzWAM,nwvn,uspcWAM,tphsWAM,twvWAM

  field_in = REFORM(unWAM[*,*,*,it])
  
  zonalwvnall,field_in,nxWAM,nyWAM,nzWAM,nwvn,uspcWAM,uphsWAM,uwvWAM

  field_in = REFORM(vnWAM[*,*,*,it])
  
  zonalwvnall,field_in,nxWAM,nyWAM,nzWAM,nwvn,uspcWAM,vphsWAM,vwvWAM

;  Grab the pressure level near 285km. WAM is bottom up

  FOR iwvn=0,nwvn DO twv1ZAllTWAM[*,*,it,iwvn] = twvWAM[*,*,142,iwvn]
  FOR iwvn=0,nwvn DO uwv1ZAllTWAM[*,*,it,iwvn] = uwvWAM[*,*,142,iwvn]
  FOR iwvn=0,nwvn DO vwv1ZAllTWAM[*,*,it,iwvn] = vwvWAM[*,*,142,iwvn]

;  Grab the specifed UT 

  IF it eq ntw THEN BEGIN 
    twv1TAllZWAM = twvWAM
    uwv1TAllZWAM = uwvWAM
    vwv1TAllZWAM = vwvWAM
    
    tphs1TAllZWAM = tphsWAM
    uphs1TAllZWAM = uphsWAM
    vphs1TAllZWAM = vphsWAM
  ENDIF

ENDFOR

;plot_waccm_wam_wvn1to12,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zgm

;plot_waccm_wam_wvn,twv1ZAllT,uwv1ZAllT,vwv1ZAllT,lat,lon,twv1ZAllTWAM,uwv1ZAllTWAM,vwv1ZAllTWAM,latWAM,lonWAM

;plot_waccm_wam_tuv_wvn_recon_12s,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM

;plot_waccm_wam_tuv_wvn_recon,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM ;zonalwvnall

;plot_waccm_wam_tuv_wvn_recon_dec,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM

;plot_waccm_wam_tuv_wvn_recon_jun,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM

;plot_waccm_wam_tuv_wvn_recon_jun_ut16,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM

;plot_waccm16ut_wam00ut_tuv_wvn_recon_jun,twv1TAllZ,twv1TAllZWAM,uwv1TAllZ,uwv1TAllZWAM,vwv1TAllZ,vwv1TAllZWAM,zplev, zplevWAM, nz, nzWAM ;zonalwvnall

;plot_waccm_wam_recon_lon,tn,tnWAM,lon,lonWAM,twv1ZAllT,twv1ZAllTWAM

plot_waccm_wam_recon_lon_eq,zgm,tn,tnWAM,lon,lonWAM,twv1ZAllT,twv1ZAllTWAM ;zonalwvnall
;plot_waccm_wam_recon_lon_15n_jun,zgm,tn,tnWAM,lon,lonWAM,twv1ZAllT,twv1ZAllTWAM ;zonalwvnall

;plot_waccm_wam_t_prof_wvnphase, zplev, zplevWAM, tphs1TAllZ, tphs1TAllZWAM

;  FOR iwv = 0, 11 DO BEGIN
;
;    nwvn = iwv + 1
;
;  FOR iwvn=1,nwvn DO uwvtotal[*,*,*,it] = uwvtotal[*,*,*,it] + uwv[*,*,*,iwvn]
;  FOR iwvn=4,nwvn DO uwvtotal_4to12[*,*,*,it] = uwvtotal_4to12[*,*,*,it] + uwv[*,*,*,iwvn]
;  uwvAll[*,*,it,iwvn] = uwv[*,*,4,iwvn]
;  uwv_2[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_3[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_4[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_5[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_6[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_7[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_8[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_9[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_10[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_11[*,*,*,it] = uwv[*,*,*,iwvn]
;  uwv_12[*,*,*,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

END