;
; Three data files are used in this script. Two contain atmospheric data,
; and the third contains global elevation data. Both data fields will be
; overlaid on a filled contour of the elevation data.
;
a = addfile("data/94082500_nP.cdf","r")
b = addfile("data/94082500_nH.cdf","r")
e = addfile("data/topo.nc","r")

;
; Sort pressure fields
;
p = a->Psl/100.0
frtime = a->frtime
tmp = new(dimsizes(p(0,:,:)),float)
ind = 0
do i = 0, dimsizes(frtime) - 2
        ind = i
        do j = i + 1, dimsizes(frtime) - 2
                if(frtime(j) .lt. frtime(ind))
                        ind = j
                end if
        end do
        if(ind.ne.i)
                tmp = p(ind,:,:)
                p(ind,:,:) = p(i,:,:)
                p(i,:,:) = tmp
                tmp2 = frtime(ind)
                frtime(ind) = frtime(i)
                frtime(i) = tmp2
        end if
end do
;
; Determine extents of data.
;
lat = a->lat
if(lat(0) .lt. lat(dimsizes(lat)-1)) then
	ymin = lat(0)
	ymax = lat(dimsizes(lat)-1)
else
	ymax= lat(0)
	ymin = lat(dimsizes(lat)-1)
end if

lon = a->lon
if(lon(0) .lt. lon(dimsizes(lon)-1)) then
	xmin = lon(0) + 360
	xmax = lon(dimsizes(lon)-1) + 360
else
	xmax= lon(0) + 360
	xmin = lon(dimsizes(lon)-1) + 360
end if
delete(lat)
delete(lon)
;
; Determine starting indexes. This is unnecessary when coordinate
; subscripting is available.
;
s_ind_lon =0
e_ind_lon =0
s_ind_lat = 0
e_ind_lat = 0
lat = e->lat
if(lat(0) .lt. lat(dimsizes(lat)-1)) then
	do i = 0,dimsizes(lat)-1,1
		if(lat(i) .ge. ymin) then
			s_ind_lat = i
			break
		end if
	end do
	do i = dimsizes(lat)-1,0,1
		if(lat(i) .le. ymax) then
			e_ind_lat = i
			break
		end if
	end do
else
	do i = 0,dimsizes(lat)-1,1
		if(lat(i) .le. ymax) then
			s_ind_lat = i
			break
		end if
	end do
	do i = dimsizes(lat)-1,0,1
		if(lat(i) .ge. ymin) then
			e_ind_lat = i
			break
		end if
	end do
end if
lon = e->lon
if(lon(0) .lt. lon(dimsizes(lon)-1)) then
	do i = 0,dimsizes(lon)-1,1
		if(lon(i) .ge. xmin) then
			s_ind_lon = i
			break
		end if
	end do
	do i = dimsizes(lon)-1,0,1
		if(lon(i) .le. xmax) then
			e_ind_lon = i
			break
		end if
	end do
else
	do i = 0,dimsizes(lon)-1,1
		if(lon(i) .le. xmax) then
			s_ind_lon = i
			break
		end if
	end do
	do i = dimsizes(lon)-1,0,1
		if(lon(i) .ge. xmin) then
			e_ind_lon = i
			break
		end if
	end do
end if
delete(lon)
delete(lat)
;
; Define initial scalarFields.
;
field1 = create "field1" scalarFieldLayerClass noparent
        "sfDataArray" : p(0,:,:) 
        "sfMissingValueV" : p@_FillValue
        "sfXCStartV" : a->lon(0) + 360
        "sfXCEndV" : a->lon(dimsizes(a->lon)- 1) + 360
        "sfYCStartV" : a->lat(0)
        "sfYCEndV" : a->lat(dimsizes(a->lat) - 1)
end create

;
; Read elevation data into variables in memory and create scalar
; field for elevation data.
;
lon = e->lon(s_ind_lon:e_ind_lon)
lat = e->lat(s_ind_lat:e_ind_lat)
ele = e->elev(s_ind_lat:e_ind_lat,s_ind_lon:e_ind_lon)
field2 = create "field2" scalarFieldLayerClass noparent
	"sfDataArray" : ele
	"sfXCStartV" : lon(0) 
        "sfXCEndV" : lon(dimsizes(lon)-1)
        "sfYCStartV" : lat(0)
        "sfYCEndV" : lat(dimsizes(lat)-1)
end create
;
; Remove values that are no longer needed.
;
delete(lon)
delete(ele)
delete(lat)

;
; Read in first timestep of geo-potential height and create scalar field.
;
Z = b->Z(:,3,:,:)
field3 = create "field3" scalarFieldLayerClass noparent
	"sfDataArray" : Z(0,:,:)
	"sfMissingValueV" : Z@_FillValue
        "sfXCStartV" : b->lon(0) + 360
        "sfXCEndV" : b->lon(dimsizes(b->lon)- 1) + 360
        "sfYCStartV" : b->lat(0)
        "sfYCEndV" : b->lat(dimsizes(b->lat) - 1)
end create

;
; Set up color map.
;
cmap = new((/25,3/),float)
cmap(0,:) = (/1.0,1.0,1.0/)
cmap(1,:) = (/0.0,0.0,1.0/)
cmap(6,:) = (/0.0,1.0,1.0/)
cmap(7,:) = (/0.0,1.0,0.0/)
cmap(20,:) = (/.66,.25,.25/)
;
; Interpolate color values from index 1 to index 6
; for ocean areas.
;
deltr = (cmap(6,0) - cmap(1,0))/6.0
deltg = (cmap(6,1) - cmap(1,1))/6.0
deltb = (cmap(6,2) - cmap(1,2))/6.0
do i = 2, 6
	cmap(i,0) = cmap(1,0) + i * deltr
	cmap(i,1) = cmap(1,1) + i * deltg
	cmap(i,2) = cmap(1,2) + i * deltb
end do
;
; Interpolate colors from green to brown for land areas.
;
deltr = (cmap(20,0) - cmap(7,0))/13.0
deltg = (cmap(20,1) - cmap(7,1))/13.0
deltb = (cmap(20,2) - cmap(7,2))/13.0
do i = 0,13
	cmap(i+7,0) = cmap(7,0) + i * deltr
	cmap(i+7,1) = cmap(7,1) + i * deltg
	cmap(i+7,2) = cmap(7,2) + i * deltb
end do

;
; Prepare additional useful colors.
;
cmap(21,:) = (/0.6,0.3,0.6/)
cmap(22,:) = (/0.0,1.0,0.0/)
cmap(23,:) = (/1.0,0.0,1.0/)
cmap(24,:) = (/.5,.5,.5/)

;
; Open NCGM workstation and assign non-default color map to it.
;
wks = create "wks" ncgmWorkstationLayerClass noparent 
"wkColorMap" : cmap
"wkMetaName" : "mapcontour03.ncgm"
end create

;
; Remove cmap array; it is no longer needed.
;
delete(cmap)

;
; Create the map plot that will serve as base of overlay.
;
map = create "map" mapPlotLayerClass wks
	"vpXF": .05
	"vpYF": .93
	"vpWidthF": .65
	"vpHeightF": .65
	"mpCenterLatF": 0.0
	"mpCenterLonF": 0.0
	"mpProjection": "MERCATOR"
	"mpMapLimitMode" : "LATLON"
	"mpMinLatF" : ymin
	"mpMaxLatF" : ymax
	"mpMinLonF" : xmin
	"mpMaxLonF" : xmax
	"mpContinentLineColor": 0
	"mpUSStateLineColor" : 0
	"mpCountryLineColor" : 0
	"mpGridLineColor" : 1
	"mpLimbLineColor" :1
	"mpPerimLineColor" :1
end create

;
; Create the first contour field overlay on the above map;
; draw the elevation data with filled contours.
;
con2 = create "con2" contourLayerClass wks
	"vpUseSegments" : "True"
	"cnScalarFieldData" : field2
	"tfOverlayPlotBase": "True"
	"trYMinF":  ymin
	"trYMaxF":  ymax
	"trXMinF":  xmin
	"trXMaxF":  xmax
	"cnLevelSelectionMode" : "EXPLICIT"
	"cnMinLevelValF" : -1250.0
        "cnMaxLevelValF" : 3250.0
;        "cnLevelSpacingF" : 250
;        "cnMaxLevelCount" : 35
	"cnBelowMinLevelColor" : 1
	"cnAboveMaxLevelColor" : 1
	"cnBelowMinLevelFillPattern" : 0
	"cnAboveMaxLevelFillPattern" : 0
	"cnLevels" : (/-1250,-1000.0,-750.0,-500.0,-250.0,0.0,250.0,500.0,750.0,1000.0,1250.0,1500.0,1750.0,2000.0,2250.0,2500.0,2750.0,3000.0,3250.0/)
	"cnFillColors" : (/2,3,4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,20/)
	"cnMonoFillPattern" : "True"
	"cnMonoLineColor" : "True"
	"cnMonoLevelFlag" : "True"
	"cnLevelFlags" : "NOLINE"
	"cnLineColors" : 1
        "cnFillPatterns" : 0
	"cnFillOn" : "True"
        "tiMainString" : "Forecast time " + frtime(0) + " hrs"
	"cnLowLabelsOn" : "off"
	"cnHighLabelsOn" : "off"
	"cnLineLabelsOn" : "off"
	"cnInfoLabelOn" : "off"
	"ovDisplayLabelBar" : "ALWAYS"
	"lbTitleString" : "Elevation"
	"lbTitleFont" : "times-roman"
	"lbLabelFont" : "times-roman"
end create
con3 = create "con3" contourLayerClass wks
	"cnScalarFieldData" : field3
	"tfOverlayPlotBase": "True"
	"cnLevelSelectionMode" : 1
	"cnMinLevelValF" : 5300.0
        "cnMaxLevelValF" : 6000.0
        "cnLevelSpacingF" : 50
        "cnMaxLevelCount" : 35
	"trYMinF": ymin
	"trYMaxF": ymax
	"trXMinF": xmin
	"trXMaxF": xmax
        "cnMonoLineColor" : "True"
	"cnLineColors" : 25
	"cnMonoLineThickness" : "True"
	"cnLineThicknesses" : 3.0
	"cnMonoLineDashPattern" : "True"
	"cnLineDashPatterns" : 0
	"cnInfoLabelOn" : "off"
;	"cnLowLabelsOn" : "off"
;	"cnHighLabelsOn" : "off"
;	"cnLineLabelsOn" : "off"
end create
con = create "con" contourLayerClass wks
	"cnScalarFieldData" : field1
	"tfOverlayPlotBase": "True"
	"cnLevelSelectionMode" : 1
	"cnMinLevelValF" : 970.0
        "cnMaxLevelValF" : 1030.0
        "cnLevelSpacingF" : 2.5
        "cnMaxLevelCount" : 35
	"trYMinF": ymin
	"trYMaxF": ymax
	"trXMinF": xmin
	"trXMaxF": xmax
        "cnMonoLineColor" : "True"
	"cnLineColors" : 22
	"cnMonoLineThickness" : "True"
	"cnLineThicknesses" : 3.0
	"cnMonoLineDashPattern" : "True"
	"cnLineDashPatterns" : 5
	"cnInfoLabelOn" : "off"
;	"cnLowLabelsOn" : "off"
;	"cnHighLabelsOn" : "off"
;	"cnLineLabelsOn" : "off"
end create
leg = create "leg" legendLayerClass wks
	"vpXF" : .05
	"vpYF" : .25
	"vpWidthF" : .65
	"vpHeightF" : .2
	"lgItemCount" : 2
	"lgLabelStrings" : (/"Surface Pressure", "500mb Height"/)
	"lgItemStrings" : (/"", ""/)
	"lgAutoManage" : "False"
	"lgItemColors" : (/22,25/)
	"lgItemThicknesses" : (/3.0,3.0/)
	"lgItemTypes" : (/0,0/)
	"lgItemIndexes" : (/5,0/)
	"lgDrawTitle" : 0
	"lgLabelFontHeightF" : .02
	"lgLabelFont" : "times-roman"
end create
overlay(map,con2)
overlay(map,con)
overlay(map,con3)
draw(map)
draw(leg)
frame(wks)
do i = 1, dimsizes(b->frtime) - 1,1
        setvalues field1
                "sfDataArray" : p(i,:,:)
        end setvalues
        setvalues field3
                "sfDataArray" : Z(i,:,:)
        end setvalues
        setvalues con2
                "tiMainString" : "Forecast time " + frtime(i) + " hrs"
        end setvalues
	draw(leg)
        draw(map)
        frame(wks)
end do
; 
; Free variables used in this script.
;
delete(p)
delete(Z)
delete(tmp)
delete(tmp2)
;delete(wks)
;delete(field2)
;delete(frtime)
;delete(a)
;delete(p)
;delete(i)
;delete(j)
;delete(ind)
;delete(lat)
;delete(lon)
;delete(xmax)
;delete(ymax)
;delete(xmin)
;delete(ymin)