PRO Streamfunction


;READ-IN DATA
;------------
fid = NCDF_OPEN('Laplace_polecorrect/psi.nc')
varid = NCDF_VARID(fid,'psi')
NCDF_VARGET, fid, varid, psi1
NCDF_CLOSE, fid

fid = NCDF_OPEN('Screened_polecorrect/Psi_rec.nc')
varid = NCDF_VARID(fid,'Psi_rec')
NCDF_VARGET, fid, varid, psi2
NCDF_CLOSE, fid



psi = (psi2 - psi1) / psi1


;SET THE LEVEL TO PLOT
;---------------------
plot_level = 11-1

;MULTIPLE PLOTS ON A PAGE
;------------------------
!P.MULTI=[0,2,2]


;OPEN POSTSCRIPT FILE
;--------------------
;SET_PLOT, 'x'
SET_PLOT, 'ps'
DEVICE, FILE='blockfill.ps', /COLOR

!P.COLOR=1   
!P.THICK=3.0
!P.FONT=1
!X.THICK=3.0
!Y.THICK=3.0
!P.POSITION = [0.1, 0.1, 0.61, 0.81]

;SELECT THE COLOUR TABLE (16 colours, plus white & black)
;--------------------------------------------------------
r = [255,0,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184]
g = [255,0,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0]
b = [255,0,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186]
TVLCT, r,g,b
mycols=INDGEN(16)+2

;WRAP AROUND GREENWICH MERIDIAN
;------------------------------
;psi1 = FLTARR(97,72)
;psi1(0:95,*) = psi(0:95,*)
;psi1(96,*)   = psi(0,*)

;SET THE MAP PROJECTION
;----------------------
MAP_SET, 0., 0., LIMIT=[-90.,0.,90.,360.], /CYLINDRICAL

;SET THE AXES
;------------
lats   = -88.75 + FINDGEN(72) * 2.5
longs  = 1.875 + FINDGEN(96) * 3.75
levels = 1. + FINDGEN(38)

;SET THE PLOTTING LEVELS
;-----------------------
min = MIN(psi(*,*,plot_level))
PRINT, 'Min value = ', min
max = MAX(psi(*,*,plot_level))
PRINT, 'Max value = ', max
range = max-min

;For the colour bar
colour_levels = (FINDGEN(16) + 0.5) * range / 16.0 + min



xinc=3.75
yinc=2.5

;FILL CELLS WITH COLOUR
;----------------------
FOR ix=0,95 DO BEGIN
  FOR iy=0,71 DO BEGIN
    ;Calculate the colour of the box 
    mycol=FIX((psi(ix,iy,plot_level)-min)*16/range)
    ;calculate the limits of the bounding box
    xmin=ix*xinc-xinc/2.0
    xmax=ix*xinc+xinc/2.0
    ymin=iy*yinc-yinc/2.0-90.0
    ymax=iy*yinc+yinc/2.0-90.0
    myxcoords=[xmin,xmax,xmax,xmin,xmin]
    myycoords=[ymin,ymin,ymax,ymax,ymin]
    ;correct limits to be +/- 90 degrees
    pts=WHERE(myycoords GT 90.0,count)
    IF (count GT 0) THEN myycoords(pts)=90.0
    pts=WHERE(myycoords LT -90.0,count)
    IF (count GT 0) THEN myycoords(pts)=-90.0
    POLYFILL, myxcoords,myycoords, COLOR=mycol+2
  ENDFOR
ENDFOR


; add the continents
MAP_CONTINENTS, /NOERASE

;Add the axes 
myticklen=-0.025
myticks=['180W', '150W', '120W','90W','60W','30W','GM','30E','60E',$
          '90E', '120E', '150E', '180E']
n=N_ELEMENTS(myticks)	  
mytickvals=FINDGEN(n)/(n-1)
AXIS, -180,-90, XAXIS=0, XTITLE='Longitude', TICKLEN=myticklen,$
         XTICKS=n-1, XTICKNAME=myticks, XTICKV=mytickvals
myticks(*)=' '
AXIS, -180,90, XAXIS=1, XTITLE=' ', TICKLEN=myticklen,$
         XTICKS=n-1, XTICKNAME=myticks, XTICKV=mytickvals
      
myticks=['90S','60S', '30S','EQ',$
         '30N','60N','90N']
n=N_ELEMENTS(myticks)	  
mytickvals=FINDGEN(n)/(n-1)
AXIS, -180,-90, YAXIS=0, YTITLE='Latitude', TICKLEN=myticklen/2.5,$
         YTICKS=n-1, YTICKNAME=myticks, YTICKV=mytickvals 
myticks(*)=' '      
AXIS, 180,-90, YAXIS=1, TICKLEN=myticklen/2.5,$
         YTICKS=n-1, YTICKNAME=myticks, YTICKV=mytickvals

;add the colour bar
COLBAR, [0.67,0.16,0.73,0.78], colour_levels, 'Streamfunction', /UP, $
        /DOWN, /TEXTPOS

;Close the postscript file
DEVICE, /CLOSE




;PLOT THE DATA
;-------------
END