;=========================================================== ; IDL Program to plot correlations ;=========================================================== PRO CovPlot ;----------------------------------------------------------- ; Set-up 16 colours + 0=black, 1=white ;----------------------------------------------------------- r = [0,255,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184] g = [0,255,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0] b = [0,255,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186] FOR col = 0, 7 DO BEGIN r[col+2] = 255*col/8 g[col+2] = 255*col/8 b[col+2] = 255 r[col+10] = 255 g[col+10] = 255*(7-col)/8 b[col+10] = 255*(7-col)/8 ENDFOR TVLCT, r,g,b ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- SET_PLOT, 'ps' DEVICE, FILE='Correlations.ps', /PORTRAIT, /COLOR, XSIZE=21.0, YSIZE=29.0, $ XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED ;----------------------------------------------------------- ; Set-up colour grades array ;----------------------------------------------------------- minval = -1.0 maxval = 1.0 range = maxval-minval grades = FINDGEN(16)*range/16 + minval ;=========================================================== ;=========================================================== ; STANDARD SCHEME ;=========================================================== ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN('StdCorrs.nc') varidy = NCDF_VARID(fileid, 'latitude') varidz = NCDF_VARID(fileid, 'level') varid = NCDF_VARID(fileid, 'Corr') NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varidz, vertlevs NCDF_VARGET, fileid, varid, Data NCDF_CLOSE, fileid ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nlats = N_ELEMENTS(y_degrees) Nz = N_ELEMENTS(vertlevs) yinc = y_degrees[1] - y_degrees[0] miny = y_degrees[0] maxy = y_degrees[Nlats-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Rearrange order of data and find rms values ;----------------------------------------------------------- ; The data array is indexed by [z,y]. More convenient to use [y,z] Data1 = FLTARR(Nlats,Nz) rms = 0.0 FOR iy = 0,(Nlats-1) DO BEGIN Data1(iy,*) = Data(*,iy) FOR iz = 0,(Nz-1) DO BEGIN rms = rms + Data1(iy,iz) * Data1(iy,iz) ENDFOR ENDFOR rms = SQRT(rms/(Nlats*Nz)) ;----------------------------------------------------------- ; Greenwich meridian plot ;----------------------------------------------------------- xinc = 1.18*(y_degrees[1] - y_degrees[0]) yinc = vertlevs[1] - vertlevs[0] minx = y_degrees[0] maxx = y_degrees[Nlats-1] midx = (maxx+minx)/2.0 miny = vertlevs[0] maxy = vertlevs[Nz-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- Px1 = 0.25 Px2 = Px1 + 0.6 Py1 = 0.7 Py2 = 0.95 Lx = Px2-Px1 Ly = Py2-Py1 Title = 'Standard scheme correlations: !C min =' + $ STRING(MIN(Data1(*,*))) + ', max =' + $ STRING(MAX(Data1(*,*))) + ', rms =' + $ STRING(rms) !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[miny,minx,maxy,maxx], TITLE=Title, /CYLINDRICAL, $ /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR iy = 0,(Nlats-1) DO BEGIN FOR iz = 0,(Nz-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(Data1(iy,iz) GE grades)) ;calculate the limits of the element to plot xmin = y_degrees[iy] - xinc/2.0 xmax = y_degrees[iy] + xinc/2.0 ymin = vertlevs[iz] - yinc/2.0 ymax = vertlevs[iz] + yinc/2.0 myxcoords = [xmin,xmax,xmax,xmin,xmin] myycoords = [ymin,ymin,ymax,ymax,ymin] ;correct limits if they go over the bounds pts = WHERE(myycoords GT maxy,count) IF (count GT 0) THEN myycoords(pts) = maxy pts = WHERE(myycoords LT miny,count) IF (count GT 0) THEN myycoords(pts) = miny pts = WHERE(myxcoords GT maxx,count) IF (count GT 0) THEN myxcoords(pts) = maxx pts = WHERE(myxcoords LT minx,count) IF (count GT 0) THEN myxcoords(pts) = minx POLYFILL, myxcoords, myycoords, COLOR=boxcolour+2 ENDFOR ENDFOR ;----------------------------------------------------------- ; Add zero contour ;----------------------------------------------------------- zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 CONTOUR, Data1(*,*), y_degrees, vertlevs, LEVELS=zerolevel, $ xstyle=1, ystyle=1, nlevels=1, $ /OVERPLOT, C_COLORS=linecolours, C_LINESTYLE=linestyle ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE=' ' AXIS, COLOR=0, XAXIS=1, TICKLEN=-0.015, XRANGE=[minx,maxx], XTICKNAME=blank AXIS, COLOR=0, YAXIS=0, TICKLEN=-0.015, YRANGE=[miny,maxy], YTITLE=' ' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- ;COLBAR, [Px2+0.02,Py1,Px2+0.04,Py2], grades, 'Term 1', /UP, /DOWN, /TEXTPOS ;=========================================================== ;=========================================================== ; PV SCHEME (USING ANTI-PV) ;=========================================================== ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN('PV1Corrs.nc') varidy = NCDF_VARID(fileid, 'latitude') varidz = NCDF_VARID(fileid, 'level') varid = NCDF_VARID(fileid, 'Corr') NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varidz, vertlevs NCDF_VARGET, fileid, varid, Data NCDF_CLOSE, fileid ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nlats = N_ELEMENTS(y_degrees) Nz = N_ELEMENTS(vertlevs) yinc = y_degrees[1] - y_degrees[0] miny = y_degrees[0] maxy = y_degrees[Nlats-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Rearrange order of data and find rms values ;----------------------------------------------------------- ; The data array is indexed by [z,y]. More convenient to use [y,z] Data1 = FLTARR(Nlats,Nz) rms = 0.0 FOR iy = 0,(Nlats-1) DO BEGIN Data1(iy,*) = Data(*,iy) FOR iz = 0,(Nz-1) DO BEGIN rms = rms + Data1(iy,iz) * Data1(iy,iz) ENDFOR ENDFOR rms = SQRT(rms/(Nlats*Nz)) ;----------------------------------------------------------- ; Greenwich meridian plot ;----------------------------------------------------------- xinc = 1.18*(y_degrees[1] - y_degrees[0]) yinc = vertlevs[1] - vertlevs[0] minx = y_degrees[0] maxx = y_degrees[Nlats-1] midx = (maxx+minx)/2.0 miny = vertlevs[0] maxy = vertlevs[Nz-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- Px1 = 0.25 Px2 = Px1 + 0.6 Py1 = 0.4 Py2 = 0.65 Lx = Px2-Px1 Ly = Py2-Py1 Title = 'PV scheme (with antiPV) correlations: !C min =' + $ STRING(MIN(Data1(*,*))) + ', max =' + $ STRING(MAX(Data1(*,*))) + ', rms =' + $ STRING(rms) !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[miny,minx,maxy,maxx], TITLE=Title, /CYLINDRICAL, $ /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR iy = 0,(Nlats-1) DO BEGIN FOR iz = 0,(Nz-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(Data1(iy,iz) GE grades)) ;calculate the limits of the element to plot xmin = y_degrees[iy] - xinc/2.0 xmax = y_degrees[iy] + xinc/2.0 ymin = vertlevs[iz] - yinc/2.0 ymax = vertlevs[iz] + yinc/2.0 myxcoords = [xmin,xmax,xmax,xmin,xmin] myycoords = [ymin,ymin,ymax,ymax,ymin] ;correct limits if they go over the bounds pts = WHERE(myycoords GT maxy,count) IF (count GT 0) THEN myycoords(pts) = maxy pts = WHERE(myycoords LT miny,count) IF (count GT 0) THEN myycoords(pts) = miny pts = WHERE(myxcoords GT maxx,count) IF (count GT 0) THEN myxcoords(pts) = maxx pts = WHERE(myxcoords LT minx,count) IF (count GT 0) THEN myxcoords(pts) = minx POLYFILL, myxcoords, myycoords, COLOR=boxcolour+2 ENDFOR ENDFOR ;----------------------------------------------------------- ; Add zero contour ;----------------------------------------------------------- zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 CONTOUR, Data1(*,*), y_degrees, vertlevs, LEVELS=zerolevel, $ xstyle=1, ystyle=1, nlevels=1, $ /OVERPLOT, C_COLORS=linecolours, C_LINESTYLE=linestyle ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE=' ' AXIS, COLOR=0, XAXIS=1, TICKLEN=-0.015, XRANGE=[minx,maxx], XTICKNAME=blank AXIS, COLOR=0, YAXIS=0, TICKLEN=-0.015, YRANGE=[miny,maxy], YTITLE=' ' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- ;COLBAR, [Px2+0.02,Py1,Px2+0.04,Py2], grades, 'Term 1', /UP, /DOWN, /TEXTPOS ;=========================================================== ;=========================================================== ; PV SCHEME (NOT USING ANTI-PV - DIRECT INVERSION) ;=========================================================== ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN('PV2Corrs.nc') varidy = NCDF_VARID(fileid, 'latitude') varidz = NCDF_VARID(fileid, 'level') varid = NCDF_VARID(fileid, 'Corr') NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varidz, vertlevs NCDF_VARGET, fileid, varid, Data NCDF_CLOSE, fileid ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nlats = N_ELEMENTS(y_degrees) Nz = N_ELEMENTS(vertlevs) yinc = y_degrees[1] - y_degrees[0] miny = y_degrees[0] maxy = y_degrees[Nlats-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Rearrange order of data and find rms values ;----------------------------------------------------------- ; The data array is indexed by [z,y]. More convenient to use [y,z] Data1 = FLTARR(Nlats,Nz) rms = 0.0 FOR iy = 0,(Nlats-1) DO BEGIN Data1(iy,*) = Data(*,iy) FOR iz = 0,(Nz-1) DO BEGIN rms = rms + Data1(iy,iz) * Data1(iy,iz) ENDFOR ENDFOR rms = SQRT(rms/(Nlats*Nz)) ;----------------------------------------------------------- ; Greenwich meridian plot ;----------------------------------------------------------- xinc = 1.18*(y_degrees[1] - y_degrees[0]) yinc = vertlevs[1] - vertlevs[0] minx = y_degrees[0] maxx = y_degrees[Nlats-1] midx = (maxx+minx)/2.0 miny = vertlevs[0] maxy = vertlevs[Nz-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- Px1 = 0.25 Px2 = Px1 + 0.6 Py1 = 0.1 Py2 = 0.35 Lx = Px2-Px1 Ly = Py2-Py1 Title = 'PV scheme (without antiPV) correlations: !C min =' + $ STRING(MIN(Data1(*,*))) + ', max =' + $ STRING(MAX(Data1(*,*))) + ', rms =' + $ STRING(rms) !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[miny,minx,maxy,maxx], TITLE=Title, /CYLINDRICAL, $ /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR iy = 0,(Nlats-1) DO BEGIN FOR iz = 0,(Nz-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(Data1(iy,iz) GE grades)) ;calculate the limits of the element to plot xmin = y_degrees[iy] - xinc/2.0 xmax = y_degrees[iy] + xinc/2.0 ymin = vertlevs[iz] - yinc/2.0 ymax = vertlevs[iz] + yinc/2.0 myxcoords = [xmin,xmax,xmax,xmin,xmin] myycoords = [ymin,ymin,ymax,ymax,ymin] ;correct limits if they go over the bounds pts = WHERE(myycoords GT maxy,count) IF (count GT 0) THEN myycoords(pts) = maxy pts = WHERE(myycoords LT miny,count) IF (count GT 0) THEN myycoords(pts) = miny pts = WHERE(myxcoords GT maxx,count) IF (count GT 0) THEN myxcoords(pts) = maxx pts = WHERE(myxcoords LT minx,count) IF (count GT 0) THEN myxcoords(pts) = minx POLYFILL, myxcoords, myycoords, COLOR=boxcolour+2 ENDFOR ENDFOR ;----------------------------------------------------------- ; Add zero contour ;----------------------------------------------------------- zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 CONTOUR, Data1(*,*), y_degrees, vertlevs, LEVELS=zerolevel, $ xstyle=1, ystyle=1, nlevels=1, $ /OVERPLOT, C_COLORS=linecolours, C_LINESTYLE=linestyle ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE=' ' AXIS, COLOR=0, XAXIS=1, TICKLEN=-0.015, XRANGE=[minx,maxx], XTICKNAME=blank AXIS, COLOR=0, YAXIS=0, TICKLEN=-0.015, YRANGE=[miny,maxy], YTITLE=' ' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- ;COLBAR, [Px2+0.02,Py1,Px2+0.04,Py2], grades, 'Term 1', /UP, /DOWN, /TEXTPOS DEVICE, /CLOSE END