;------------------------------------------------------------- ;+ ; NAME: ; CTM_PLOT ; ; PURPOSE: ; General plotting tool for CTM data that is stored in the ; GAMAP datainfo and fileinfo structures. CTM_PLOT can ; handle everything from 0-D to 3-D and has many different ; plot options, especially for 2-D plots for which it was ; originally designed. E.g. full map support. ; ; CATEGORY: ; Plotting ; ; CALLING SEQUENCE: ; CTM_PLOT, [ Diagn [ [ ,Keywords ] ] ; ; INPUTS: ; DIAGN -> The diagnostic number (e.g. ND22, ND45, etc or ; category name (e.g. "JV-MAP", "IJ-AVG") for which to ; read data from the punch file. ; ; KEYWORD PARAMETERS: ; Keyword parameters passed to CTM_GET_DATABLOCK: ; =============================================== ; TRACER -> Number of tracer to read from the punch file. ; ; ALTRANGE -> A vector specifying the min and max altitude ; values to be included in the extracted data block. ; ; AVERAGE -> If = 0, will not average the data ; = 1, will average data longitudinally ; = 2, will average data latitudinally ; = 4, will average data vertically ; These are cumulative (e.g. AVERAGE=3 will average over ; both lat and lon, and AVERAGE=7 will average over lat, ; lon, and vertical levels to produce 1 data point). ; ; FILENAME -> Name of the punch file to read data from. ; FILENAME is passed to CTM_GET_DATABLOCK. You can also ; use a file mask, in which case FILENAME will return the ; full filename if it is a named variable. If an empty ; filename is provided, the default search mask from ; gamap.defaults (see gamap.cmn) will be used. If no ; filename is given, ctm_plot will try ; to find the records from data already loaded. ; ; /INDEX -> If set, will interpret LAT, LEV, and LON as CTM ; indices. If not set, will interpret LAT, LEV, and LON ; as the actual values of latitude, level, and longitude. ; ; LAT -> If /INDEX is set, LAT denotes the CTM latitude index ; of the box to plot. Otherwise, LAT denotes the actual ; latitude value of that box. ; ; LEV -> An index array of sigma levels *OR* a two-element ; vector specifying the min and max sigma levels to be ; included in the plot. Default is [ 1, GRIDINFO.LMX ]. ; ; LON -> If /INDEX is set, LON denotes the CTM longitude index ; of the box to plot. Otherwise, LON denotes the actual ; longitude value of that box. ; ; PRANGE -> A vector specifying the min and max pressure ; levels to be included in the extracted data block. ; ; TOTAL -> If = 0, will not total data ; = 1, will total data longitudinally ; = 2, will total data latitudinally ; = 4, will total data vertically ; These are cumulative (e.g. TOTAL=3 will total over both ; lat and lon, and TOTAL=7 will total over lat, lon, and ; vertical levels to produce 1 data point). ; ; USE_FILEINFO -> (optional) If provided, CTM_GET_DATABLOCK will ; restrict its search to only the files that are ; contained in USE_FILEINFO which must be a FILEINFO ; structure array. Default is to use the global information ; (see gamap_cmn.pro). ; If an undefined named variable is provided in USE_FILEINFO, ; it will either contain the global FILEINFO structure array ; or the FILEINFO record of the specified file. ; USE_FILEINFO and USE_DATAINFO must be consistent, and should ; either both be used or omitted. However, it is ; possible to provide the global FILEINFO structure ; (or nothing) with a local subset of DATAINFO. ; ; USE_DATAINFO -> (optional) Restrict search to records contained ; in USE_DATAINFO which must be a DATAINFO structure array. ; If an undefined named variable is provided in USE_DATAINFO, ; it will either contain the global DATAINFO structure array ; or all DATAINFO records of the specified file. ; See also USE_FILEINFO. ; ; Keywords passed to CTM_CONVERT_UNIT: ; ==================================== ; UNIT -> Name of the unit that the DATA array will be converted ; to. If not specified, then no unit conversion will be done. ; ; Keywords passed to TVMAP: ; ========================= ; NOERASE -> Do not erase previous plot. ; ; POSITION -> A four-element array of normalized coordinates ; that specifies the location of map on the plot. POSITION ; has the same form as the POSITION keyword on a plot. ; Default is [0.1, 0.05, 0.9, 0.08]. (Passed to TVMAP). ; ; NCOLORS -> This is the maximum color index that will be used. ; ; BOTTOM -> The lowest color index of the colors to be loaded ; used in the color map and color bar. ; ; /NOCBAR -> If set, will not plot the colorbar below the map ; in the position specified by CBPOSITION. Default is to ; plot a colorbar. ; ; CBCOLOR -> Color index of the colorbar outline and ; characters. Defaults to BLACK (from colors_default.pro). ; ; CBPOSITION -> A four-element array of normalized coordinates ; that specifies the location of the colorbar. BARPOSITION ; has the same form as the POSITION keyword on a plot. ; Default is [0.1, 0.05, 0.9, 0.08]. ; ; CBUNIT -> Passes the Unit string to COLORBAR, which will be ; plotted to the right of the color bar. ; ; ; CBFORMAT -> format to use in call to colorbar. Default is I12 ; if abs(max(data)) < 1e4, else e12.2 (strings get trimmed) ; ; COLOR -> Color index of the map outline and title characters. ; Defaults to BLACK (from colors_default.pro). Also used as ; plot color for 1-D (line) plots. ; ; MPARAM -> A 3 element vector containing values for ; [ P0Lat, P0Lon, Rot ]. Default is [ 0, 0, 0 ]. ; ; POSITION -> A four-element array of normalized coordinates ; that specifies the location of the map. POSITION has ; the same form as the POSITION keyword on a plot. ; Default is [0.1, 0.1, 0.9, 0.9]. ; ; TITLE -> The title string that is to be placed atop the ; map window. ; ; /NOCONTINENTS -> If set, will suppress adding continent lines ; to the map. Default is to call MAP_CONTINENTS to plot ; continent outlines or filled boundaries. ; ; CCOLOR -> The color index of the continent outline or fill ; region. Default is BLACK (from colors_default.pro). ; ; CFILL -> Value passed to FILL_CONTINENTS keyword of ; MAP_CONTINENTS. If CFILL=1 then will fill continents ; with a solid color (as specified in CCOLOR above). ; If CFILL=2 then will fill continents with hatching. ; ; /NOGRID -> If set, will suppress printing of grid lines. ; Default is to call MAP_GRID to overlay grid lines. ; ; GCOLOR -> The color index of the grid lines. ; Default is BLACK (from colors_default.pro) ; ; /NOGLABELS -> Will suppres printing of labels for each grid ; line in TVMAP.PRO. Default is to print grid labels ; for each grid line. ; ; /NOISOTROPIC -> Will suppress plotting of an isotropic map ; (i.e. one with the same X and Y scale). Default is to ; print an isotropic map. ; ; /CONTOUR -> Will produce a line-contour map instead of the ; default color-pixel image map. ; ; /FCONTOUR -> Will produce a filled contour map instead ; of the default color-pixel image map. ; ; C_LEVELS -> Vector containing the contour levels. If not ; specified, TVMAP will use quasi-logarithmic levels. ; ; C_COLORS -> Index array of color levels for each line (or ; each fill section) of the contour map. If not ; specified, TVMAP will select default colors from the ; colortable. ; ; C_ANNOTATION -> Vector containing the contour labels. ; Default is to use string representations of C_LEVELS. ; ; C_FORMAT -> Format string used in converting C_LEVELS to ; the default C_ANNOTATION values. Default is '(f8.1)'. ; ; C_LABELS -> Specifies which contour levels should be labeled. ; By default, every other contour level is labeled. C_LABELS ; allows you to override this default and explicitly ; specify the levels to label. This parameter is a vector, ; converted to integer type if necessary. If the LEVELS ; keyword is specified, the elements of C_LABELS ; correspond directly to the levels specified, otherwise, ; they correspond to the default levels chosen by the ; CONTOUR procedure. Setting an element of the vector to ; zero causes that contour label to not be labeled. A ; nonzero value forces labeling. ; ; NOTE: If C_LABELS is given as a scalar, then it will be ; expanded to a vector with all elements having the same value. ; ; /NOLINES -> Will suppress overplotting a filled-contour map ; with lines. Default is to overlay filled-contour maps ; with lines. ; ; /NOLABELS -> Will suppress printing contour labels on both ; line-contour and filled-contour maps. ; ; OVERLAYCOLOR -> Color of the solid lines that will be ; overlaid atop a filled contour map. Default is BLACK. ; ; /OVERPLOT -> Add an additional line plot to an existing one. ; You should specify LINE=n and/or COLOR=n to distinguish ; between curves in this case. You can manually add a legend ; with legend.pro after the plot(s) are produced. ; Note that the title string will contain information on ; your first selection only. Use an explicit TITLE for ; best results. ; ; /SAMPLE -> Will cause REBIN (in TVMAP) to use nearest- ; neighbor sampling rather than bilinear interpolation. ; ; /LOG -> Will create a color-pixel plot with logarithmic ; scaling. /LOG has no effect on line-contour or ; filled-contour plots, since the default contour levels ; are quasi-logarithmic. ; ; ; Other Keywords: ; =============== ; USTR -> Unit string to be plotted to the right of the colorbar. ; If not specified, then CTM_PLOT will construct a unit ; string based on the value of TRACERINFO.UNIT. ; ; THISDATAINFO -> Returns to the calling program the THISDATAINFO ; structure obtained by CTM_GET_DATABLOCK. ; ; LABELSTRU -> Returns to the calling program the structure ; of label information obtained by CTM_LABEL. ; ; YRANGE -> range of y values for color scaling (default: ; scale each plot seperately with data min and max) ; ; _EXTRA=e -> Picks up extra keywords for CTM_GET_DATABLOCK, ; CTM_LABEL, TVMAP, and TVPLOT. ; ; OUTPUTS: ; None ; ; SUBROUTINES: ; CHKSTRU (function) COLORS_DEFAULT.PRO ; CTM_CONVERT_UNIT CTM_GET_DATABLOCK (function) ; CTM_LABEL (function) CTM_TRACERINFO ; GAMAP_CMN.PRO REPLACE_TOKEN (function) ; CTM_BOXSIZE (function) TVMAP ; CTM_DRAW_CUBE CONVERT_LON ; ; REQUIREMENTS: ; Uses GAMAP package subroutines. ; ; NOTES: ; (1) Some keywords are saved in local variables with ; slightly different names (e.g. MCONTOUR for /CONTOUR) ; to prevent confusion with routine names and/or keywords ; that are picked up by the _EXTRA=e facility. ; ; (2) Not every possible combination of keywords has been thoroghly ; tested. *PLEASE* report reproducible errors to mgs@io.harvard.edu!! ; ; (3) As of 11/17/98, CTM_PLOT can only produce X-Y plots with ; either PRESSURE or ALTITUDE along the left Y-Axis. The right ; Y-Axis is left disabled (but will fix that later on...) ; ; (4) Now define X-axis labels for longitude. The labels are ; defined correctly for data blocks that span the ; ; EXAMPLE: ; (0) To plot an ozone surface map (default) of a user-selected ; file, simply call ; CTM_PLOT ; ; (1) ; FileName = '~/terra/CTM4/run_code/ctm.pch.sep1' ; CTM_PLOT, 22, 1, FileName=FileName, Lev=[1,14], $ ; Total=4, /NoErase ; ; ; Plots vertically-summed map for tracer 1 of ND22 ; ; (J-Values map). Will erase screen before plotting map. ; ; (2) ; CTM_PLOT, 'JV-MAP-$', 1, FileName=FileName, Lev=[1,14], $ ; Total=4, /NoErase ; ; ; Same as above, but uses category name instead of ; ; diagnostic number. ; ; ; MODIFICATION HISTORY: ; bmy, 21 Sep 1998: VERSION 1.00 ; bmy, 22 Sep 1998: VERSION 1.01 ; - added AVERAGE and TOTAL keywords ; bmy, 22 Sep 1998: VERSION 1.10 ; - Modified for use with new versions of ; CTM_GET_DATABLOCK, CTM_EXTRACT, ; CTM_LABEL, REPLACE_TOKEN, and TVMAP ; bmy, 25 Sep 1998: VERSION 1.11 ; - modified for TVMAP v. 2.0 ; bmy, 28 Sep 1998: VERSION 2.00 ; - modified for TVMAP v. 2.01 ; - renamed LONSHIFT to LSHIFT ; bmy, 29 Sep 1998: - added ALTRANGE and PRANGE keywords ; (which had been previously omitted) ; bmy, 30 Sep 1998: VERSION 2.01 ; - added call to CTM_CONVERT_UNIT ; - added LABELSTRU keyword ; bmy, 07 Oct 1998 VERSION 2.02 ; - now works with TVMAP 3.0 ; - added /CONTOUR, /FCONTOUR, and ; /COLORBAR keywords ; - removed I/O error handling (that ; is already done in CTM_GET_DATABLOCK) ; bmy, 08 Oct 1998: VERSION 2.03 ; - now works with CTM_GET_DATABLOCK v. 1.03 ; and CTM_EXTRACT v. 1.04 ; - added DATA and THISDATAINFO keywords so ; that an external data block can be ; passed. ; - another bug fix for UNITSTR ; bmy, 03 Nov 1998: VERSION 2.04 ; - works with new CTM_GET_DATA, ; CTM_GET_DATABLOCK and CTM_LABEL routines ; - Now pass the FILEINFO and DATAINFO ; structures via USE_FILEINFO and ; USE_DATAINFO keywords ; - removed DATA keyword ; - changed %NAME% token to %TRACERNAME% ; - Now can pass an explicit unit string ; via the USTR keyword ; mgs, 10 Nov 1998: - adapted for use with new CTM_GET_DATABLOCK ; - changed keyword Erase to NoErase ; - defaults set to produce an OX surface map ; from IJ-AVG-$ diagnostics ; - allow for vertical cross section plots ; (interface to TVPLOT) : ** still needs work! ** ; - changed CBAR to NOCBAR ; bmy, 12 Nov 1998: - TRACER is now a keyword instead of ; an argument ; - Changed keyword CONTINENTS to ; NOCONTINENTS and GRID to NOGRID ; - added NOISOTROPIC, SAMPLE and ; keywords ; bmy, 13 Nov 1998: - VERSION 3.00 ; - ***** RENAMED to CTM_PLOT ***** ; - updated documentation header ; - renamed C_LABELS to C_ANNOTATION to ; prevent keyword name confusion ; - added NOLINES, NOLABELS, C_LABELS, ; and OVERLAYCOLOR keywords for tvmap ; - now gets default colors from ; DEFAULT_COLORS.PRO ; - Error checking for LIMIT keyword ; (OK for now...fix this later on...) ; bmy, 16 Nov 1998: - Now use %DATE% instead of %YMD1% for ; all default plot titles ; - now enhanced for TVPLOT v. 2.0 ; - now only convert units for a tracer ; if the default unit is different from ; the desired unit!! ; bmy, 17 Nov 1998: - added /PRESSURE keyword to plot pressure ; instead of altitude on the left Y-axis ; mgs, 17 Nov 1998: - messed around quite a bit, because of ; (unfortunate) changes in default_range !@#$!@ ; - added CBFormat keyword ; bmy, 23 Nov 1998: - eliminated %SCALE% token from UNITSTR, ; to be consistent with the latest ; upgrade to COLORBAR.PRO. ; - now pass SECONDS to CTM_CONVERT_UNIT ; bmy, 13 Jan 1999: - add support for line plots. Also, if ; the DATA array is averaged down to a ; single point, will print the value of ; that point and exit. ; - use NEWXTITLE and NEWYTITLE as token-replaced ; strings for XTITLE and YTITLE ; bmy, 15 Jan 1999: - add support for 3-D visualization plots ; - added unit string for contour plots ; - now place CTM_LABEL call after the case ; statement for PLOTTYPE, so that we can ; suppress printing of special characters ; in plot titles. ; bmy, 19 Jan 1999: - improve 0-D output ; - fixed [XY]Titles for line plots ; - "unitless" is now a unit string option ; - now use new default color names from ; DEFAULT_COLORS.PRO ; bmy, 20 Jan 1999: - Updated comments ; mgs, 22 Jan 1999: - couple bug fixes, some code cleaning ; - added OverPlot keyword to allow multiple ; line plots ; bmy, 19 Feb 1999: - now pass DEBUG (from GAMAP_CMN) to ; CTM_GET_DATABLOCK via DEBUG keyword ; - Rename XIND, YIND, ZIND keywords to ; XMID, YMID, ZMID, in call to CTM_GET_DATABLOCK ; bmy, 23 Feb 1999: - Add XTICKNAME, XTICKS, XTICKV in call to ; TVPLOT...fix for map regions smaller than ; the globe ; - bug fix.../NOGRID was listed as GRID!!! ; - added keyword /NOGLABELS, which ; suppresses grid line labels in MAP_SET ; - updated comments ; bmy, 26 Feb 1999: - now calls MAP_LABELS to get latitude labels ; for X, XZ, Y, YZ plot types. ; - updated comments ; bmy, 04 Mar 1999: - now pass DEBUG keyword to TVMAP ; - now use GRIDINFO.XEDGE, GRIDINFO.YEDGE ; to compute the LIMIT keyword for MAP_SET ; mgs, 18 Mar 1999: - minor cleaning ; mgs, 23 Mar 1999: - added ILun to keyword list to prevent retrieval ; of two otherwise identical records from two ; different files ; bmy, 25 Mar 1999: - now line plots use MULTIPANEL ; - if NPANELS >=2 then place the plot title ; higher above the window, to allow for ; carriage-returned lines ; for X, Y, Z, XY, XZ, YZ plots ; - if XZ or YZ plots, print unit string ; as X-axis title. The X-axis labels will ; denote latitude or longitude. ; mgs, 25 Mar 1999: -no unit conversion if not necessary ; ;- ; Copyright (C) 1998, 1999, Bob Yantosca and Martin Schultz, ; Harvard University. ; This software is provided as is without any warranty ; whatsoever. It may be freely used, copied or distributed ; for non-commercial purposes. This copyright notice must be ; kept with any copy of this software. If this software shall ; be used commercially or sold as part of a larger package, ; please contact the author to arrange payment. ; Bugs and comments should be directed to bmy@io.harvard.edu ; with subject "IDL routine ctm_plot" ;------------------------------------------------------------- pro CTM_Plot, DiagN, $ Use_FileInfo=Use_FileInfo, Use_DataInfo=Use_DataInfo, $ ThisDataInfo=ThisDataInfo, Tracer=Tracer, $ ILun=ILun, $ Average=Average, Total=FTotal, $ FileName=FileName, Index=Index, $ Lat=Lat, Lev=Lev, $ Lon=Lon, AltRange=AltRange, $ PRange=PRange, YRange=YRange, $ Unit=Unit, UStr=UStr, $ DLon=DLon, LShift=LShift, $ NoErase=NoErase, Pressure=Pressure, $ MaxData=MaxData, MinData=MinData, $ NColors=NColors, Bottom=Bottom, $ NoCBar=NoCBar, CBPosition=CBPosition, $ CBColor=CBColor, CBUnit=CBUnit, $ CBFormat=CBFormat, $ Divisions=Divisions, Sample=Sample, $ Color=MColor, MParam=MParam, $ Title=Title, Position=Position, $ NoContinents=NoContinents, CColor=CColor, $ CFill=CFill, $ NoGrid=NoGrid, GColor=GColor, $ NoGLabels=NoGLabels, $ LabelStru=LabelStru, $ Contour=MContour, FContour=FContour, $ C_Annotation=C_Annotation, C_Colors=C_Colors, $ C_Labels=C_Labels, NoLines=NoLines, $ NoLabels=NoLabels, OverLayColor=OverLayColor, $ OverPlot=OverPlot, $ _EXTRA=e ;==================================================================== ; Pass external functions ;==================================================================== FORWARD_FUNCTION ChkStru, CTM_Get_DataBlock, CTM_Label, $ Replace_Token, Get_GridSpacing, USSA_Press ; Include global common block (for DEBUG flag) @gamap_cmn ; Include symbolic names for MYCT drawing colors @default_colors ;==================================================================== ; Keyword settings ;==================================================================== if ( N_Elements( DiagN ) eq 0 ) then begin Message,'No DIAGN passed!', /Continue Message,'Will use IJ-AVG-$ ...',/INFO,/NONAME DiagN = 'IJ-AVG-$' if (n_elements(Lev) eq 0 AND n_elements(AltRange) eq 0 $ AND n_elements(PRange) eq 0 AND n_elements(Lon) eq 0 $ AND n_elements(Lat) eq 0) then Lev = 1 endif if ( N_Elements( Tracer ) eq 0 ) then begin Message,'No TRACER passed!', /Continue Message,'Will use tracer 2 as default ...',/INFO,/NONAME Tracer = 2 endif if ( N_Elements( Bottom ) eq 0 ) $ then Bottom = DEFAULT_MYCT_BOTTOM < ( !D.N_Colors-1 ) if ( N_Elements( NColors ) eq 0 ) $ then NColors = DEFAULT_MYCT_NCOLORS < ( !D.N_Colors-Bottom ) if ( N_Elements( MColor ) eq 0 ) then MColor = BLACK Pressure = Keyword_Set( Pressure ) MContour = Keyword_Set( MContour ) FContour = Keyword_Set( FContour ) ; Number of plot panels per page NPanels = !P.MULTI[1] * !P.MULTI[2] ;=================================================================== ; Call CTM_GET_DATABLOCK to read in a data block ; ; NOTE: If CTM_PLOT is being called from GAMAP.PRO, then: ; - LAT and LON contain latitude and longitude EDGES ; - LEV contains level numbers ; - XMID, YMID, ZMID contain one of the following: ; latitude CENTERS, longitude CENTERS, level CENTERS ; The ordering corresponds to the order of the dimensions ; in the data block. ;=================================================================== Success = CTM_Get_DataBlock( Data, DiagN, $ XMid=XMid, YMid=YMid, ZMid=ZMid, $ Use_FileInfo=Use_FileInfo, $ Use_DataInfo=Use_DataInfo, $ ThisDataInfo=ThisDataInfo, $ ILun=ILun, $ Average=Average, Total=FTotal, $ Tracer=Tracer, FileName=FileName, $ GridInfo=GridInfo, ModelInfo=ModelInfo, $ Index=Index, Lev=Lev, $ Lat=Lat, Lon=Lon, $ AltRange=AltRange, PRange=PRange, $ WE=WE, SN=SN, $ UP=UP, Debug=Debug, $ _EXTRA=e ) if (not Success) then return ;=================================================================== ; Store the number of dimensions of DATA in SDATA. ; If SDATA is > 3 then print error message and return. ;=================================================================== SData = Size( Data, /N_DIMENSIONS ) if ( SData gt 3 ) then begin Message, 'DATA must have <= 3 dimensions!', /Continue returnv endif ;**** Debug output (bmy, 3/4/99) if ( DEBUG ) then begin print, '### CTM_PLOT : Lon : ', Lon print, '### CTM_PLOT : Lat : ', Lat print, '### CTM_PLOT : Lev : ', Lev print, '### CTM_PLOT : WE : ', we print, '### CTM_PLOT : SN : ', sn print, '### CTM_PLOT : UP : ', up print, '### CTM_PLOT : XMid : ', Xmid print, '### CTM_PLOT : YMid : ', Ymid print, '### CTM_PLOT : ZMid : ', Zmid endif ;=================================================================== ; determine plot type: XY, XZ, YZ, etc. ; for XY plots, /CONTINENTS will be defaulted to true ;=================================================================== NLon = N_Elements( WE ) NLat = N_Elements( SN ) NAlt = N_Elements( UP ) PlotTypeList = [ '0', 'X', 'Y', 'XY', 'Z', 'XZ', 'YZ', 'XYZ' ] PlotType = PlotTypeList[ (NLon gt 1) + 2*(NLat gt 1) + 4*(NAlt gt 1) ] if (DEBUG) then print,'##DEBUG: PlotType , NLon, NLat, NAlt = ', $ PlotType,NLon,NLat,NAlt ;=================================================================== ; Call CTM_TRACERINFO to get the carbon number ; Call CTM_CONVERT_UNIT to do unit conversion ; ; NOTE: Only call CTM_CONVERT_UNIT if the default ; unit and the desired unit are different!!! ;=================================================================== CTM_TracerInfo, ThisDataInfo.Tracer, Name=Name, MolC=CNum, MolWt=MolWt ; *** QUICK FIX: ALWAYS OVERWRITE TRACER NAME IN DATAINFO STRUCTURE ThisDataInfo.TracerName = Name if ( N_Elements( Unit ) gt 0 ) then begin From_Unit = ThisDataInfo.Unit print, '### converting units from ', from_unit, ' to >', unit, '< ?' if ( From_Unit ne Unit AND Unit ne '') then begin ; Get the grid box areas in cm/2 for mol/cm2/s -> Tg ; *** HERE'S A WEAKNESS: CTM_Box_Size should only be called if necessary, ; I.E. if the unit conversion needs the box sizes ! ************ ; (perhaps we can do this in ctm_convert_units.pro altogether ?) AreaCm2 = CTM_BoxSize( GridInfo, /cm2, $ GISS_Radius = ( ModelInfo.Family eq 'GISS' )) if (n_elements(WE) gt 0) then $ if (WE[0] ge 0) then AreaCm2 = AreaCm2[ WE, * ] if (n_elements(SN) gt 0) then $ if (SN[0] ge 0) then AreaCm2 = AreaCm2[ *, SN ] ; Convert the units! ; Elapsed seconds during diagnostic interval Seconds = ( ThisDataInfo.Tau1 - ThisDataInfo.Tau0 ) * 3600.0 CTM_Convert_Unit, Data, From_Unit=From_Unit, $ To_Unit=Unit, CNum=CNum, Result=Result, $ MolWt=MolWt, AreaCm2=AreaCm2, Seconds=Seconds, $ _EXTRA=e ; Reassign Unit and Scale numbers in THISDATAINFO if (Result) then ThisDataInfo.Unit = Unit endif endif ;=================================================================== ; If YRANGE is specified, then use those values for the min and ; max of the byte scaling. Otherwise, use the min and max values ; of the DATA array. ;=================================================================== if ( N_Elements( YRange ) gt 0 ) $ then MaxData = Max( YRange, Min=MinData ) $ else MaxData = Max( Data, Min=MinData ) ;=================================================================== ; Settings (titles, etc) for the different plot types ;=================================================================== NoTit = ( n_elements(Title) eq 0 ) NoXTit = ( n_elements(XTitle) eq 0 ) NoYTit = ( n_elements(YTitle) eq 0 ) case ( PlotType ) of '0' : begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE%' if (NoXTit) then XTitle = '%LON% %LAT% %LEV% (%ALT%)' if (NoYTit) then YTitle = '' ; suppress super/subscript characters No_Special = 1 end 'X' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE% %LAT% %LEV% (%ALT%)' endif else begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE%!C!C%LAT% %LEV% (%ALT%)' endelse if (NoXTit) then XTitle = 'Longitude' if (NoYTit) then YTitle = '' ; fill in YTitle below XRange = [ min( XMid, max=M ), M ] YRange = [ min( Ymid, max=M ), M ] end 'Y' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE% %LON% %LEV% (%ALT%)' endif else begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE%!C!C%LON% %LEV% (%ALT%)' endelse if (NoXTit) then XTitle = 'Latitude' if (NoYTit) then YTitle = '' ; Fill in YTitle below XRange = [ min( XMid, max=M ), M ] YRange = [ min( Ymid, max=M ), M ] end 'Z' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE% %LAT% %LON%' endif else begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE%!C!C%LAT% %LON%' endelse if (NoXTit) then XTitle = '' ; fill in XTitle below if ( Pressure ) then begin if (NoYTit) then YTitle = 'Pressure (mb)' endif else $ if (NoYTit) then YTitle = 'Altitude (km)' XRange = [ min( XMid, max=M ), M ] YRange = [ min( Ymid, max=M ), M ] end 'XY' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE% %LEV% (%ALT%)' endif else begin if (NoTit) then Title = $ '%MODEL% %TRACERNAME% %DATE%!C!C%LEV% (%ALT%)' endelse if (NoXTit) then XTitle = '' if (NoYTit) then YTitle = '' if (n_elements(continents) eq 0) then continents = 1 end 'XZ' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE% %LAT%' endif else begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE%!C!C%LAT%' endelse if (NoXTit) then XTitle = 'Longitude' if ( Pressure ) then begin if (NoYTit) then YTitle = 'Pressure (mb)' ;YA_Title = 'Altitude (km)' endif else begin if (NoYTit) then YTitle = 'Altitude (km)' ;YA_Title = 'Pressure (mb)' endelse YRange = [ min(Ymid), max(Ymid) ] if (n_elements(LON) eq 2) then XRange = Lon end 'YZ' : begin if ( NPanels lt 2 ) then begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE% %LON%' endif else begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE%!C!C%LON%' endelse if (NoXTit) then XTitle = 'Latitude' if ( Pressure ) then begin if (NoYTit) then YTitle = 'Pressure (mb)' ;YA_Title = 'Altitude (km)' endif else begin if (NoYTit) then YTitle = 'Altitude (km)' ;YA_Title = 'Pressure (mb)' endelse YRange = [ min(Ymid), max(Ymid) ] if (n_elements(LAT) eq 2) then XRange = Lat end 'XYZ' : begin if (NoTit) then Title = '%MODEL% %TRACERNAME% %DATE%' if (NoXTit) then XTitle = '' if (NoYTit) then YTitle = '' ; suppress super/subscript characters No_Special = 1 end else : begin Message, 'Invalid Plot Type!', /Continue return end endcase ;=================================================================== ; Call CTM_LABEL to return the LABELSTRU structure ; Also pass altitude and pressure information ;=================================================================== if ( not ChkStru(GridInfo,'LMX')) then begin Prs = ModelInfo.PSurf Alt = USSA_Alt(Prs) endif else if ( N_Elements( Lev ) eq 1) then begin Alt = GridInfo.ZMid[ Lev[0] - 1 ] Prs = GridInfo.PMid[ Lev[0] - 1 ] endif else begin MinLev = ( Min( Lev, Max=MaxLev ) > 0 ) MaxLev = MaxLev < ( GridInfo.LMX - 1) Alt = [ GridInfo.ZMid( MinLev-1 ), GridInfo.ZMid( MaxLev - 1 ) ] Prs = [ GridInfo.PMid( MinLev-1 ), GridInfo.PMid( MaxLev - 1 ) ] endelse LabelStru = CTM_Label( ThisDataInfo, ModelInfo, $ Lat=Lat, Lon=Lon, Lev=Lev, Alt=Alt, Prs=Prs, $ Average=Average, Total=FTotal, $ No_Special=No_Special, _Extra=e ) ; Debug output if (DEBUG gt 1) then help,labelstru,/Stru if (DEBUG gt 1) then help,gridinfo,/stru ;=================================================================== ; Call REPLACE_TOKEN to replace tokens in the title ; strings with text from the LABELSTRU structure ;=================================================================== NewTitle = Replace_Token( Title, LabelStru ) NewXTitle = Replace_Token( XTitle, LabelStru ) NewYTitle = Replace_Token( YTitle, LabelStru ) ;=================================================================== ; Unit string for colorbar. If the unit is PPBC or PPBV ; then we don't need to print the scale factor as well. ;=================================================================== if ( ChkStru( LabelStru, [ 'UNIT' ] ) gt 0 ) then begin if ( N_Elements( UStr ) eq 0 ) then begin case ( StrUpCase( StrTrim( ThisDataInfo.Unit, 2 ) ) ) of '' : UnitStr = ' ' 'UNDEFINED' : UnitStr = '' 'UNITLESS' : UnitStr = '[unitless]' 'NONE' : UnitStr = '[unitless]' else : begin UnitStr = '[%UNIT%]' ; ### mgs, 23 Nov 1998: commented out because of problems with scale *** ; if ( ThisDataInfo.Scale eq 1.0 ) $ ; then UnitStr = '[%UNIT%]' $ ; else UnitStr = '[%SCALE% %UNIT%]' end endcase endif else begin UnitStr = UStr endelse endif ;=================================================================== ; Call REPLACE_TOKEN to replace tokens in the unit strings with ; text from the LABELSTRU structure ; ; For contour plots also add the tracername to the unit string, ; which will be placed below the plot window. ; ; For Line Plots, set the X-axis or Y-axis title to the tracer ; name plus the unit string. ;=================================================================== if ( ( PlotType eq 'XY' OR PlotType eq 'YZ' OR PlotType eq 'XZ' ) AND $ ( MContour OR FContour ) ) $ then UnitStr = 'Unit: ' + UnitStr if ( PlotType eq 'X' OR PlotType eq 'Y' OR PlotType eq 'Z' ) $ then UnitStr = '%TRACERNAME% ' + UnitStr NewUnitStr = Replace_Token( UnitStr, LabelStru ) case ( PlotType ) of 'X' : NewYTitle = NewUnitStr 'Y' : NewYTitle = NewUnitStr 'Z' : NewXTitle = NewUnitStr else : ;Null command endcase ;=================================================================== ; Call MAP_LABELS to compute the latitude and longitude ; labels that will go on the X-axis ; ; Double the default longitude spacing if there are more ; than 2 plot panels per page (bmy, 3/25/99) ;=================================================================== LonRange = [ Lon[0], Lon[1] ] LatRange = [ Lat[0], Lat[1] ] Map_Labels, LatLabel, LonLabel, $ LatRange=LatRange, LonRange=LonRange, $ Lats=Lats, Lons=Lons, $ DLat=DLat, DLon=DLon, $ _EXTRA=e ;=================================================================== ; Data block is zero-dimensional ; ; SData = 0: We are printing a single point, or have averaged ; (or totaled) a larger data cube down to a single point. ; Print value and exit (bmy, 1/13/99) ;=================================================================== if ( SData eq 0 ) then begin LineStr = $ '---------------------------------------------------------------' DataStr = 'Value of Data : ' + StrTrim( String( Data ), 2 ) Print, ' ', LineStr Print, ' ', NewTitle Print, ' ', NewXTitle Print, ' ', DataStr Print, ' ', LineStr return endif ;=================================================================== ; Data block is 1-dimensional ; ; SDATA = 1: Produce a line plot and exit ; Now use MULTIPANEL to get the right plot position ; ; NOTE: It is important to set position here and not at the top ; of the program, since position is also passed to TVMAP ; and TVPLOT. Thus, we do not let settings for line plots ; affect settings for contour or pixel plots. ;=================================================================== if ( SData eq 1 ) then begin ; For MULTIPANEL if (!D.NAME eq 'PS') then PSOffset = 0.02 else PSOffset = 0. if ( N_Elements( Position ) eq 0 ) then begin if ( MContour or FContour ) then $ Position = [ 0.05, 0.03, 1.0, 1.0 ] $ else $ Position = [ 0.05, 0.15+psoffset, 1.0, 1.0 ] ; room for colorbar endif else print,'Position passed: ',Position ;================================================================ ; Get actual position of current plot panel ; (we may want to add margin parameters at some point??) ; Use extra-large left margin, since the Y-ticks may be in ; scientific notation (bmy, 3/25/99) ;================================================================ MultiPanel, Position=PlotPosition, Margin=[0.10,0.02,0.02,0.07] NPanels = !P.MULTI[1] * !P.MULTI[2] ;============================================================== ; Calculate true window position from position and MPosition ; Here we don't need to add a colorbar ... ;============================================================== ; get width of plot window wx = ( PlotPosition[2] - PlotPosition[0] ) wy = ( PlotPosition[3] - PlotPosition[1] ) Position[0] = PlotPosition[0] + wx * Position[0] Position[1] = PlotPosition[1] + wy * Position[1] Position[2] = PlotPosition[0] + wx * Position[2] Position[3] = PlotPosition[1] + wy * Position[3] ;============================================================== ; Scale factor for charsize ;============================================================= OldXcs = !X.CHARSIZE OldYcs = !Y.CHARSIZE CsFac = 1.0 if ( NPanels gt 1 ) then CsFac = 0.9 if ( NPanels gt 4 ) then CsFac = 0.75 if ( NPanels gt 9 ) then CsFac = 0.6 if ( !D.NAME ne 'PS' ) then CsFac = CsFac * 1.2 ; Set X and Y charsize appropriate to the plot window !X.CHARSIZE = CsFac !Y.CHARSIZE = CsFac if ( PlotType eq 'Z' ) then begin if ( Pressure ) then begin ;========================================================== ; Vertical profile with pressure on Y-axis ; Keep altitude gridding, but relabel w/ pressure values ;========================================================== NewXmid = USSA_Press( XMid ) YTickName = StrTrim( String( NewXMid, Format='(i4)' ), 2 ) YTicks = N_Elements( YTickName ) - 1 if (not keyword_set(OverPlot)) then $ Plot, Data, XMid, $ Color=BLACK, XTitle=NewXTitle, $ YTitle=NewYTitle, /XStyle, /YStyle, $ YTickName=YTickName, YTicks=YTicks, YMinor=3, $ Position=Position, CharSize=CsFac, _EXTRA=e OPlot, Data, XMid, Color=MColor, _EXTRA=e endif else begin ;========================================================== ; Vertical profile with altitude on Y-axis ;========================================================== if (not keyword_set(OverPlot)) then $ Plot, Data, XMid, $ Color=BLACK, XTitle=NewXTitle, $ YTitle=NewYTitle, /XStyle, /YStyle, $ Position=Position, CharSize=CsFac, _EXTRA=e OPlot, Data, XMid, Color=MColor, _EXTRA=e endelse endif else begin ; Select X-axis labels and number of ticks if ( PlotType eq 'X' ) then begin XTickName = LonLabel XTickV = Lons Delta = DLon endif else begin XTickName = LatLabel XTickV = Lats Delta = DLat endelse XTicks = N_Elements( XTickV ) - 1 ;**** Choose appropriate minor X-tick interval (bmy, 3/25/99) case ( Delta ) of 5 : XMinor = 5 10 : XMinor = 2 else : XMinor = 3 endcase ;**** For the Y-Axis ticks YTickV = ( dindgen(10) / 10.0 ) * (MaxData-MinData) + MinData YTicks = N_Elements( YTickV ) - 1 Y_Format = get_defaultformat( YTickV[0],YTickV[YTicks], $ DefaultLen=['14.2','9.2'], Log=0) YTickName = StrTrim( String( YTickV, Format=Y_Format ), 2 ) YMinor = 2 ;============================================================= ; Horizontal profile with tracer on the X-axis and ; either latitude or longitude on the Y-axis ;============================================================= if (not keyword_set(OverPlot)) then $ Plot, Xmid, Data, $ Color=BLACK, XTitle=NewXTitle, $ YTitle=NewYTitle, /XStyle, /YStyle, $ Position=Position, CharSize=CsFac, XTickName=XTickName, $ XTicks=XTicks, XTickV=XTickV, XMinor=XMinor, $ YTickName=YTickName, YTicks=YTicks, YTickV=YTickV, $ YMinor=YMinor, _EXTRA=e OPlot, XMid, Data, Color=BLACK, _EXTRA=e endelse xpmid = (!x.window[1]+!x.window[0])/2. ;**** Place a little higher, for carriage return lines (bmy, 3/25/99) if ( NPanels lt 2 ) $ then yptop = !y.window[1]+0.025 $ else yptop = !y.window[1]+0.040 xyouts,xpmid,yptop,NewTitle,color=MColor,/norm,align=0.5, $ charsize=1.2*csfac ; Advance to next frame w/o erasing screen MultiPanel, /Advance, /NoErase ; restore old charsize for !x and !y !X.CHARSIZE = OldXcs !Y.CHARSIZE = OldYcs return endif ;=================================================================== ; Data block is 2-dimensional ; ; SDATA = 2: Plot 2-D map either as a pixel plot or as a ; contour plot, and then exit. ;=================================================================== if ( SData eq 2 ) then begin ;================================================================ ; Define some plot variables ;================================================================ Erase = 1 - Keyword_Set( NoErase ) CBar = 1 - Keyword_Set( NoCBar ) Continents = 1 - Keyword_Set( NoContinents ) Grid = 1 - Keyword_Set( NoGrid ) Isotropic = 1 - Keyword_Set( NoIsotropic ) ;================================================================ ; Make sure that the values specified for the LIMIT keyword ; for MAP_SET do not exceed the latitude and longitude ; extents specified in the GRIDINFO structure. ; This error checking will probably be moved to TVMAP later. ; ##### DON'T LIKE THIS !! MGS !!!!! ############### ; ; NOTE: Once again, if CTM_PLOT is called from GAMAP ; then LAT and LON are latitude and longitude EDGES. ;================================================================ if ( PlotType eq 'XY' ) then begin ;**** We have to go back to the end points of XEDGE ;**** and YEDGE since some grids are edged on -180 and ;**** some grids are centered on -180. Thus the placement ;**** of the map needs to be adjusted to the grid (bmy, 3/4/99) MinLat = Lat[0] > GridInfo.YEDGE[ 0 ] MinLon = Lon[0] > GridInfo.XEDGE[ 0 ] MaxLat = Lat[1] < GridInfo.YEDGE[ GridInfo.JMX ] MaxLon = Lon[1] < GridInfo.XEDGE[ GridInfo.IMX ] Limit = [ MinLat, MinLon, MaxLat, MaxLon ] LShift = 0. if (MaxLon lt MinLon) then begin Limit[3] = Limit[3] + 360. LShift = (Limit[3]+Limit[1])/2. endif DLon = GridInfo.DI endif if ( MContour OR FContour ) then Cbar = 0 if ( N_Elements( Divisions ) eq 0 ) then Divisions=4 print,'Min and Max of selected data : ',min(data),max(data) ;================================================================ ; Call TVMAP to plot the color-pixel map, ; line-contour map, or filled contour map... ; ;**** Now pass DEBUG flag to TVMAP (bmy, 3/4/99) ;================================================================ if ( PlotType eq 'XY' ) then begin TVMap, Data, XMid, YMid, $ BLACK=BLACK, DLon=DLon, $ LShift=LShift, $ ; Erase=Erase, $ MaxData=MaxData, MinData=MinData, $ NColors=NColors, Bottom=Bottom, $ CBar=CBar, CBPosition=CBPosition, $ CBColor=CBColor, CBUnit=NewUnitStr, $ CBFormat=CBFormat, $ Divisions=Divisions, Log=Log, $ Sample=Sample, Isotropic=Isotropic, $ Limit=Limit, Color=MColor, $ MParam=MParam, Title=NewTitle, $ Position=Position, $ Continents=Continents, CColor=CColor, $ CFill=CFill, $ Grid=Grid, GColor=GColor, $ NoGLabels=NoGLabels, $ Contour=MContour, FContour=FContour, $ C_Levels=C_Levels, C_Colors=C_Colors, $ C_Annotation=C_Annotation, C_Format=C_Format, $ C_Labels=C_Labels, NoLines=NoLines, $ NoLabels=NoLabels, OverLayColor=OverLayColor, $ Debug=Debug, _EXTRA=e endif $ ;================================================================ ; ... or call TVPLOT to plot the color-pixel plot, ; line contour plot, or filled-contour plot ;================================================================ else begin ; Get levels for y ticks if PRESSURE keyword set if ( Pressure ) then begin YTickV = Ymid YTicks = N_Elements( YTickV ) - 1 YTickName = StrTrim( String( GridInfo.PMid[UP], $ Format='(i4)' ), 2 ) ;YA_TickName = StrTrim( String( GridInfo.ZMid[UP], $ ; Format='(i4)' ), 2 ) ; Blank out a couple of elements of YTICKNAME, for clarity if ( N_Elements( YTickName ) gt 10 AND Min( UP ) lt 3 ) $ then YTickName[1:2] = ' ' ;if ( N_Elements( YA_TickName ) gt 10 AND Min( UP ) lt 3 ) $ ; then YA_TickName[1:2] = ' ' endif XStyle = 1 YStyle = 1 ;============================================================= ; For X or XZ plots, set XTICKV, XTICKNAME, etc ; to the quantities for longitude ;============================================================= if ( StrPos( PlotType, 'X' ) ge 0 ) then begin ; Convert to 0..360 degrees for plotting purposes if ( XRange[0] gt XRange[1] ) then begin Convert_Lon, XMid, /Pacific Convert_Lon, XRange, /Pacific endif XTickV = Lons XTickName = LonLabel XTicks = N_Elements( XTickV ) - 1 ; DLon is always 5, 10, 15, or 30, so choose an ; appropriate minor tick interval (bmy, 2/26/99) case ( DLon ) of 5 : XMinor = 5 10 : XMinor = 2 else : XMinor = 3 endcase ;**** For X-Z plots, print the unit string as the X-axis ;**** title, since the map labels already denote latitude ;**** (bmy, 3/25/99) if ( PlotType eq 'XZ' ) then NewXTitle = NewUnitStr ;============================================================= ; For Y or YZ plots, set XTICKV, XTICKNAME, etc ; to the quantities for latitude ;============================================================= endif else if ( StrPos( PlotType, 'Y' ) ge 0 ) then begin ; Number of ticks XTickName = LatLabel XTickV = Lats XTicks = N_Elements( XTickV ) - 1 ; DLon is always 5, 10, 15, or 30, so choose an ; appropriate minor tick interval (bmy, 2/26/99) case ( DLat ) of 5 : XMinor = 5 10 : XMinor = 2 else : XMinor = 3 endcase ;**** For X-Z plots, print the unit string as the X-axis ;**** title, since the map labels already denote latitude ;**** (bmy, 3/25/99) if ( PlotType eq 'YZ' ) then NewXTitle = NewUnitStr endif ;**** Debug output (bmy, 2/26/99) if ( Debug ) then begin print, '### CTM_PLOT -- before TVPLOT!' print, '### XTickV : ', XTickV print, '### XTickName : ', XTickName print, '### XTicks : ', XTicks print, '### XRange : ', XRange print, '### XMid : ', XMid print, '### YMid : ', YMid endif ;============================================================= ; Call TVPLOT with the proper keyword settings!!! ;============================================================= TVPlot, Data, XMid, YMid, $ BLACK=BLACK, $ ; Erase=Erase, $ MaxData=MaxData, MinData=MinData, $ NColors=NColors, Bottom=Bottom, $ CBar=CBar, CBPosition=CBPosition, $ CBColor=CBColor, CBUnit=NewUnitStr, $ CBFormat=CBFormat, $ Divisions=Divisions, Log=Log, $ Sample=Sample, Color=MColor, $ Title=NewTitle, Position=Position, $ XRange=XRange, YRange=YRange, $ XTitle=NewXTitle, YTitle=NewYTitle, $ XStyle=XStyle, YStyle=YStyle, $ Contour=MContour, FContour=FContour, $ C_Levels=C_Levels, C_Colors=C_Colors, $ C_Annotation=C_Annotation, C_Format=C_Format, $ C_Labels=C_Labels, NoLines=NoLines, $ NoLabels=NoLabels, OverLayColor=OverLayColor, $ YTickV=YTickV, YTicks=YTicks, $ YTickName=YTickName, YA_Title=YA_Title, $ YA_TickName=YA_TickName, $ ;**** fix for longitudes spanning the date line (bmy, 2/23/99) XTickName=XTickName, XTicks=XTicks, $ XTickV=XTickV, XMinor=XMinor, $ _EXTRA=e endelse return endif ;=================================================================== ; Data block is 3-dimensional ; ; SDATA = 3: visualize the 3-D data cube and then exit. ; Since the data slicer modifies the color tables, restore ; the original coloratable after quitting CTM_DRAW_CUBE. ;=================================================================== if ( SData eq 3 ) then begin CTM_Draw_Cube, Data, _EXTRA=e return endif ;=================================================================== ; If we have reached here, then this is an invalid plot type. ; Print error message and return. ;=================================================================== Message, 'Data block cannot have more than 3 dimensions.', /Continue return end