; plot_wam_tuvamp_wn5frq5_latlev_seasonal
;
;  Gets seasonal WAM data and plots migrating diurnal wavenumber 5 latitude cross-section
;

SET_PLOT, 'PS'

;psfile = '/local/i/joemci/SW/WAM/tuvphs_WN5Frq5_latlev_seasonal_WAM.ps'
psfile = 'tuvphs_WN5Frq5_latlev_seasonal_WAM.ps'

DEVICE, /COLOR, /LAND, FILENAME = psfile
LOADCT, 39

xTitle = 'Latitude (Degrees)'
yTitle = 'Pressure (hPa)'
!P.CHARSIZE = 1.2

pTop = 1.E-9
;pBot = 1.E-1
pBot = 1000.0

nsea = 4

for isea = 0, nsea-1 do begin

  if isea eq 0 then IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_mar.nc', /NOWRITE)
  if isea eq 0 then IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_mar.nc', /NOWRITE)
  if isea eq 0 then IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_mar.nc', /NOWRITE)
  if isea eq 1 then IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_jun.nc', /NOWRITE)
  if isea eq 1 then IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_jun.nc', /NOWRITE)
  if isea eq 1 then IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_jun.nc', /NOWRITE)
  if isea eq 2 then IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_sep.nc', /NOWRITE)
  if isea eq 2 then IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_sep.nc', /NOWRITE)
  if isea eq 2 then IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_sep.nc', /NOWRITE)
  if isea eq 3 then IdWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_t_dec.nc', /NOWRITE)
  if isea eq 3 then IdUWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_u_dec.nc', /NOWRITE)
  if isea eq 3 then IdVWAM = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WAM/wam_v_dec.nc', /NOWRITE)

  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

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

  NCDF_CLOSE, IdWAM
  NCDF_CLOSE, IdUWAM
  NCDF_CLOSE, IdVWAM

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

  ntime2WAM=ntimeWAM/2+1
  nshiftWAM=(nxWAM+1)/2-1

  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
;
; WAM temperature westward wavenumber 5 and zonal frequency 5
;
  tmppWAM5 = REFORM(tphsWAM[100,5,*,0:149])
 
  if isea eq 0 then Title = ' WAM Temperature Phase Migrating 4.8 Hour Zonal Wavenumber 5 March Mean '
  if isea eq 1 then Title = ' WAM Temperature Phase Migrating 4.8 Hour Zonal Wavenumber 5 June Mean '
  if isea eq 2 then Title = ' WAM Temperature Phase Migrating 4.8 Hour Zonal Wavenumber 5 September Mean '
  if isea eq 3 then Title = ' WAM Temperature Phase Migrating 4.8 Hour Zonal Wavenumber 5 December Mean '
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10, /fill,c_colors=BYTSCL(indgen(10)),xrange=[-90.0,90.0],yrange=[pBot,pTop], xstyle=1,ystyle=1,/ylog,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3]
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10,/follow,xrange=[-90.0,90.0], yrange=[pBot,pTop],xstyle=1,ystyle=1,/ylog,/overplot,c_charsize=1.2
;
; WAM zonal wind westward wavenumber 5 and zonal frequency 5
;
  tmppWAM5 = REFORM(uphsWAM[100,5,*,0:149])
 
  if isea eq 0 then Title = ' WAM Zonal Wind Phase Migrating 4.8 Hour Wavenumber 5 March Mean '
  if isea eq 1 then Title = ' WAM Zonal Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 June Mean '
  if isea eq 2 then Title = ' WAM Zonal Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 September Mean '
  if isea eq 3 then Title = ' WAM Zonal Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 December Mean '
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10, /fill,c_colors=BYTSCL(indgen(10)),xrange=[-90.0,90.0],yrange=[pBot,pTop], xstyle=1,ystyle=1,/ylog,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3]
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10,/follow,xrange=[-90.0,90.0], yrange=[pBot,pTop],xstyle=1,ystyle=1,/ylog,/overplot,c_charsize=1.2
;
; WAM meridional wind westward wavenumber 5 and zonal frequency 5
;
  tmppWAM5 = REFORM(vphsWAM[100,5,*,0:149])
 
  if isea eq 0 then Title = ' WAM Meridional Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 March Mean '
  if isea eq 1 then Title = ' WAM Meridional Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 June Mean '
  if isea eq 2 then Title = ' WAM Meridional Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 September Mean '
  if isea eq 3 then Title = ' WAM Meridional Wind Phase Migrating 4.8 Hour Zonal Wavenumber 5 December Mean '
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10, /fill,c_colors=BYTSCL(indgen(10)),xrange=[-90.0,90.0],yrange=[pBot,pTop], xstyle=1,ystyle=1,/ylog,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3]
  CONTOUR, tmppWAM5[*,*],latWAM, zplevWAM, nlev=10,/follow,xrange=[-90.0,90.0], yrange=[pBot,pTop],xstyle=1,ystyle=1,/ylog,/overplot,c_charsize=1.2

endfor

DEVICE, /CLOSE

SET_PLOT, 'x'

END