;=========================================================== ; IDL Program to plot output terms from hydrostatic equation ;=========================================================== PRO PlotStuff ;----------------------------------------------------------- ; 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] TVLCT, r,g,b ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- SET_PLOT, 'ps' DEVICE, FILE='Hydrostatic.ps', /PORTRAIT, /COLOR, XSIZE=18.0, YSIZE=27.0, $ XOFFSET=2.0, YOFFSET=2.0 FOR plotdata = 1, 5 DO BEGIN CASE plotdata OF 1: BEGIN file1 = 'Tend_lin.nc' data1 = 'Tend_lin' title1 = 'Linear tendency terms' file2 = 'Tend_lin_FT.nc' &$ data2 = 'Tend_lin_FT' &$ title2 = 'Fourier transform of linear tendency terms' file3 = 'Tend_lin_FT.scale' title3 = 'Scale dependency of linear tendency terms' END 2: BEGIN file1 = 'Tend_nonlin.nc' data1 = 'Tend_nonlin' title1 = 'Non-linear tendency terms' file2 = 'Tend_nonlin_FT.nc' data2 = 'Tend_nonlin_FT' title2 = 'Fourier transform of non-linear tendency terms' file3 = 'Tend_nonlin_FT.scale' title3 = 'Scale dependency of non-linear tendency terms' END 3: BEGIN file1 = 'Adv_lin.nc' data1 = 'Adv_lin' title1 = 'Linear advection terms' file2 = 'Adv_lin_FT.nc' data2 = 'Adv_lin_FT' title2 = 'Fourier transform of linear advection terms' file3 = 'Adv_lin_FT.scale' title3 = 'Scale dependency of linear advection terms' END 4: BEGIN file1 = 'Adv_nonlin.nc' data1 = 'Adv_nonlin' title1 = 'Non-linear advection terms' file2 = 'Adv_nonlin_FT.nc' data2 = 'Adv_nonlin_FT' title2 = 'Fourier transform of non-linear advection terms' file3 = 'Adv_nonlin_FT.scale' title3 = 'Scale dependency of non-linear advection terms' END 5: BEGIN file1 = 'Residual.nc' data1 = 'Residual' title1 = 'Hydrostatic residual terms' file2 = 'Residual_FT.nc' data2 = 'Residual_FT' title2 = 'Fourier transform of hydrostatic residual terms' file3 = 'Residual_FT.scale' title3 = 'Scale dependency of hydrostatic residual terms' END ENDCASE ;----------------------------------------------------------- ; PLOT 1: REAL-SPACE DATA ;----------------------------------------------------------- PRINT, 'Dealing with ',file1, data1 ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN(File1) varidx = NCDF_VARID(fileid, 'x') varidy = NCDF_VARID(fileid, 'y') varid = NCDF_VARID(fileid, Data1) NCDF_VARGET, fileid, varidx, x_degrees NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varid, data NCDF_CLOSE, fileid x_degrees = TEMPORARY(x_degrees) - 360 ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nx = N_ELEMENTS(x_degrees) Ny = N_ELEMENTS(y_degrees) xinc = x_degrees[1] - x_degrees[0] yinc = y_degrees[1] - y_degrees[0] minx = x_degrees[0] maxx = x_degrees[Nx-1] midx = (maxx+minx)/2.0 miny = y_degrees[0] maxy = y_degrees[Nx-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up colour grades array ;----------------------------------------------------------- minval = MIN(data) maxval = MAX(data) grades = FINDGEN(16)*(maxval-minval)/16 + minval ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- ; Px, Py, etc are bounds of on-page plotting area in normalised co-ords Px1 = 0.05 Px2 = 0.8 Py1 = 0.7 Py2 = 1.0 Lx = Px2-Px1 Ly = Py2-Py1 !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[miny,minx,maxy,maxx], TITLE=Title1, /CYLINDRICAL ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR ix = 0,(Nx-1) DO BEGIN FOR iy = 0,(Ny-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(data(ix,iy) GE grades)) ;calculate the limits of the element to plot xmin = x_degrees[ix] - xinc/2.0 xmax = x_degrees[ix] + xinc/2.0 ymin = y_degrees[iy] - yinc/2.0 ymax = y_degrees[iy] + 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 contour plot ; /FILL ; C_COLORS=INDGEN(16)+2 ; /FOLLOW ; /OVERPLOT ;----------------------------------------------------------- linecolours = INDGEN(16) & linecolours = 0 ;All contours black ;linecolours = INDGEN(16)+2 ;Contours many colours CONTOUR, data, x_degrees, y_degrees, xstyle=1, ystyle=1, nlevels=16, $ /OVERPLOT, C_COLORS=linecolours ;----------------------------------------------------------- ; Add continents ;----------------------------------------------------------- ;MAP_CONTINENTS, /NOERASE ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE='degrees' 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='degrees' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- COLBAR, [0.85,Py1,0.95,Py2], grades, 'Term 1', /UP, /DOWN, /TEXTPOS ;----------------------------------------------------------- ; PLOT 2: 2d FOURIER TRANSFORM DATA ;----------------------------------------------------------- PRINT, 'Dealing with ',file2, data2 ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN(File2) varidx = NCDF_VARID(fileid, 'kx') varidy = NCDF_VARID(fileid, 'ky') varid = NCDF_VARID(fileid, Data2) NCDF_VARGET, fileid, varidx, kx NCDF_VARGET, fileid, varidy, ky NCDF_VARGET, fileid, varid, data NCDF_CLOSE, fileid ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nx = N_ELEMENTS(kx) Ny = N_ELEMENTS(ky) xinc = kx[1] - kx[0] yinc = ky[1] - ky[0] minx = kx[0] maxx = kx[Nx-1] midx = (maxx+minx)/2.0 miny = ky[0] maxy = ky[Nx-1] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up colour grades array ;----------------------------------------------------------- minval = MIN(data) maxval = MAX(data) grades = FINDGEN(17)*(maxval-minval)/17 + minval ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- ; Px, Py, etc are bounds of on-page plotting area in normalised co-ords Px1 = 0.05 Px2 = 0.8 Py1 = 0.325 Py2 = 0.625 Lx = Px2-Px1 Ly = Py2-Py1 !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[0.0,0.0,1.0,1.0], TITLE=Title2, /CYLINDRICAL, /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR ix = 0,(Nx-1) DO BEGIN FOR iy = 0,(Ny-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(data(ix,iy) GE grades)) ;calculate the limits of the element to plot xmin = (ix - 0.5)*Lx/Nx + Px1 xmax = (ix + 0.5)*Lx/Nx + Px1 ymin = (iy - 0.5)*Ly/Ny + Py1 ymax = (iy + 0.5)*Ly/Ny + Py1 myxcoords = [xmin,xmax,xmax,xmin,xmin] myycoords = [ymin,ymin,ymax,ymax,ymin] POLYFILL, myxcoords, myycoords, COLOR=boxcolour+1, /NORMAL ENDFOR ENDFOR ;----------------------------------------------------------- ; Add contour plot ; /FILL ; C_COLORS=INDGEN(16)+2 ; /FOLLOW ; /OVERPLOT ;----------------------------------------------------------- ;linecolours = INDGEN(16) & linecolours = 0 ;All contours black ;linecolours = INDGEN(16)+2 ;Contours many colours ;CONTOUR, data, kx, ky, xstyle=1, ystyle=1, nlevels=16, $ ; /OVERPLOT, C_COLORS=linecolours ;----------------------------------------------------------- ; Add continents ;----------------------------------------------------------- ;MAP_CONTINENTS, /NOERASE ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE='kx' 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='ky' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- COLBAR, [0.85,Py1,0.95,Py2], grades, 'Term 1', /UP, /DOWN, /TEXTPOS ;----------------------------------------------------------- ; PLOT 3: FOURIER TRANSFORM IN TERMS OF TOTAL WAVENUMBER ;----------------------------------------------------------- PRINT, 'Dealing with ',file3 ;----------------------------------------------------------- ; Read data from text file ;----------------------------------------------------------- tempstring = '' tempfloat1 = 0.0 tempfloat2 = 0.0 tempfloat3 = 0.0 tempint = 0 lengthscale1 = 0.0 weight1 = 0.0 lengthscale = FLTARR(49) weight = FLTARR(49) OPENR, unit, File3, /GET_LUN READF, unit, tempstring ;null first line FOR line = 0, 48 DO BEGIN READF, unit, tempfloat1, lengthscale1, weight1, tempint, tempfloat2 lengthscale(line) = lengthscale1/1000.0 ;Convert to km weight(line) = weight1 ENDFOR CLOSE, unit FREE_LUN, unit ;----------------------------------------------------------- ; Set-up the plot position ;----------------------------------------------------------- ; Px, Py, etc are bounds of on-screen plotting area in normalised co-ords Px1 = 0.05 Px2 = 0.8 Py1 = 0.05 Py2 = 0.25 Lx = Px2-Px1 Ly = Py2-Py1 !P.POSITION = [Px1, Py1, Px2, Py2] PLOT, lengthscale, weight, XSTYLE=5, YSTYLE=5, PSYM=-1, COLOR=13, TITLE=Title3, /NOERASE AXIS, XAXIS=0, XTITLE='Scale', COLOR=0 AXIS, XAXIS=1, XTICKNAME=blank, COLOR=0 AXIS, YAXIS=0, YTITLE='Weight', COLOR=0 AXIS, YAXIS=1, YTICKNAME=blank, COLOR=0 ENDFOR DEVICE, /CLOSE END