;*******************************************
; lcnative_1.ncl
;*******************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
begin
;*******************************************
; open file and read in data
;*******************************************
  f    = addfile ("MMOUT_DOMAIN1.nc", "r")
;  P = f->pslv
  P = f->t2m
  x = f->lon
  y = f->lat
  time = f->time
;
; these are the values that are interactively input to the 'invproj' command documented below
;
  print(min(x) + " " + min(y))
  print(max(x) + " " + max(y))
;********************************************
;********************************************
; get some parameters
;********************************************
  nx = dimsizes(x)
  ny = dimsizes(y)
  x2d = conform_dims((/ny,nx/),x,1)
  y2d = conform_dims((/ny,nx/),y,0)
  outstr = sprintf("%9.3f",x2d) + " " + sprintf("%9.3f",y2d)
  asciiwrite("xy.txt",outstr)

  nt = dimsizes(P(:,0,0))
;********************************************
; create plot
;********************************************
  wks  = gsn_open_wks ("png", "lcnative")          ; open workstation
  gsn_define_colormap (wks,"ncl_default")         ; choose color map
  
  res                        = True               ; plot mods desired
  res@cnFillOn               = True               ; color fill  
  res@cnLinesOn              = False              ; no contour lines
  res@cnLineLabelsOn         = False              ; no contour labels
  res@gsnSpreadColors        = True               ; use total colormap
  res@gsnSpreadColorStart    = 4 
  res@gsnSpreadColorEnd      = -1
  res@cnInfoLabelOn          = False              ; no contour info label
  res@cnMaxLevelCount = 200
  res@cnFillMode = "rasterfill"
  res@gsnMaximize = True
  res@lbBoxLinesOn = False
  
  res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
  
  res@tiMainString           = "Native Lambert Conformal Grid"
  res@tiMainFontHeightF      = 0.020              ; smaller title
  res@gsnAddCyclic           = False              ; regional data
  res@mpGridAndLimbOn    = True
  res@mpDataBaseVersion = "mediumres"
  res@mpOutlineSpecifiers = "mexico : states"
  res@mpLimitMode       = "Corners"            ; choose range of map
;
; !!!!! any plot of data that is on a native grid, must use the "corners"
; method of zooming in on map.
;
; corner points derived using the PROJ4 tool with the following parameters 
; invproj +proj=lcc  +lon_0=-124.17 +lat_1=0 +lat_2=0.01 +y_0=-27000 -f "%.4f"  +units=km +ellips=sphere 
; The min and max metric values output above serve as the interactive input to this command.
; Note that the whole 2d coordinate arrays could be generated using the invproj command as well,
; but this script uses NCL to generate the arrays after the first plot
;
  res@mpLeftCornerLatF = 6.7422
  res@mpLeftCornerLonF = -124.0487 
  res@mpRightCornerLatF  = 39.3387 
  res@mpRightCornerLonF  = -76.0217

; The following 4 pieces of information are REQUIRED to properly display
; data on a native lambert conformal grid.  This data should be specified
; somewhere in the model itself.
  
  res@mpProjection        = "LambertConformal"
  res@mpLambertParallel1F = 0.0
  res@mpLambertParallel2F = 0.01    ; apparently the LC projection cannot be tangent to the equator or cross the equator
  res@mpLambertMeridianF  = -124.17
  
; usually, when data is placed onto a map, it is TRANSFORMED to the specified
; projection. Since this model is already on a native lambert conformal grid,
; we want to turn OFF the tranformation.
  
  res@tfDoNDCOverlay = True
  
  
;  do i = 0,nt-1
  do i = 16,16
    res@gsnRightString = time(i) + " " + time@units
    plot = gsn_csm_contour_map(wks,P(i,:,:),res)     ; Draw contours over a map.
  end do

;
; once the data has been correctly plotted in its native projection, the plot can be sampled at the grid point locations (which are linearly
; spaced in the native projection) to generate 2d lat/lon coordinate arrays.
;
; first get the edges of the viewport in the NDC coordinate space
;
  getvalues plot
    "vpXF" : xn
    "vpYF" : yn
    "vpWidthF" : wn
    "vpHeightF" : hn
  end getvalues
;  print("x: " + xn + " y: " + yn + " w: " + wn + " h: " + hn)
; 
; generate the grid point locations in x and y 
;
  ndcx = fspan(xn, xn + wn,nx)
  ndcy = fspan(yn-hn,yn,ny)
;
; nudge the end points just a bit to ensure they are inside the boundary
;
  eps = 5e-7
  ndcx(0) = ndcx(0) + eps
  ndcx(nx-1) = ndcx(nx -1) - eps
  ndcy(0) = ndcy(0) + eps
  ndcy(ny-1) = ndcy(ny -1) - eps
;
  print( "calculating lat/lon coordinates")
  tndcy = new(nx, float)
  outlat = new((/ny,nx/), float)
  outlon = new((/ny,nx/), float)
;
; ndctodata generates the data coordinates (lat/lon) from the NDC locations
;
  do i = 0, ny-1
     tndcy = ndcy(i)
     ndctodata(plot,ndcx,tndcy,outlon(i,:),outlat(i,:))
  end do

;
; write out the 2d lat/lon cooridates in a NetCDF file
;

  system("rm -f latlon.nc")
  out = addfile("latlon.nc","c")
  outlon@long_name = "longitude"
  outlon@units = "degrees_east"
  outlat@long_name = "latitude"
  outlat@units = "degrees_north" 
  out->lon = outlon
  out->lat = outlat
;
; redraw the data in an orthographic projection to demonstrate the correctness of the generated lat/lon arrays
; 
  res@mpProjection		      = "orthographic"
  res@mpCenterLonF                    = -100
  res@mpCenterLatF                    = 20
  res@mpLimitMode		      = "Corners"
  res@mpLeftCornerLatF                = 2
  res@mpLeftCornerLonF                = -125
  res@mpRightCornerLatF               = 40
  res@mpRightCornerLonF               = -65
  res@tfDoNDCOverlay		      = False
  res@trGridType                      = "triangularmesh"

  res@sfXArray = outlon
  res@sfYArray = outlat
  res@tiMainString           = "Reprojected to orthographic using generated coordinates"

  res@gsnStringFontHeightF = 0.015
  plot = gsn_csm_contour_map(wks,P(16,:,:),res)     ; Draw contours over a map.

end