;------------------------------------------------------------- ;+ ; NAME: ; TSDIAG ; ; PURPOSE: ; Reads and plots CTM time series data ; ; CATEGORY: ; CTM Tools ; ; CALLING SEQUENCE: ; TSDIAG [,RESULT] [,keywords] ; ; INPUTS: ; None ; ; KEYWORD PARAMETERS: ; UNIT -> select output unit (default "pptv") ; ; /PS -> Selects the PostScript device. ; ; FILENAME -> Path of the input file containing CTM ; Data data. ; ; OUTFILENAME -> Path for the PostScript file. ; ; PRINTER -> Name of the PostScript printer to send output to. ; ; WINPARAM -> A 3-element array: ; [ Window Number, Window X-Size, Window Y-Size] ; WINPARAM is used as input to subroutine OPEN_DEVICE. ; ; /SMALLCHEM -> Denotes simple NOX-OX-CO chemistry mechanism. ; The tracer indices are different for the small ; mechanism than for the full mechanism. ; ; /SEQUENCE -> Pauses after each plot, so that the user ; has a chance to look at the screen. ; ; /LOCALTIME -> Plots local time along the X-axis. Default is ; to plot number of time steps along the X-axis. ; ; COLS -> Number of COLS per page, used to set !P.MULTI. ; Default is 4. ; ; ROWS -> Number of rows per page, used to set !P.MULTI. ; Default is 4. ; ; /PPBC -> Will report concentration of certain Non-methane ; hydrocarbons (ALK4, ISOP, PRPE, C3H8, C2H6) ; as ppb Carbon. The default is to report these ; concentrations as ppbv (mol/mol) tracer. ; ; INTERVAL -> If /LOCALTIME is set then INTERVAL specifies ; the time interval in hours between major tick ; marks. Default is 24 hrs. ; ; _EXTRA -> Picks up any extra keywords for the PLOT command. ; ; OUTPUTS: ; RESULT -> A structure containing LON, LAT, ALT and Data data ; together with TRACER and a "global" TAU0 and TAU1 value. ; ; SUBROUTINES: ; OPEN_DEVICE CLOSE_DEVICE ; CTM_TRACERINFO FILE_EXIST (function) ; MYCT CLEANPLOT ; STRSCINOT (function) ; ; REQUIREMENTS: ; (1) The CTM Data files are produced by FORTRAN subroutines ; diag48.f and diag5.f. ; ; (2) If /LOCALTIME is set, then Tracer #11 (Local Time) must ; be printed out to the CTM Data file for each ; station. ; ; NOTES: ; ONLinesy the Data data (Header = 'TB' or 'DV') will be ; read from disk. Statistics are ignored for now. ; ; EXAMPLE: ; tsdiag, FILENAME='ctm.ts', /LOCALTIME, /SEQUENCE ; ; MODIFICATION HISTORY: ; bmy, 06 May 1998: VERSION 1.00 ; bmy, 07 May 1998: VERSION 1.01 ; - added PPBC and INTERVAL keywords ; - now calls FILE_EXIST to make sure ; the input file exists ; bmy, 27 May 1998 - now uses CTM_DIAGINFO to return ; the proper tracer offset ; bmy, 28 May 1998 - now uses SCALE, UNIT, and MOLC information ; as returned by CTM_TRACERINFO. ; - now uses EXP_STR to compute a ; power-of-ten string for the plot title ; bmy, 29 May 1998 - now calls CLEANPLOT to initialize ; all system variables ; bmy, 02 Jun 1998 - now uses STRSCI and STRCHEM string ; formatting functions ; mgs, 11 Jun 1998: - couple of bug fixes ; mgs, 15 Jun 1998: - default tick interval now 48 h ; mgs, 22 Jun 1998: - added Data and CSTEPS keywords ; mgs, 20 Nov 1998: - now uses convert_unit ; hsu, 22 Mar 1999: - MaxSteps now 15000 (prev. 5400) line 148 ; hsu, 23 Mar 1999: - Changed line 218 due to formatting change ;- ; Copyright (C) 1998, Bob Yantosca, 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 tsdiag" ;------------------------------------------------------------- pro tsdiag, Result, FileName=FileName, Scale=Scale, Plot=DoPlot Result = -1L ; default if (n_elements(Scale) eq 0) then Scale = 1.e12 ; pptv as default ; Default Filename DefaultFile = '~/terra/CTM4/*.ts.*' if ( n_elements(FileName) eq 0 ) then $ FileName = DefaultFile open_file,FileName,Ilun,default=DefaultFile,FileName=TrueFile, $ title='Choose a time-series file' ; Check to see if the input file exists before proceeding if ( ilun le 0 ) then begin message,'Could not open input file '+FileName+'!',/Cont return end ; return true filename as FileName FileName = TrueFile ; Inititalize array variables MaxStations = 700 MaxSteps = 20500 Data = fltarr( MaxSteps, MaxStations ) Header = strarr( MaxStations ) Lat = intarr( MaxStations ) Lon = intarr( MaxStations ) Alt = intarr( MaxStations ) Tracer = intarr( MaxStations ) ; MS = intarr( MaxStations ) CSteps = lonarr( MaxStations ) ; Initialize scalar variables Tau0 = 0L Tau1 = 0L ; Call CTM_DIAGINFO to get information about ; this diagnostic (ND48 = Data) CTM_DiagInfo, 48, Offset=Offset ;==================================================== ; READ Data DATA ; * Loop through file ; * Whenever a time series label is detected (TS, DV in ; old format, 'TIME-SERIES', 'DEPVEL' in new format) ; parse header line and read in data block ; * add as new station if location and tracer number ; were not encountered before, append if data block ; is continuation. ;==================================================== First = 1 NStations = 0 while ( not EOF( Ilun ) ) do begin Line = '' readf, Ilun, Line ; ---------------------------------------------------------------- ; #### NEW FORMAT ##### mgs, 25 Nov 1998 ; TIME-SERIES indicates time series concentration data ; DEPVEL denoted deposition velocities etc. ; ---------------------------------------------------------------- on_ioerror,try_old_format thislabel = '' I = 0 & J = 0 & L = 0 & N = 0 & MS = 0 AScale = 0.0 ATau0 = 0L ATau1 = 0L reads,Line,thislabel,NLines,I,J,L,N,MS,AScale,ATau0,ATau1, $ format='(1X,A12,6I5,E10.3,2I8)' goto,read_data ; ---------------------------------------------------------------- ; #### OLD FORMAT ##### ; ---------------------------------------------------------------- try_old_format: Line = strtrim(Line,2) ; test if this is a valid header line if ( strmid( Line, 0, 2 ) ne 'TB' AND $ strmid( Line, 0, 2 ) ne 'DV' ) then goto,Next_Iteration ; define additional dummy variables sdum = '' fdum1 = 0.0 & fdum2 = 0.0 & fdum3 = 0.0 on_ioerror,Next_Iteration reads, Line, sdum, $ I, J, L, MS, N, NLines, ATau0, ATau1, AScale, $ format='( a2, 6I5, 2I10, e14 )' print,'Scaling factor for this record:',AScale ; ---------------------------------------------------------------- ; Analyse input from header and loop through data lines ; ---------------------------------------------------------------- read_data: on_ioerror,NULL if (First) then begin Tau0 = ATau0 First = 0 endif Tau1 = ATau1 ; will always be at least previous value ; check whether that same station has already been read test = where(Lon eq I AND Lat eq J AND Alt eq L AND Tracer eq N) if (test[0] ge 0) then begin test = test[0] step = CSteps[test] Station = test print,'Continue station ',strtrim(Station,2) endif else begin step = 0L Station = NStations NStations = NStations + 1 if (NStations ge MaxStations) then begin message,'Maximum number of stations exceeded!',/Cont return endif Header[Station] = strtrim(thislabel,2) Lon[Station] = I Lat[Station] = J Alt[Station] = L Tracer[Station] = N ; MS[Station] = MS print,'New Station (lon,lat,alt,tracer):',I,J,L,N ; if (NStations eq 120) then goto,Abort_Reading endelse for i=1,NLines do begin if (eof(ilun)) then goto,Next_Iteration ; will end loop readf,Ilun,Line DumArr = fltarr(12) - 9.99E31 on_ioerror,ignore_it ; takes care of incomplete lines reads, Line, DumArr, format='(12f6.2)' ignore_it: on_ioerror,NULL ind = where(dumarr gt -9.99E31,COUNT) for J = 0, COUNT-1 do begin Data(Step, Station) = DumArr(J) * AScale Step = Step+1 endfor CSteps[Station] = step endfor Next_Iteration: on_ioerror,NULL endwhile Abort_Reading: ; Close file free_lun,ilun ; Test if any data was read successfully if (NStations eq 0) then begin message,'No data was read!',/Continue return endif print,'Stations:',Nstations,' maximum number of STEPS: ', $ max(csteps[0:NStations-1]) print,'Tau0, Tau1, delta-Tau :',Tau0,Tau1,Tau1-Tau0, $ ' = ',Tau2YYMMDD(Tau0,/SHORT),'-',Tau2YYMMDD(Tau1,/SHORT) ; Strip arrays MaxSteps = max(CSteps) Data = Data[*,0:NStations-1] Lon = Lon[0:NStations-1] Lat = Lat[0:NStations-1] Alt = Alt[0:NStations-1] Tracer = Tracer[0:NStations-1] CSteps = CSteps[0:NStations-1] Data = Data[0:MaxSteps-1,*] ; Rescale Data Data = Scale * Data ; Create Time vector Time = findgen(MaxSteps)/MaxSteps*(Tau1-Tau0)/24. ; ===================================================== ; Prepare result structure ; ===================================================== result = { Lon:Lon, Lat:Lat, Alt:Alt, Tracer:Tracer, $ TAU0:tau0, TAU1:tau1, $ Time:Time, DATA:Data } if (keyword_set(DoPlot)) then begin ; quick and dirty quality control plots ; get unique station identifiers ulon = lon(uniq(lon,sort(lon))) ulat = lat(uniq(lat,sort(lat))) utra = tracer(uniq(tracer,sort(tracer))) nplots = n_elements(ulon)*n_elements(ulat) ; plots per page nx = (sqrt(nplots-1)+1) < 5 ny = (sqrt(nplots-1)+1) < 4 !p.multi=[0,nx,ny] !x.margin=[2,1] !y.margin=[2,1] for n=0,n_elements(utra)-1 do begin for i=0,n_elements(ulon)-1 do begin for j=0,n_elements(ulat)-1 do begin ind = where(Lon eq ulon[i] AND Lat eq ulat[j] AND Tracer eq utra[n]) if (ind[0] lt 0) then goto,skip_this ; sort altitudes salt = sort(Alt[ind]) nalt = n_elements(salt) ; get max data maxdat = max(data[*,ind]) pw = 10.^(-fix(alog10(maxdat)-0.9999 ) ) maxdat = fix(maxdat*pw)/pw ; draw coordinate frame plot,time,data[*,0],color=1,yrange=[0.,(nalt+1)*maxdat*0.5], $ title=strtrim(ulon[i],2)+','+strtrim(ulat[j],2)+':'+ $ strtrim(utra[n],2) for l=0,nalt-1 do begin oplot,time,Data[*,ind[salt[l]]]+0.5*(l+1)*maxdat, $ color=1,line=(l mod 2)*2 endfor skip_this: endfor endfor endfor endif return end