;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;load "../ncl-beta/contributed_beta.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

;
;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/shea_util.ncl"



begin
; We generate plots, but what kind do we prefer?
;  type = "x11"
;type = "png"
 type = "pdf"
; type = "ps"
; type = "ncgm"
  lats = (/  38.00,  42.50 /)
  lons = (/ -76.0, -68.75 /)


; Set some basic resources
  res = True
; res@MainTitle                   = "12 km WRF: "+model+"-"+PBL+"-"+micro
  res@MainTitle                   = "NAM"
;wks = gsn_open_wks(type,"fronto_NAM_02_08_2013")
 wks = gsn_open_wks(type,"fronto_test")

;
;
;
; http://www.ncl.ucar.edu/Document/Functions/Built-in/systemfunc.shtml
; to retrieve a list of files from another directory 
 Date = "20130208/"
 DATADir = "/data1/wrf/NAM/"
 diri = DATADir+Date
 all_files = systemfunc("ls " + diri + "nam*grb")
 n_files = dimsizes(all_files)
 print("n_files = " + n_files)
;do n=0,n_files -1
 do n=0,1,1
   print("all_files(n) = " + all_files(n))
  filename = all_files(n)
  print ("filename" + filename)
  a = addfile(filename,"r")
  if (n.eq.0)then
; loc = wrf_user_ll_to_ij(a, lons, lats, True)

; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
; x_start = loc(0,0) - 1
; x_end   = loc(0,1) - 1
; y_start = loc(1,0) - 1
; y_end   = loc(1,1) - 1
  end if
;
  nc_file = a
  MY_LVLS = nc_file->lv_ISBL2 
  print(MY_LVLS)
  N_LEVELS = dimsizes(MY_LVLS)
  MY_LAT = nc_file->gridlat_218
  MY_LON = nc_file->gridlon_218
  MY_PRES = nc_file->PRMSL_218_MSL
  print("MY_PRES" + MY_PRES(5,5))
  MY_HGT = nc_file->HGT_218_ISBL
  dims = dimsizes(MY_HGT) 
  print(dims)
  print(MY_HGT(:,20,20))
  MY_TMP = nc_file->TMP_218_ISBL
  dims = dimsizes(MY_TMP) 
  print(dims)
  print(MY_TMP(:,20,20))
  MY_RH  = nc_file->R_H_218_ISBL
  dims = dimsizes(MY_RH)
  print(dims)
  print(MY_RH(:,20,20))
; th_plane = wrf_user_intrp3d(th,p,"h",pressure,0.,False)
  pressure = 700
  do k=0,N_LEVELS-1
   if (MY_LVLS(k).eq.pressure)then
    k_700 = k
   end if
  end do
  print("k = " + k_700)
  pltres = True
  mpres = True
  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpNationalLineColor    = "Black"
  mpres@mpUSStateLineColor     = "Black"
  mpres@mpNationalLineColor    = "Black"
  mpres@mpGridLineColor        = 0
  mpres@mpLimbLineColor        = "Black"
  mpres@mpPerimLineColor       = "Black"
  mpres@mpUSStateLineThicknessF = 1.0
  mpres@mpGeophysicalLineThicknessF = 1.0
;  mpres@ZoomIn = True
;  mpres@Xstart = x_start
;  mpres@Ystart = y_start
;  mpres@Xend = x_end
;  mpres@Yend = y_end

;     z_dec@units="dam"
;     z_dec@short_name="Height"
      z_plane = MY_HGT(k_700,:,:)
      z_plane = z_plane/10.

;

;---Resources common to both plots
  res                = True
  res@gsnRightString = ""
  res@gsnLeftString  = ""
  res@gsnDraw        = False
  res@gsnFrame       = False
  res@gsnAddCyclic   = False

; Contour fill resources
  cnres_fill = res

  cnres_fill@gsnMaximize           = True

  cnres_fill@cnLevelSelectionMode  = "ExplicitLevels"
  cnres_fill@cnFillOn              = True
; cnres_fill@cnFillMode            = "RasterFill"     ; handles memory better, but can produce blocky contours

  cnres_fill@cnLinesOn             = False
  cnres_fill@cnLineLabelsOn        = False
; cnres_fill@cnLevels              = (/ 0, 1,2.,3, 4., 5., 6, 7, 8, 9, \
;                                      10,11,12,13,14,15,16,17,18,19,20/)
  cnres_fill@cnLevels              = (/ 0, 5,10.,15, 20., 25., 30, 35, 40, 45, \
                                       50,55,60,65,70,75,80,85,90,95,100/)

  cnres_fill@cnFillColors          = (/"White","White","AntiqueWhite","AntiqueWhite1","AntiqueWhite2",\
                                       "AntiqueWhite3","olivedrab","olivedrab1","olivedrab2",\
                                       "olivedrab3","SpringGreen","SpringGreen2","SpringGreen4",\
                                       "Yellow","Yellow2","Yellow4","Orange","Red","Red3","HotPink3",\
                                       "HotPink1","HotPink","Violet"/)


  cnres_fill@mpFillOn              = False
  cnres_fill@mpDataBaseVersion     = "MediumRes"
  cnres_fill@mpDataSetName         = "Earth..2"
  cnres_fill@mpOutlineBoundarySets = "GeophysicalAndUSStates"

  cnres_fill@mpMinLatF             = min(MY_LAT)-1
  cnres_fill@mpMinLonF             = min(MY_LON)-1
  cnres_fill@mpMaxLatF             = max(MY_LAT)+1
  cnres_fill@mpMaxLonF             = max(MY_LON)+1

; cnres_fill@tiMainString          = pett_fgen@short_name + " (" + pett_fgen@units + ")"
  cnres_fill@pmTickMarkDisplayMode = "Always"     ; nicer map tickmark labels
  pett_fgen = MY_RH(k_700,:,:)

  contour_fgen = gsn_csm_contour_map(wks,pett_fgen,cnres_fill)


;


    ; Plotting options for Z
; Plotting options for Z
  cnres_line                                = res
  cnres_line@cnConstFLabelPerimOn           = False
  cnres_line@cnLevelSelectionMode           = "ManualLevels"
  cnres_line@cnMinLevelValF                 = 210.
  cnres_line@cnMaxLevelValF                 = 390
  cnres_line@cnLevelSpacingF                = 3
  cnres_line@cnLineLabelAngleF              = 0.0
  cnres_line@cnLineLabelFontHeightF         = 0.01
  cnres_line@cnLineLabelBackgroundColor     = -1
  cnres_line@gsnContourLineThicknessesScale = 3
  cnres_line@cnLabelMasking                 = True
  cnres_line@cnInfoLabelOn                  = False
  cnres_line@cnLineLabelPerimOn             = False

  print("z_plane" + z_plane(:,250))
  contour_z = gsn_csm_contour(wks,z_plane,cnres_line)

; delete(cnres_line)
  overlay(contour_fgen,contour_z)
  draw(contour_fgen)   ; this will draw both the line and filled contours over a map
; draw(contour_z)   ; this will draw both the line and filled contours over a map

  frame(wks)



  end do
end