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

begin
;
;  Open netCDF files containing the grid and data from FV3
;

  fg1   = addfile("nnr_terra.nc","r")        ; lons and lats
  fg2   = addfile("nnr_aqua.nc","r")        ; lons and lats

  kt=0 ;18 utc
;  klev=3 ; land: 440, 470, 500, 550, 660, 870 - for new data
  klev=0 ; for old data


  if (klev .eq. 0)
     name="440"
  end if	

  if (klev .eq. 1)
     name="470"
  end if	

  if (klev .eq. 2)
     name="500"
  end if	


  if (klev .eq. 3)
     name="550"
  end if	

  if (klev .eq. 4)
     name="660"
  end if	

  if (klev .eq. 5)
     name="870"
  end if	


  lat  =                                 ndtooned(fg1->lat(:))    ; Pull off lat/lon variables

  lon  =                                 ndtooned(fg1->lon(:))    ; Pull off lat/lon variables

;  print(dimsizes(lon))
;  print(dimsizes(lat))

  nx=dimsizes(lon)
  ny=dimsizes(lat)

  grid_latt  = new(nx*ny,float)
  grid_lont  = new(nx*ny,float)
  
  ij = 0
  j = 0
  do while(j.le.ny-1)

     i = 0
     do while(i.le.nx-1)
     	  grid_latt(ij)=doubletofloat(lat(j))
     	  i=i+1
          ij=ij+1
     end do	

     j=j+1

  end do

  ij = 0
  j = 0
  do while(j.le.ny-1)

     i = 0
     do while(i.le.nx-1)
     	  grid_lont(ij)=doubletofloat(lon(i))
     	  i=i+1
          ij=ij+1
     end do	

     j=j+1

  end do

;  print(min(grid_lont))  
;  print(max(grid_lont))  

;  print(min(grid_latt))  
;  print(max(grid_latt))  

; Glue all 6 panels together

  aqua        =  ndtooned(fg1->tau(kt,klev,:,:))        ; Read some data; convert to 1D.

  terra       =  ndtooned(fg2->tau(kt,klev,:,:))

; aod
  aqua = where(aqua.gt.1.3,1.3,aqua)
  terra = where(terra.gt.1.3,1.3,terra)


  str = str_concat((/"nnr_aqua_terra_",name/))

  wks = gsn_open_wks("png",str)

;
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
;  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
;  gsn_define_colormap(wks,"BlWhRe")
;  gsn_define_colormap(wks,"NCV_blue_red")

;  gsn_reverse_colormap(wks)

;
; Set some plot resources.  Order doesn't matter here.
;
; Any time you don't want a plot drawn, just set the resources:
;
;  res@gsnDraw  = False
;  res@gsnFrame = False
;
; Of course, you'll have to remember to delete these resource
; settings or change them to True if you want later plots to
; be drawn.
;
  res                     = True
  res@gsnMaximize         = True         ; Maximize size of plot in frame
  res@gsnSpreadColors     = True         ; Use full colormap, but startblw
;  res@gsnSpreadColorStart = 24           ; at color index 24.

  res@sfXArray            = grid_lont        ; Required to tell NCL where to
  res@sfYArray            = grid_latt        ; overlay data on globe.

  res@cnFillOn            = True         ; Turn on contour fill
  res@cnFillMode          = "AreaFill"   ; Style of fill. You can also
                                         ; use "CellFill" and "RasterFill"
  res@cnLinesOn           = False        ; Turn off contour lines
  res@cnLineLabelsOn      = False        ; Turn off contour line labels

  res@lbLabelAutoStride   = True         ; Clean up labelbar labels.
  res@lbBoxLinesOn        = False        ; Turn of labelbar box lines.

;  res@gsnSColors  = True

;  res@tiMainString        = "FV3 AOD 550nm"    ; Title

;  contour = gsn_csm_contour(wks,temp,res)     ; Draw a contour plot.

  res@mpProjection          = "Orthographic"
;  res@mpProjection = "LambertEqualArea"
  res@mpDataBaseVersion     = "MediumRes"      ; Improve the resolution
                                               ; of the map outlines 
                                               ; Default is "LowRes".

  res@mpPerimOn             = False       ; Turn off map perimeter
  res@mpCenterLatF          =  45.         ; Rotate map.
  res@mpCenterLonF          =  -100.

  res@pmTickMarkDisplayMode = "Always"    ; Turn on map tickmarks

  map = gsn_csm_contour_map_overlay(wks,aqua,terra,res,res)	

end