;#####################################################################
;---------------------------------------------------------------------
;_________ SCRIPT OF STUFF TO GET GOING WITH "XPROTO" ________________
;---------------------------------------------------------------------
;	There are 4 basic steps:
;		1. Up-front housekeeping
;		2. Read in your data
;		3. Perform some actions based on your data
;		4. Plot and customize
;---------------------------------------------------------------------
;...Here we go!.........
;
;....Enter the following bits of NCL into the NCL Editor window:...........
;
;#####################################################################
;#####################  UP-FRONT HOUSEKEEPING  #######################
;#####################################################################
;...You can just copy-and-paste all this stuff.......
;
;..Get some useful procedures (colormap generation, map bkgd options,...):....
;	procedure map1 (wid,mapid)
;	procedure map2 (wid,mapid)
;	procedure map3 (wid,mapid)
;	procedure mapfill (mapid)
;	procedure set_mplims(mapid, lat1,lat2, lon1,lon2)
;	procedure adjust_colormap(wid)
;	procedure set_cnfill(wid, cnid)
;	procedure black_on_white (wid)
;	procedure white_on_black (wid)
;	procedure cnlabels(cnid)
;	procedure contour(appid,wid,cnid,dataid, ARRAY[*][*],lat1,lat2,lon1,lon2)
;	function min(a,b)
;	function max(a,b)
 load "/home/jps/ncarg/ncarg4.0/ncl/test/basic_fns.ncl"

;__________
;...adjust colormap to have 10 elements at start for ocean, land, etc.
;....(see "mpFillColors")....
 adjust_colormap(xprotoX)

;__________
 delete(MapPlotObj0)
 MapPlotObj0 = create "MapPlotObj0" mapPlotClass xprotoX
     "vpYF" 		: 0.9
     "vpHeightF" 	: 0.8
     "vpXF" 		: 0.1
     "vpWidthF" 	: 0.7
     "mpGeophysicalLineDashPattern" : 5
     "mpUSStateLineDashPattern"     : 2
     "mpLabelsOn"       : False
 end create
;...Set up "mpFillColors" for land, water........
 mapfill (MapPlotObj0)

;__________
 delete(ContourObj0)
 ContourObj0 = create "ContourObj0" contourPlotClass xprotoX
     "vpYF" 		: 0.85
     "vpHeightF" 	: 0.7
     "vpXF" 		: 0.12
     "vpWidthF" 	: 0.7
     "pmLabelBarDisplayMode" : "Conditional" ;
     "pmLabelBarWidthF"      : 0.07         ;
     "pmLabelBarHeightF"     : 0.55         ;
     "lbPerimOn"     : False		   ;
     "lbLabelFont"   : 21 		   ;-Eventually, put all these in a proc
     "tmXBLabelFont"   : 21                ;
     "tmYLLabelFont"   : 21               ;
     "tiXAxisOn"   : True                ;
     "tiYAxisOn"   : True                ;
     "tiXAxisFont"   : 21                ;
     "tiYAxisFont"   : 21                ;
     "tiXAxisFontHeight"   : 0.022       ;
     "tiYAxisFontHeight"   : 0.022       ;
 end create
;...Set a bunch of contour label specs....
  cnlabels(ContourObj0)
;...Pick colorfill colors from the wid colormap...
  set_cnfill(xprotoX, ContourObj0)

;__________
;...Filename label...
 delete(TextObj0)	
 TextObj0 = create "TextObj0" textItemClass xprotoX
    "txFuncCode"    : "&"
    "txPosXF"       : 0.012
    "txPosYF"       : 0.97
    "txFontHeightF" : 0.011
    "txJust"        : 1
    "txFont"        : 4
 end create
;__________
;...Variable name label...
 delete(TextObj1)	
 TextObj1 = create "TextObj1" textItemClass xprotoX
    "txFuncCode"    : "&"
    "txPosXF"       : 0.012
    "txPosYF"       : 0.95
    "txFontHeightF" : 0.011
    "txJust"        : 1
    "txFont"        : 4
 end create
;__________
;...Slice info label...
 delete(TextObj2)	
 TextObj2 = create "TextObj2" textItemClass xprotoX
    "txFuncCode"    : "&"
    "txPosXF"       : 0.012
    "txPosYF"       : 0.93
    "txFontHeightF" : 0.011
    "txJust"        : 1
    "txFont"        : 4
 end create

;__________
;...CREATE "TEMPLATE" DATA OBJECT.....
;....(- use a bogus data array to initialize with)...
 delete (dataid)
 ARRAY=new((/2,2/), float)
 dataid = create "data_1" scalarFieldClass defaultapp
        "sfDataArray" 	: ARRAY
 end create

 epsilon=.0001

;#####################################################################
;#######################  READ IN YOUR DATA  #########################
;#####################################################################
;...READ IN DATA FILE.....
 filename="/net/jps/archive/b/jps/archive/b/ljp/unpack/ECD_to_NMC/Dec85/EC.85120100.nc"
 file1=addfile(filename,"r")
 setvalues TextObj0
    "txString"      : "File : "+filename
 end setvalues
; 
;...To find out what's in your file, enter this:...
 print (file1)
;
;...EXTRACT A HORIZ. 2-D SLICE OF A SPECIFIC VARIABLE......................
 setvalues TextObj1
    "txString"      : "Variable : "+file1->T@long_name
 end setvalues
 setvalues dataid
    "sfMissingValueV"      : file1->T@missing_value
 end setvalues
;
;slice_method=1
 slice_method=2
;
 if (slice_method.eq.1) then
;  .................HERE'S THE PAINFUL PART:..............
;  ..Method #1..Preferred, since no question of data limits...
   it1=0
   k1=5
     slice="Slice : "+file1->T!1 + "="+file1->pres(k1)+", "+file1->T!0 + "="+file1->time(it1)
     iaxlab=file1->lon@long_name
     jaxlab=file1->lat@long_name
   j1=0
   j2=70
   i1=0
   i2=143
     lata=file1->lat(j1)
     latb=file1->lat(j2)
     lon1=file1->lon(i1)+epsilon
     lon2=file1->lon(i2)-epsilon
     lat1=amin(lata,latb)+epsilon			;Note: "amin" and "amax" are MY functions!
     lat2=amax(lata,latb)-epsilon
     print("minimum longitude = "+lon1)
     print("maximum longitude = "+lon2)
     print("minimum latitude = "+lat1)
     print("maximum latitude = "+lat2)
   delete(TEMP)
;   TEMP=file1->T(it1, k1, j1:j2, i1:i2)
   TEMP=file1->T(it1, k1, {lat1:lat2}, i1:i2)
 else
;
;...Method #2....
   it1=0
   pres1=500.
   pres2=500.
     slice="Slice : "+file1->T!1 + "="+pres1+", "+file1->T!0 + "="+file1->time(it1)
     iaxlab=file1->lon@long_name
     jaxlab=file1->lat@long_name
   lat1=-80.
   lat2=80.
   lon1=10.
   lon2=350.
   delete(TEMP)
   TEMP=file1->T(it1, {pres1:pres2}, {lat1:lat2}, {lon1:lon2} )
 end if

;...for either method, set the slice label....
 setvalues TextObj2
    "txString"      : slice
 end setvalues
 setvalues ContourObj0
    "tiXAxisString"  : iaxlab
    "tiYAxisString"  : jaxlab
 end setvalues


;#####################################################################
;#####################  DATA-DEPENDENT ACTIONS  ######################
;#####################################################################
;.............................................................................
;...Set the map limits based on what data was selected....
 set_mplims(MapPlotObj0, lat1,lat2, lon1,lon2)

;...do the contouring...
 contour(xproto,xprotoX,ContourObj0,dataid, TEMP,lat1,lat2,lon1,lon2)

;...to overlay:...
 overlay(MapPlotObj0,ContourObj0)		; mapid is in charge of the vp setup and coords

;...to undo overlay:....
 restore = True
 NhlRemoveOverlay(MapPlotObj0,ContourObj0,restore)

;#####################################################################
;#######################  PLOT AND CUSTOMIZE  ########################
;#####################################################################
;  The following can be done all within the GUI (xproto), but sometimes 
;    they can be handy cut-and-paste shortcuts.....

 draw(ContourObj0)

 clear(xprotoX)

 draw(MapPlotObj0)