;----------------------------------------------------
; read_sim.ncl
; Used to read data file sim.dat created by bal.f90
; In particular, reads 6 data variables: uc,vc,wc,pc,tc,and rc
; Christina Bonfanti 26 June 2015
;----------------------------------------------------

;====================================================
; load libraries and read data file
;====================================================

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
   
  nlyr = 40			; number of layers
  nlon = 52			; number of longitude points
  nlat = 42			; number of latitude points
  k_lyr = 0			; dummy layer used for the data loading
						; conditional tests
  npa = nlon*nlat

  print("")
  print("nlon = " + nlon)
  print("nlat = " + nlat)
  print("nlyr = " + nlyr)
  print("npa = " + npa)

; ======================================================================
; Read in u, v, w, p, t, and r fields from binary data file.
; ======================================================================

; ======================================================================
; Read the file into one 1-D array of length (6*nlyr+1)*npa.  npa is 
; the number of elements on each horizontal layer of each array.  There
; are 5 arrays with nlyr horizontal layers and 1 array (w) with nlyr+1 
; horizontal layers, so there are a total of 5*nlyr+(nlyr+1) = 6*nlyr+1
; horizontal layers.  
; ======================================================================

  file_name = "nimlim_int.dat"
  simDat = fbinrecread(file_name, 0, (6*nlyr+1)*npa, "float")
;  printVarSummary(simDat)

  nsize = nlyr*npa
  u_1D = simDat(0*nsize:1*nsize-1)
  v_1D = simDat(1*nsize:2*nsize-1)
  w_1D = simDat(2*nsize:3*nsize-1+npa)	; w has one more vertical layer than the other variables.
  r_1D = simDat(3*nsize+npa:4*nsize-1+npa)
  t_1D = simDat(4*nsize+npa:5*nsize-1+npa)
  p_1D = simDat(5*nsize+npa:6*nsize-1+npa)

;print("Here's p  " + p_1D(0:2*npa))	; debugging
;print("Here's t  " + t_1D(0:2*npa))	; debugging

;********************************************
; conversions
;********************************************
  
  t_1D = t_1D - 273.15				; convert temp K to C

; ======================================================================
; Change arrays from 1-D to 3-D.  Also, name the dimensions of each ar-
; ray.  This allows the dimensions to be easilty reordered.
; ======================================================================

  lat_pts = ispan(20,61,1)
  lat_pts@units = "degrees_east"
  lon_pts = ispan(0,51,1)
  lon_pts@units = "degress_north"
  alt_pts = ispan(0,19500,500)
  alt_pts@units = "m"
  altw_pts = ispan(0,20000,500)
  altw_pts@units = "m"
 
  u := onedtond(u_1D, (/nlat,nlon,nlyr/))	; works in previous
  u!0 = "lat"
  u!1 = "lon"
  u!2 = "alt"
  u&lon = lon_pts
  u&lat = lat_pts
  u&alt = alt_pts
  u&lon@units = "degrees_east"			; units to appear on axes
  u&lat@units = "degrees_north"
  u&alt@units = "m"


  v := onedtond(v_1D, (/nlat,nlon,nlyr/))
  copy_VarCoords(u, v)

  w := onedtond(w_1D, (/nlat,nlon,nlyr+1/))
  w!0 = "lat"
  w!1 = "lon"
  w!2 = "alt"
  w&lon = lon_pts
  w&lat = lat_pts
  w&alt = altw_pts
  w&lon@units = "degrees_east"
  w&lat@units = "degrees_north"
  w&alt@units = "m"

  r := onedtond(r_1D, (/nlat,nlon,nlyr/))
  copy_VarCoords(u, r)

  t := onedtond(t_1D, (/nlat,nlon,nlyr/))
  copy_VarCoords(u, t)

  p := onedtond(p_1D, (/nlat,nlon,nlyr/))
  copy_VarCoords(u, p)

;print(p(:,:,0))

; ======================================================================
; Reorder the array dimensions so that the first index is longitude, 
; the second index is latitude, and the third index is altitude.
; ======================================================================

;  u_3D := u(lon|:,lat|:,alt|:)
  uflip_3D := u(alt|:,lat|:,lon|:)
;  v_3D := v(lon|:,lat|:,alt|:)
  vflip_3D := v(alt|:,lat|:,lon|:)
;  w_3D := w(lon|:,lat|:,alt|:)
  wflip_3D := w(alt|:,lat|:,lon|:)
;  r_3D := r(lon|:,lat|:,alt|:)
  rflip_3D := r(alt|:,lat|:,lon|:)
;  t_3D := t(lon|:,lat|:,alt|:)
  tflip_3D := t(alt|:,lat|:,lon|:)
;  p_3D := p(lon|:,lat|:,alt|:)
  pflip_3D := p(alt|:,lat|:,lon|:)

  u_3D := u(lat|:,lon|:,alt|:)
  v_3D := v(lat|:,lon|:,alt|:)
  w_3D := w(lat|:,lon|:,alt|:)
  r_3D := r(lat|:,lon|:,alt|:)
  t_3D := t(lat|:,lon|:,alt|:)
  p_3D := p(lat|:,lon|:,alt|:)

;***********************************************
; create plot page for sim
;***********************************************

  wks3b = gsn_open_wks("pdf","w_cross_plot_40N")

;***********************************************
; create a contour plot for w cross section
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "W (m/s) cross-section at 40N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

;  wcontour = gsn_csm_contour(wks3b,wflip_3D(:,15,:),cnres)	; plot 35N
  wcontour = gsn_csm_contour(wks3b,wflip_3D(:,20,:),cnres)	; plot 40N
;  wcontour = gsn_csm_contour(wks3b,wflip_3D(:,25,:),cnres)	; plot 45N

;***********************************************
; create plot page for sim
;***********************************************

  wks4b = gsn_open_wks("pdf","p_5250_plot")

;***********************************************
; create a contour plot for pressure
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "P (pascals) at 5250m"		; title
 
  cnres@tiYAxisSide = "Left"
;  cnres@tiYAxisString = " Latitude"
;  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

 
;  pcontour = gsn_csm_contour(wks4b,p_3D(:,:,0),cnres)		; plot surface
  pcontour = gsn_csm_contour(wks4b,p_3D(:,:,11),cnres)	; plot 5250m
;  pcontour = gsn_csm_contour(wks4b,p_3D(:,:,3),cnres)	; plot 1250m
;  pcontour = gsn_csm_contour(wks4b,p_3D(:,:,24),cnres)	; plot 11750m

;***********************************************
; create plot page for sim
;***********************************************

  wks5b = gsn_open_wks("pdf","p_cross_plot_40N")

;***********************************************
; create a contour plot for pressure cross section
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "P (pascals) cross-section at 40N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
 
;  pcontour = gsn_csm_contour(wks5b,pflip_3D(:,15,:),cnres)	; plot 35N
  pcontour = gsn_csm_contour(wks5b,pflip_3D(:,20,:),cnres)	; plot 40N
;  pcontour = gsn_csm_contour(wks5b,pflip_3D(:,25,:),cnres)	; plot 45N

;***********************************************
; create a plot page for sim
;***********************************************

  wks6b = gsn_open_wks("pdf","t_5250_plot")

;***********************************************
; create a contour plot for temp
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "Temperature (C) at 5250m"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Latitude "
  cnres@tiXAxisString = " Longitude "

  cnres@tmXBMode = "Explicit"
  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
;  tcontour = gsn_csm_contour(wks6b,t_3D(:,:,0),cnres)	; plot surface
  tcontour = gsn_csm_contour(wks6b,t_3D(:,:,11),cnres)	; plot 5250m
;  tcontour = gsn_csm_contour(wks6b,t_3D(:,:,3),cnres)	; plot 1250m
;  tcontour = gsn_csm_contour(wks6b,t_3D(:,:,21),cnres)	; plot 11750m

;***********************************************
; create a plot page for sim
;***********************************************

  wks7b = gsn_open_wks("pdf","t_cross_plot_40N")

;***********************************************
; create a contour plot for temp cross section
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "Temperature (C) cross-section at 40N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "

;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
;  tcontour = gsn_csm_contour(wks7b,tflip_3D(:,15,:),cnres)	; plot 35N
  tcontour = gsn_csm_contour(wks7b,tflip_3D(:,20,:),cnres)	; plot 40N
;  tcontour = gsn_csm_contour(wks7b,tflip_3D(:,25,:),cnres)	; plot 45N


;***********************************************
; create plot page for sim
;***********************************************

  wks8b = gsn_open_wks("pdf","r_5250_plot")

;***********************************************
; create a contour plot for rho
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "Rho at 5250m"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Latitude "
  cnres@tiXAxisString = " Longitude "

  cnres@tmXBMode = "Explicit"
  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
;  rcontour = gsn_csm_contour(wks8b,r_3D(:,:,0),cnres)	; plot surface
  rcontour = gsn_csm_contour(wks8b,r_3D(:,:,11),cnres)	; plot 5250m
;  rcontour = gsn_csm_contour(wks8b,r_3D(:,:,3),cnres)	; plot 1250m
;  rcontour = gsn_csm_contour(wks8b,r_3D(:,:,24),cnres)	; plot 11750m

;***********************************************
; create plot page for sim
;***********************************************

  wks9b = gsn_open_wks("pdf","r_cross_plot_40N")

;***********************************************
; create a contour plot for rho cross section
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "Rho at 40N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "

;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
;  rcontour = gsn_csm_contour(wks9b,rflip_3D(:,15,:),cnres)	; plot 35N
  rcontour = gsn_csm_contour(wks9b,rflip_3D(:,20,:),cnres)	; plot 40N
;  rcontour = gsn_csm_contour(wks9b,rflip_3D(:,25,:),cnres)	; plot 45N


;***********************************************
; create plot page for sim
;***********************************************

  wks10b = gsn_open_wks("pdf","t_5250_map_plot")

;***********************************************
; create a contour plot on map of region at 
;   fixed height for temp
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "Temperature (C) at 5250m"		; title
  cnres@gsnAddCyclic = False 

  cnres@mpMaxLatF = max(t&lat)
  cnres@mpMinLatF = min(t&lat)
  cnres@mpMinLonF = min(t&lon)
  cnres@mpMaxLonF = max(t&lon)
  
  cnres@tiYAxisString = " "
  cnres@tiXAxisString = " "
;  cnres@gsnTickMarksOn = False

 
;  tmapcontour = gsn_csm_contour_map_ce(wks10b,t_3D(:,:,0),cnres)	; plot surface
  tmapcontour = gsn_csm_contour_map_ce(wks10b,t_3D(:,:,11),cnres); plot 5250m
;  tmapcontour = gsn_csm_contour_map_ce(wks10b,t_3D(:,:,24),cnres); plot 11750m
;  tmapcontour = gsn_csm_contour_map_ce(wks10b,t_3D(:,:,3),cnres); plot 1250m


;***********************************************
; createa plot for sim
; **********************************************

  wksc = gsn_open_wks("pdf","barb_cross_plot_45N")                  ; open pdf file

;************************************************
; create wind barb plot cross section
;************************************************

  res                         = True               ; plot mods desired

  res@vcRefMagnitudeF         = 10.                ; make vectors larger
  res@vcRefLengthF            = 0.050              ; ref vec length
  res@vcGlyphStyle            = "WindBarb"         ; select wind barbs 
  res@vcMinDistanceF          = 0.025              ; thin out windbarbs

  res@tiYAxisSide = "Left"
  res@tiYAxisString = " Latitude "
  res@tiXAxisString = " Longitude "
 
;  res@tmXBMode = "Explicit"
;  res@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  res@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

  res@tiMainString            = "Wind Barbs cross-section at 45N" ; title 

;  barb=gsn_csm_vector(wksc,uflip_3D(15,:,:),vflip_3D(15,:,:),res) ; cross 35N
;  barb=gsn_csm_vector(wksc,uflip_3D(15,:,:),vflip_3D(20,:,:),res) ; cross 40N
  barb=gsn_csm_vector(wksc,uflip_3D(15,:,:),vflip_3D(25,:,:),res) ; cross 45N

;***********************************************
; createa plot for sim
; **********************************************

  wks1c = gsn_open_wks("pdf","barb_11750_plot")                  ; open pdf file

;************************************************
; create wind barb plot
;************************************************

  res                         = True               ; plot mods desired

  res@vcRefMagnitudeF         = 10.                ; make vectors larger
  res@vcRefLengthF            = 0.050              ; ref vec length
  res@vcGlyphStyle            = "WindBarb"         ; select wind barbs 
  res@vcMinDistanceF          = 0.025              ; thin out windbarbs

  res@tiYAxisSide = "Left"
  res@tiYAxisString = " Latitude "
  res@tiXAxisString = " Longitude "
 
  res@tmXBMode = "Explicit"
  res@tmXBValues = (/0, 10, 20, 30, 40, 50/)
  res@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

  res@tiMainString            = "Wind Barbs at 11750m" ; title 

;  barb=gsn_csm_vector(wks1c,u_3D(:,:,0),v_3D(:,:,0),res)          ; plot surface
;  barb=gsn_csm_vector(wks1c,u_3D(:,:,11),v_3D(:,:,11),res)        ; plot 5250m
;  barb=gsn_csm_vector(wks1c,u_3D(:,:,3),v_3D(:,:,3),res)          ; plot 1250m
  barb=gsn_csm_vector(wks1c,u_3D(:,:,21),v_3D(:,:,24),res)        ; plot 11750m

;***********************************************
; create plot page for sim
;***********************************************

  wks2c = gsn_open_wks("pdf","w_11750_plot")

;***********************************************
; create a contour plot for w
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "W (m/s) at 11750m"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Latitude"
  cnres@tiXAxisString = " Longitude "
  
  cnres@tmXBMode = "Explicit"
  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

;  wcontour = gsn_csm_contour(wks2c,w_3D(:,:,0),cnres)	; plot surface
;  wcontour = gsn_csm_contour(wks2c,w_3D(:,:,11),cnres)	; plot 5250m
;  wcontour = gsn_csm_contour(wks2c,w_3D(:,:,3),cnres)	; plot 1250m
  wcontour = gsn_csm_contour(wks2c,w_3D(:,:,24),cnres)	; plot 11750m

;***********************************************
; create plot page for sim
;***********************************************

  wks3c = gsn_open_wks("pdf","w_cross_plot_45N")

;***********************************************
; create a contour plot for w cross section
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "W (m/s) cross-section at 45N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

;  wcontour = gsn_csm_contour(wks3c,wflip_3D(:,15,:),cnres)	; plot 35N
;  wcontour = gsn_csm_contour(wks3c,wflip_3D(:,20,:),cnres)	; plot 40N
  wcontour = gsn_csm_contour(wks3c,wflip_3D(:,25,:),cnres)	; plot 45N

;***********************************************
; create plot page for sim
;***********************************************

  wks4c = gsn_open_wks("pdf","p_11750_plot")

;***********************************************
; create a contour plot for pressure
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "P (pascals) at 11750m"		; title
 
  cnres@tiYAxisSide = "Left"
;  cnres@tiYAxisString = " Latitude"
;  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)

 
;  pcontour = gsn_csm_contour(wks4c,p_3D(:,:,0),cnres)		; plot surface
;  pcontour = gsn_csm_contour(wks4c,p_3D(:,:,11),cnres)	; plot 5250m
;  pcontour = gsn_csm_contour(wks4c,p_3D(:,:,3),cnres)	; plot 1250m
  pcontour = gsn_csm_contour(wks4c,p_3D(:,:,24),cnres)	; plot 11750m

;***********************************************
; create plot page for sim
;***********************************************

  wks5c = gsn_open_wks("pdf","p_cross_plot_45N")

;***********************************************
; create a contour plot for pressure
;***********************************************

  cnres = True					; plot contour
  cnres@gsnFrame = False
  cnres@cnFillOn = False
  cnres@cnLinesOn = True
  cnres@tiMainString = "P (pascals) cross-section at 45N"		; title
 
  cnres@tiYAxisSide = "Left"
  cnres@tiYAxisString = " Height (m)"
  cnres@tiXAxisString = " Longitude "
  
;  cnres@tmXBMode = "Explicit"
;  cnres@tmXBValues = (/0, 10, 20, 30, 40, 50/)
;  cnres@tmXBLabels = (/"0", "10E", "20E", "30E", "40E", "50E"/)
 
 
;  pcontour = gsn_csm_contour(wks5c,pflip_3D(:,15,:),cnres)	; plot 35N
;  pcontour = gsn_csm_contour(wks5c,pflip_3D(:,20,:),cnres)	; plot 40N
  pcontour = gsn_csm_contour(wks5c,pflip_3D(:,25,:),cnres)	; plot 45N

end