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

IdT = 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)
;IdT = 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)
;IdUV = NCDF_OPEN('/aim/d/waccm/sminyear/03/wa319_2x_sminyear.cam2.h1.1951-03_MeanUV.nc', /NOWRITE)
IdOO2 = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-03_MeanOO2.nc', /NOWRITE)
;IdOO2 = NCDF_OPEN('/net/aim/d/waccm/wa319_TIPHYS/WACCM/WA319JanSMinHourly_2x.cam2.h1.1993-06_MeanOO2.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, IdT, 'lon', lon         ; longitude
ncdf_varget, IdT, 'lat', lat         ; latitude
ncdf_varget, IdT, 'lev', zplev       ; sigma pressure
ncdf_varget, IdT, 'time', time_min   ; time in minites

;; Field variables

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

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

ncdf_varget, IdT, '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

ncdf_varget, IdOO2, 'O', O1      ;atomic oxygen 

ncdf_varget, IdOO2, 'O2', O2     ;molecular oxygen

NCDF_CLOSE, IdT
NCDF_CLOSE, IdUV
NCDF_CLOSE, IdOO2

mO1 = 16.E-03
mO2 = 32.E-03
mNO2 = 28.E-03

RStar = 8.314 ; Universal gas constant(N*m/mol/K)

NO2 = 1. - O1 - O2
;
; (N*m/mol/K) / ((kg/mol) x (mol/mol)) = N*m/kg/K
;
RV = RStar / (mO1 * O1 + mO2 * O2 + mNO2 * NO2)

;RV(WHERE(RV GT 1.E35)) = 1.E+35

; calculates density at every grid pt.

dens = FLTARR(nx,ny,nz,ntime)

for i = 0,nx-1 do begin 
  for j = 0,ny-1 do begin
    for k = 0,nz-1 do begin 
;(N/m**2)/(N*m/kg/K)*K=kg/m**3    
     dens[i, j, k, *] = zplev[k] * 100. / (RV[i, j, k, *] * tn[i, j, k, *])
;     dens[i, j, k] = zplev[k] * 100. / MEAN((RV[i, j, k, *] * tn[i, j, k, *]))
    endfor
  endfor
endfor

print, 'MIN(dens) ', MIN(dens)
print, 'MAX(dens) ', MAX(dens)

ntime2=ntime/2+1
nshift=(nx+1)/2-1
;frq=findgen(ntime2)
;xwvn=-findgen(nx)+nshift

;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
;
;tspcdens = FLTARR(nx,ntime2,ny,nz)
;uspcdens = FLTARR(nx,ntime2,ny,nz)
;vspcdens = FLTARR(nx,ntime2,ny,nz)
;for j = 0,ny-1 do begin
;    for it = 0,ntime2-1 do begin
;	 for i=0,nx-1 do begin 
;	    tspcdens[i, it, j, *] = tspc[i, it, j, *] / tn[i, j, *, it] * SQRT(dens[i, j, *, it])
;	    uspcdens[i, it, j, *] = uspc[i, it, j, *] * SQRT(dens[i, j, *, it])
;	    vspcdens[i, it, j, *] = vspc[i, it, j, *] * SQRT(dens[i, j, *, it])
;	 endfor
;    endfor
;endfor
;
;plot_waccm_tuv_wnfrq1to12_dens_prof, zplev, tspcdens, uspcdens, vspcdens ;uses zonalfreq0 output
;plot_waccm_tuv_wnfrq1to12_dens_prof_jun, zplev, tspcdens, uspcdens, vspcdens ;uses zonalfreq0 output
;plot_waccm_tuv_wnfrq1to12w_dens_prof, zplev, tspcdens, uspcdens, vspcdens ;uses zonalfreq0 output
;plot_waccm_tuv_wnfrq1to12wm_dens_prof, zplev, tspcdens, uspcdens, vspcdens ;uses zonalfreq0 output
;
;  Need to convert from phase in radians from zonalfreq0 to phase in degrees
;
;tphsWACCM = tphsWACCM / !dtor
;uphsWACCM = uphsWACCM / !dtor
;vphsWACCM = vphsWACCM / !dtor

;plot_waccm_tuv_wnphs1to6wm_prof, zplev, tphs, uphs, vphs 
;
;  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,tspc,nwv,nfrq,tcompd

zonalfreq1,un,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompd

zonalfreq1,vn,ntime,nx,ny,nz,nshift,ntime2,vspc,nwv,nfrq,tcompd
;
; Normalize profiles by density
;
tspcdens = FLTARR(nx,ntime2,ny,nz)
uspcdens = FLTARR(nx,ntime2,ny,nz)
vspcdens = FLTARR(nx,ntime2,ny,nz)
for j = 0,ny-1 do begin
    for it = 0,ntime2-1 do begin
	 for i=0,nx-1 do begin 
	    tspcdens[i, it, j, *] = tspc[i, it, j, *] / tn[i, j, *, it] * SQRT(dens[i, j, *, it])
	    uspcdens[i, it, j, *] = uspc[i, it, j, *] * SQRT(dens[i, j, *, it])
	    vspcdens[i, it, j, *] = vspc[i, it, j, *] * SQRT(dens[i, j, *, it])
;
;  Weight the amplitude profiles by the top data point (minimum value in density weighted profiles)
;
            tspcdens0 = tspcdens[i, it, j, 0]
            tspcdens[i, it, j, *] = tspcdens[i, it, j, *] / tspcdens0
            uspcdens0 = uspcdens[i, it, j, 0]
            uspcdens[i, it, j, *] = uspcdens[i, it, j, *] / uspcdens0
            vspcdens0 = vspcdens[i, it, j, 0]
            vspcdens[i, it, j, *] = vspcdens[i, it, j, *] / vspcdens0

	 endfor
    endfor
endfor

plot_waccm_tuv_prof_dens, zplev, tspcdens, uspcdens, vspcdens

;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 ter-diurnal wavenumber 3
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcomptd
;
;nfrq = -4
;nwv = 4 ;migrating 4th diurnal wavenumber 4
;
;zonalfreq1,tn,ntime,nx,ny,nz,nshift,ntime2,uspc,nwv,nfrq,tcompqd
;
;tcompd_sd_td_qd = tcompd + tcompsd + tcomptd + tcompqd

;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

;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

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

;  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_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

END