;=========================================================== ; IDL Program to plot output terms from hydrostatic equation ;=========================================================== PRO PVPlots ;----------------------------------------------------------- ; 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 TestDir = 'Test12/' nlevs = 50 ;----------------------------------------------------------- ; Read-in some data from the stdout file ; (convergence and pressure profile information) ;----------------------------------------------------------- linestring = '' relerr = FLTARR(45000) itercount = -1 p_before = FLTARR(nlevs) & p_before(*) = 0 & done_p_before = 0 p_after = FLTARR(nlevs) & p_after(*) = 0 & done_p_after = 0 OPENR, unit, TestDir+'stdout.html', /GET_LUN WHILE (EOF(unit) EQ 0) DO BEGIN READF, unit, linestring IF (STRMID(linestring, 1, 4) EQ 'Step') THEN BEGIN itercount = itercount + 1 relerr(itercount) = STRMID(linestring, 25) ;PRINT, itercount, relerr(itercount) ENDIF IF ( (STRMID(linestring, 1, 25) EQ 'ROSS DARC PBMEAN (before)') AND (done_p_before EQ 0) ) THEN BEGIN FOR iz = 0, nlevs-1 DO BEGIN READF, unit, linestring p_before(iz) = linestring ENDFOR done_p_before = 1 ENDIF IF ( (STRMID(linestring, 1, 24) EQ 'ROSS DARC PBMEAN (after)') AND (done_p_after EQ 0) ) THEN BEGIN FOR iz = 0, nlevs-1 DO BEGIN READF, unit, linestring p_after(iz) = linestring ENDFOR done_p_after = 1 ENDIF ENDWHILE CLOSE, unit FREE_LUN, unit SET_PLOT, 'ps' ;----------------------------------------------------------- ; Loop over each plot to output ;----------------------------------------------------------- FOR plotdata = 1, 13 DO BEGIN CASE plotdata OF 1: BEGIN dothis = 1 psfile = 'PsiBal.eps' file = 'CntrlVars.nc' dataname = 'unspecified_1' title = 'Bal. Psi' long = 'longitude_1' lat = 'latitude_1' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 2: BEGIN dothis = 1 psfile = 'PV.eps' file = 'PV_Residuals.nc' dataname = 'unspecified' title = 'PV' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 3: BEGIN dothis = 1 psfile = 'PVres.eps' file = 'PV_Residuals.nc' dataname = 'unspecified_1' title = 'Residual PV' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 4: BEGIN dothis = 0 psfile = 'aPV.eps' file = 'aPV_Residuals.nc' dataname = 'unspecified' title = 'Anti-PV' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 5: BEGIN dothis = 0 psfile = 'aPVres.eps' file = 'aPV_Residuals.nc' dataname = 'unspecified_1' title = 'Residual anti-PV' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 6: BEGIN dothis = 1 psfile = 'Unbalp.eps' file = 'CntrlVars.nc' dataname = 'unspecified_2' title = 'Unbal p' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 7: BEGIN dothis = 1 psfile = 'Balp.eps' file = 'Pressures.nc' dataname = 'unspecified' title = 'Bal p' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 8: BEGIN dothis = 1 psfile = 'TotalPsi.eps' file = 'PsiGuess.nc' dataname = 'unspecified' title = 'Total and first guess Psi' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 9: BEGIN dothis = 1 psfile = 'UnbalPsi.eps' file = 'UnbalPsi.nc' dataname = 'unspecified' title = 'Unbal Psi' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 10: BEGIN dothis = 1 psfile = 'ResPsi.eps' file = ' ' dataname = ' ' title = 'Residual Psi' long = ' ' lat = ' ' lev = ' ' plusminus = 1 readdata = 0 END 11: BEGIN dothis = 1 psfile = 'Totalp.eps' file = 'CntrlVars.nc' dataname = 'unspecified' title = 'Total pressure' long = 'longitude' lat = 'latitude' lev = 'hybrid_ht' plusminus = 1 readdata = 1 END 12: BEGIN dothis = 1 psfile = 'Resp.eps' file = ' ' dataname = ' ' title = 'Residual p' long = ' ' lat = ' ' lev = ' ' plusminus = 1 readdata = 0 END 13: BEGIN dothis = 1 psfile = 'pfg.eps' file = ' ' dataname = ' ' title = 'First guess p' long = ' ' lat = ' ' lev = ' ' plusminus = 1 readdata = 0 END ENDCASE IF (dothis EQ 1) THEN BEGIN ;-------------------------- ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- DEVICE, FILE=TestDir+psfile, /PORTRAIT, /COLOR, XSIZE=34.0, YSIZE=8.0, $ XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED PRINT, 'Dealing with ',title ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- IF (readdata EQ 1) THEN BEGIN fileid = NCDF_OPEN(TestDir+file) varidx = NCDF_VARID(fileid, long) varidy = NCDF_VARID(fileid, lat) varidz = NCDF_VARID(fileid, lev) varid = NCDF_VARID(fileid, dataname) NCDF_VARGET, fileid, varidx, x_degrees NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varidz, vertlevs NCDF_VARGET, fileid, varid, data NCDF_CLOSE, fileid ENDIF ;----------------------------------------------------------- ; Decide which height levels to plot ;----------------------------------------------------------- lookatlevels = INDGEN(3) lookatlevels(0) = 1 lookatlevels(1) = 11 lookatlevels(2) = 30 ; Take away 1 due to IDL array indexing FOR levelloop = 0, 2 DO BEGIN lookatlevels(levelloop) = lookatlevels(levelloop)-1 ENDFOR ;----------------------------------------------------------- ; Find limits of domain, etc ;----------------------------------------------------------- Nx = N_ELEMENTS(x_degrees) Ny = N_ELEMENTS(y_degrees) Nz = N_ELEMENTS(vertlevs) 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[Ny-1] midy = (maxy+miny)/2.0 IF (Nz NE nlevs) THEN BEGIN PRINT, 'Nz and nlevs do not match. Please set nlevs at start of program to ',Nz STOP ENDIF ;----------------------------------------------------------- ; Do we need to store this field for later processing ; or make some calculations on the field? ;----------------------------------------------------------- CASE plotdata OF 1: BEGIN ;balanced Psi psibal = data END 6: BEGIN ;Unbalanced pressure punbal = data END 7: BEGIN ;Balanced pressure pbal = data FOR ix = 0,(Nx-1) DO BEGIN FOR iy = 0,(Ny-1) DO BEGIN FOR iz = 0,(nlevs-1) DO BEGIN pbal(ix,iy,iz) = data(ix,iy,iz) + p_after(iz) ENDFOR ENDFOR ENDFOR END 8: BEGIN ;Psi first guess psitotal = data END 9: BEGIN ;unbalanced Psi psiunbal = data END 10: BEGIN ;residual Psi data = psitotal - (psibal + psiunbal) END 11: BEGIN ;total p ptotal = data END 12: BEGIN ;residual p data = ptotal - (pbal + punbal) END 13: BEGIN ;first guess p data = ptotal - pbal END ELSE : PRINT, 'No calculations for this variable' ENDCASE ;----------------------------------------------------------- ; Set-up colour grades array ;----------------------------------------------------------- minval1= FINDGEN(3) maxval1= FINDGEN(3) FOR levelloop = 0, 2 DO BEGIN minval1(levelloop) = MIN(data(*,*,lookatlevels(levelloop))) maxval1(levelloop) = MAX(data(*,*,lookatlevels(levelloop))) ENDFOR minval = MIN(minval1) maxval = MAX(maxval1) IF (plusminus EQ 1) THEN BEGIN ;We expect the data to be roughly balanced in +/- extreme values IF (ABS(minval) GT ABS(maxval)) THEN BEGIN maxval = ABS(minval) ENDIF ELSE BEGIN minval = -ABS(maxval) ENDELSE ENDIF range = maxval-minval grades = FINDGEN(16)*range/16 + minval ;----------------------------------------------------------- ; Loop round the levels to be plotted ;----------------------------------------------------------- FOR levelloop = 0, 2 DO BEGIN PRINT, 'Level ',lookatlevels(levelloop) ;----------------------------------------------------------- ; 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.25 ;Px2 = 0.7 ;Py1 = 1 - levelloop*0.24 - 0.16 ;Py2 = 1 - levelloop*0.24 Px1 = levelloop*0.22 + 0.03 Px2 = Px1 + 0.2 Py1 = 0.1 Py2 = 0.8 Lx = Px2-Px1 Ly = Py2-Py1 !P.POSITION = [Px1, Py1, Px2, Py2] ; Develop a new title that includes other information (max, min) TitleA = Title + ' level' + STRING(lookatlevels(levelloop)+1) + $ '!C Min = ' + STRING(minval1(levelloop)) + $ '!C Max = ' + STRING(maxval1(levelloop)) MAP_SET, 0.0, 180.0, LIMIT=[miny,minx,maxy,maxx], TITLE=TitleA, /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,lookatlevels(levelloop)) 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 zero contour ;----------------------------------------------------------- zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 ;linecolours = INDGEN(16) & linecolours = 0 ;All contours black ;linecolours = INDGEN(16)+2 ;Contours many colours CONTOUR, data(*,*,lookatlevels(levelloop)), x_degrees, y_degrees, LEVELS=zerolevel, $ xstyle=1, ystyle=1, nlevels=1, $ /OVERPLOT, C_COLORS=linecolours, C_LINESTYLE=linestyle ;----------------------------------------------------------- ; Add continents ;----------------------------------------------------------- MAP_CONTINENTS, /NOERASE ;----------------------------------------------------------- ; 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 ENDFOR ;Different levels ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- ;COLBAR, [0.03, Py1-0.1, 2*0.22+0.22, Py1-0.05], grades, ' ' ;----------------------------------------------------------- ; Greenwich meridian plot ;----------------------------------------------------------- PRINT, 'Greenwich' xinc = y_degrees[1] - y_degrees[0] yinc = vertlevs[1] - vertlevs[0] minx = y_degrees[0] maxx = y_degrees[Ny-1] midx = (maxx+minx)/2.0 miny = vertlevs[0] maxy = vertlevs[Nz-1] midy = (maxy+miny)/2.0 minval = MIN(data(0,*,*)) maxval = MAX(data(0,*,*)) ; Develop a new title that includes other information (max, min) TitleA = Title + ' Greenwich' + $ '!C Min = ' + STRING(minval) + $ '!C Max = ' + STRING(maxval) IF (plusminus EQ 1) THEN BEGIN ;We expect the data to be roughly balanced in +/- extreme values IF (ABS(minval) GT ABS(maxval)) THEN BEGIN maxval = ABS(minval) ENDIF ELSE BEGIN minval = -ABS(maxval) ENDELSE ENDIF range = maxval-minval grades = FINDGEN(16)*range/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.25 ;Px2 = 0.7 ;Py1 = 1 - 3*0.24 - 0.16 ;Py2 = 1 - 3*0.24 Px1 = 3*0.22 + 0.03 Px2 = Px1 + 0.2 Py1 = 0.1 Py2 = 0.8 Lx = Px2-Px1 Ly = Py2-Py1 !P.POSITION = [Px1, Py1, Px2, Py2] MAP_SET, 0.0, 0.0, LIMIT=[miny,minx,maxy,maxx], TITLE=TitleA, /CYLINDRICAL, $ /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR ix = 0,(Ny-1) DO BEGIN FOR iy = 0,(Nz-1) DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(data(0,ix,iy) GE grades)) ;calculate the limits of the element to plot xmin = y_degrees[ix] - yinc/2.0 xmax = y_degrees[ix] + yinc/2.0 ymin = vertlevs[iy] - yinc/2.0 ymax = vertlevs[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 zero contour ;----------------------------------------------------------- zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 ;linecolours = INDGEN(16) & linecolours = 0 ;All contours black ;linecolours = INDGEN(16)+2 ;Contours many colours CONTOUR, data(0,*,*), 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 ENDIF ;---- ENDFOR ;----------------------------------------------------------- ; PLOT CONVERGENCE DATA ;----------------------------------------------------------- ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- DEVICE, FILE=TestDir+'Conv.eps', /PORTRAIT, /COLOR, XSIZE=18.0, YSIZE=10.0, $ XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED !P.POSITION = [0.05, 0.05, 0.95, 0.95] PLOT, INDGEN(itercount+1), relerr(0:itercount), /YLOG, YRANGE=[0.001,1.0] DEVICE, /CLOSE ;----------------------------------------------------------- ; PLOT GLOBAL AVERAGE PRESSURE DATA ;----------------------------------------------------------- ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- DEVICE, FILE=TestDir+'GlAvP.eps', /PORTRAIT, /COLOR, XSIZE=10.0, YSIZE=18.0, $ XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED !P.POSITION = [0.1, 0.1, 0.9, 0.9] PLOT, p_before(0:nlevs-1), INDGEN(nlevs) PLOT, p_after(0:nlevs-1), INDGEN(nlevs), LINESTYLE=1, /NOERASE DEVICE, /CLOSE END