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

Name = "12_27_00"
Date = "0000 UTC 27 December 1993"
hour=0

ncid3 = addfile("q2700.nc","r")
time3 = ncid3->time(0, 5)
pressure3 = ncid3->pressure(:)
lat3 = ncid3->latitude(:)
lon3 = ncid3->longitude(:)
PV = ncid3->Q(0,2:9,:,:)
printMinMax(PV,True)



ncid = addfile("analysis_fields.nc", "r")
time = ncid->time(:, 5) ; First number is the time, second number is the character length
SLP = ncid->pressure(1)
plotdata = ncid->hght(hour,{500},:,:)
height_500 = ncid->hght(hour,{500},:,:)
;time = ncid->time(1, 2) ; First number is the time, second number is the character length
temp1 = ncid->tmpk(hour,1:8,:,:)
pressure1 = ncid->pressure(1:8)
lat1 = ncid->latitude(:) 
long1 = ncid->longitude(:)

ncid2 = addfile("psml_fields.nc", "r")
time2 = ncid2->time(:,2)
lat2 = ncid2->latitude(:)
lon2 = ncid2->longitude(:)
mslp = ncid2->pmsl(hour,0,:,:)
none = ncid2->none(0)

;Calculate theta for each grid point. After, convert pressure levels to isentropic levels
p0=1000. ; Reference pressure level to calculate theta

print(pressure1)
print(PV&pressure)
print(temp1&pressure)

;theta = temp1(1:8)*(p0/pressure3(2:9))^0.286  ; calculate potential temperature
pressure_4d=conform(PV,pressure1,0)

theta = temp1*(p0/pressure_4d)^0.286  ; calculate potential temperature

printMinMax(temp1,True)
printMinMax(theta,True)



lvl = ispan(390,300,30)*1.  ; specify desired isentropic levels
xlvl = int2p_n (theta, PV, lvl, 1, 0) ;Change to isentropic coordinates (xlvl)



 xlvl!0    = "theta"                     ; name dimensions
 xlvl!1    = "lat1"                      
 xlvl!2    = "lon1"

 xlvl&theta = lvl
 xlvl&lat1 = lat2
 xlvl&lon1 = lon2

 xlvl@_FillValue = 9999.9

 xlvl = xlvl / 100

;  xlvl&time =  x&time                    ; assign coordinates
;  xlvl&lvl  =  lvl                       ; isentropic levels
;  xlvl&lat  =  x&lat
;  xlvl&lon  =  x&lon
;  xlvl@long_name = x@long_name           ; attributes
;  xlvl@units     = x@units


printVarSummary(xlvl)
printMinMax(xlvl,True)

resmap                        = True
resmap@gsnMaximize            = True
resmap@gsnDraw                = False
resmap@gsnFrame               = False
resmap@gsnPaperOrientation    = "Portrait"

resmap@mpProjection           = "CylindricalEquidistant"
resmap@gsnAddCyclic           = False
resmap@mpLimitMode            = "LatLon"

resmap@mpFillOn               = False
resmap@mpOutlineOn            = True

resmap@mpGeophysicalLineColor      = "Black"
resmap@mpGeophysicalLineThicknessF = 5

resmap@mpGridAndLimbOn        = "True"
resmap@mpGridAndLimbDrawOrder = "Draw"
resmap@mpFillDrawOrder        = "PreDraw"
resmap@mpGridMaskMode         = "MaskNone"
resmap@mpGridSpacingF         = 15.0
resmap@mpGridLineColor        = "Grey37"

resmap@tmXBLabelFontHeightF  = 0.015
resmap@tmYLLabelFontHeightF  = 0.015
resmap@tmXBMajorLengthF      = -0.001
resmap@cnInfoLabelOn         = False
resmap@cnLineLabelsOn        = False ; EDITED: True to False
resmap@pmTickMarkDisplayMode = "Never"

resmap@gsnContourNegLineDashPattern = 16
resmap@cnLineThicknessF       = 0.5
resmap@cnFillOn               = False
resmap@cnLevelSelectionMode = "AutomaticLevels"
;resmap@cnLevelSelectionMode   = "ManualLevels"
;resmap@cnMinLevelValF         = 4800.
;resmap@cnMaxLevelValF         = 6000.
;resmap@cnLevelSpacingF        = 60.

resmap@mpMinLatF            = 25.0
resmap@mpMaxLatF            = 65.0
resmap@mpMinLonF            = -120.0
resmap@mpMaxLonF            = -30.0
resmap@mpRelativeCenterLon = True
resmap@cnFillOn = True
resmap@cnFillLines = False
resmap@tiMainString = Date
;gsn_define_colormap(wks,"MPL_rainbow")


;Resources for MSLP stuff
  res2=False
  res2@cnFillOn = False
  res2@cnFillLines = True
  res2@cnFillColor="green"
resmap@cnLevelSelectionMode   = "ManualLevels"
  res2@cnMinLevelValF         = 950.
  res2@cnMaxLevelValF         = 1050.
  res2@cnLevelSpacingF        = 5.
  res2@gsnDraw                = False
  res2@gsnFrame               = False
  res2@cnLineLabelsOn        = True ; EDITED: True to false
  res2@cnLineThicknessF = 6.          ; line thickness
  res2@cnLineLabelsOn = True
  res2@cnInfoLabelOn         = False
  res2@cnFillOn               = False
  res2@cnLevelFlags = "LineAndLabel"
  res2@cnLineLabelDensityF      = 1.0             ; <1.0 = less, >1.0 = more
  res2@cnLineLabelBackgroundColor = -1
  res2@cnLineLabelPlacementMode = "constant"
;  res2@cnLineLabelInterval        = 1             ; default = 2

;ncid = addfile("../output/analysis_fields.nc","r")
;plotdata = ncid->hght(5,{500},:,:)

filename = "test_pv_invert"
wks = gsn_open_wks("png", Name)
gsn_define_colormap(wks,"MPL_rainbow")

;wrf_smooth_2d(height_500,2)
;wrf_smooth_2d(mslp,2)

plot1 = gsn_csm_contour_map_ce(wks,xlvl({330},:,:),resmap)  

;plot1 = gsn_csm_contour_map(wks,height_500,resmap)

;plot2 = gsn_csm_contour(wks,mslp,res2)

;overlay(plot1,plot2)

draw(plot1)
frame(wks)

format = "png"
system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+Name+".png -trim +repage "+Name+"."+format)

;draw(plot1)
;frame(wks)

;delete(plot)
;delete(wks)

end