mirror of
https://github.com/KevinMidboe/linguist.git
synced 2025-10-29 17:50:22 +00:00
116 lines
4.5 KiB
Plaintext
116 lines
4.5 KiB
Plaintext
;*****************************************************
|
|
; mcsst_1.ncl
|
|
;
|
|
; Concepts illustrated:
|
|
; - Plotting NAVO MCSST data
|
|
; - Using fbindirread to read in fortran binary data
|
|
; - Converting "byte" data to "float"
|
|
; - Adding meta data (attributes and coordinates) to a variable
|
|
; - Adding gray to an existing color map
|
|
; - Spanning all but the last two colors in a color map for contour fill
|
|
; - Drawing raster contours
|
|
;
|
|
;*****************************************************
|
|
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
|
|
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
|
|
;***************************************
|
|
; type of data available on file
|
|
;***************************************
|
|
; ipar=0 Weekly Binned Sea Surface Temperature
|
|
; ipar=1 Number of Points in Bin
|
|
; ipar=2 Weekly Binned Sea Surface Temperature Anomaly
|
|
; ipar=3 Interpolated Sea Surface Temperature
|
|
; ipar=4 Interpolated Sea Surface Temperature Anomaly
|
|
;***************************************
|
|
begin
|
|
ipar = 3
|
|
fname = "2001311d18N16.dat"
|
|
tmp = fbindirread(fname,ipar,(/1024,2048/),"byte")
|
|
;***************************************
|
|
; convert to float and then change to true SST
|
|
;***************************************
|
|
xslope = 0.15
|
|
if(ipar.eq.4.or.ipar.eq.2)then ; anom has different intercept
|
|
yint = -20.0
|
|
end if
|
|
if(ipar.eq.3.or.ipar.eq.0)then
|
|
yint = -3.0
|
|
end if
|
|
sst = new((/1024,2048/),"float") ; create float var
|
|
sst = tmp*xslope+yint ; convert to float
|
|
delete(tmp) ; delete unecessary array
|
|
;***************************************
|
|
; assign missing values. The original missing value was zero, but since it was
|
|
; not assigned in NCL, it was not recognized. The new missing values are
|
|
; listed below. These will be changed later.
|
|
;***************************************
|
|
if(ipar.eq.4)then
|
|
sst@_FillValue = -20
|
|
end if
|
|
if(ipar.eq.3.or.ipar.eq.0)then
|
|
sst@_FillValue = -3
|
|
end if
|
|
;***************************************
|
|
; create coordinate variables
|
|
;***************************************
|
|
nlat = 1024
|
|
dy = 180./nlat
|
|
lat = (90. -(ispan(0,1023,1)*dy))-dy/2
|
|
lat!0 = "lat"
|
|
lat&lat = lat
|
|
lat@units = "degrees_north"
|
|
|
|
nlon = 2048
|
|
dx = 360./nlon
|
|
lon = (ispan(0,2047,1)*dx)+dx/2-180. ; note -180. added by sjm to align
|
|
lon!0 = "lon"
|
|
lon&lon = lon
|
|
lon@units = "degrees_east"
|
|
;***************************************
|
|
; fill out the netCDF data model
|
|
;***************************************
|
|
sst!0 = "lat" ; name dimensions
|
|
sst!1 = "lon" ; ditto
|
|
sst = sst(::-1,:) ; reverse lat orientation
|
|
sst@long_name = "NAVO MCSST" ; assign long_name
|
|
sst@units = "deg C" ; assign units
|
|
sst&lat = lat ; assign lat cv
|
|
sst&lon = lon ; assign lon cv
|
|
sst@_FillValue = -999. ; assign missing value
|
|
;***************************************
|
|
; get year and day from filename
|
|
;***************************************
|
|
res = True ; plot mods desired
|
|
title = stringtochar(fname) ; parse file name to get date
|
|
year = title(0:3)
|
|
jday = title(4:6)
|
|
res@gsnCenterString = year+" "+jday ; create center string
|
|
;***************************************
|
|
; create plot
|
|
;***************************************
|
|
wks = gsn_open_wks("ps","mcsst") ; open workstation (plot destination)
|
|
gsn_define_colormap(wks,"BlGrYeOrReVi200") ; choose colormap
|
|
;
|
|
; This will not be necessary in V6.1.0 and later. Named colors can
|
|
; be used without having to first add them to the color map.
|
|
;
|
|
d = NhlNewColor(wks,0.8,0.8,0.8) ; add gray to colormap
|
|
|
|
|
|
res@cnFillOn = True ; turn on color
|
|
res@gsnSpreadColors = True ; use full range of colormap
|
|
res@gsnSpreadColorStart = 2 ; start at color 2
|
|
res@gsnSpreadColorEnd = -3 ; don't use added gray
|
|
res@cnLinesOn = False ; no contour lines
|
|
res@cnFillDrawOrder = "PreDraw" ; draw contours before continents
|
|
res@gsnMaximize = True ; maximize plot
|
|
|
|
|
|
; For a grid this size, it is better to use raster mode. It will be
|
|
; significantly faster, and will not go over NCL's 16mb default plot size.
|
|
res@cnFillMode = "RasterFill" ; turn on raster mode
|
|
|
|
plot = gsn_csm_contour_map_ce(wks,sst,res) ; contour the variable
|
|
|
|
end
|