; plotcrutsbulk.pro
; IDL program written by Tim Mitchell on 07.11.02

pro PlotCRUtsBulk, $
	DataFile=DataFile, $
	Filter=Filter, $
	FileType=FileType, $
	SetDivisor=SetDivisor, $
	NoMetaPlot=NoMetaPlot, $	; remove meta info from plot
	XSizeMulti=XSizeMulti, $	; default = 1.0
	YSizeMulti=YSizeMulti, $	; default = 1.0
	ForceScale=ForceScale, $	; 61(absT),62(anoT),63(absP),64(anoP)
	ForceMin=ForceMin, $		; force minimum of scale
	ForceMax=ForceMax, $		; force maximum of scale
	Manual=Manual, $		; set for manual mode
	ManYear0=ManYear0, $		; in manual, force first yr of plot to this year AD
	ManYear1=ManYear1		; in manual, force last yr of plot to this year AD

if (keyword_set(Manual) or keyword_set(DataFile) or keyword_set(filetype)) then begin
  QManualAuto = 0
endif else if (keyword_set(Filter)) then begin
  QManualAuto =-1
endif else begin
  QManualAuto = 1
endelse

if not(keyword_set(DataFile)) then DataFile = ""
if not(keyword_set(FileType)) then FileType = 1

; *************************************** intro

print, ""
print, "  > ***** PlotCRUts.pro: plots station from CRU time-series file *****"
if (QManualAuto EQ 1) then print, "  > ***** AUTOMATIC mode: you many minimise this terminal *****"
print, ""

Black = 0 & Grey = 254 & White = 255
MissVal  = -999.0 & MissCol = White
ScaleSegN = 10 & BatchFileN = 0 & GivenYear = 0 & SuspectBreak = -999
InputText = "" & GivenFile = "" & LastDumpText="standard"
RangeFract = 0.0D & ColSelect  = 0.0D
Int0=0L & Int1=0L & Int2=0L & Int3=0L & Int4=0L & Int5=0L
Int6=0L & Int7=0L & Int8=0L & Int9=0L & Int10=0L & Int11=0L
DumpWMO=0L & DumpName="" & DumpCty="" & DumpYearN=0 & DumpText="" & DumpTitle=""
MemYearAD0=-999 & MemYearAD1=-999 & QFound=-1

Months = [0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5]
XLabel = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

if (keyword_set(NoMetaPlot)) then begin
  PlotCorners  = [0.10,0.05,0.80,0.95]
  ScaleCorners = [0.85,0.05,0.95,0.95]
endif else begin
  PlotCorners  = [0.10,0.12,0.80,0.88]
  ScaleCorners = [0.85,0.08,0.95,0.88]
endelse

XPixels=12 & YPixels=12
if (keyword_set(XSizeMulti)) then XPixels = XPixels * XSizeMulti
if (keyword_set(YSizeMulti)) then YPixels = YPixels * YSizeMulti

; *************************************** load CRU ts file

if (QManualAuto EQ 1) then begin		; automatic
    spawn, 'rm /cru/tyn1/f014/scratch/f90toidl.pipe', StOut
    spawn, 'mkfifo -m a=rw /cru/tyn1/f014/scratch/f90toidl.pipe'
    spawn, 'tail < /cru/tyn1/f014/scratch/f90toidl.pipe', F90info
    DataFile = strtrim(F90info(0),2)
    QNoData = strpos(DataFile,".nrm")
    
    if (QNoData GT -1) then begin
      ; nothing required
    endif else begin
      LoadCRUtsHeads, Code,Lat,Lon,Elv,YearAD0,YearAD1,NameStn,NameCty, $
		CallFile=DataFile,Data=Data,YearAD=YearAD,TmpPre=TmpPre,Divisor=Divisor
    endelse
endif else if (QManualAuto EQ -1) then begin	; bulk
    GetFileBatch, Batch, BatchN, CallFilter=Filter
    XBatch=-1 & QNoData=-1
endif else begin				; manual
    QNoData = -1
    LoadCRUtsHeads, Code,Lat,Lon,Elv,YearAD0,YearAD1,NameStn,NameCty, $
		CallFile=DataFile,Data=Data,YearAD=YearAD,TmpPre=TmpPre,Divisor=Divisor
endelse

if (QNoData EQ -1 AND QManualAuto EQ 1) then begin
  DataSize = size(Data) & YearN = DataSize[3] & MonthN = DataSize[2] & StationN = DataSize[1]
endif else begin
  TmpPre=0 & Divisor=1.0
endelse
if (keyword_set(SetDivisor)) then begin
  Divisor=SetDivisor & print, "  > Divisor set to ", Divisor
endif

AbsScale=61 & AnoScale=62 & StDevRange=10
if (TmpPre EQ 1) then begin
  AbsScale=63 & AnoScale=64 & StDevRange=5
endif

if (keyword_set(ForceScale)) then begin
  loadanyct, ForceScale
endif else begin
  loadanyct, AbsScale
endelse

DumpMissVal = -9999.0/Divisor

; *************************************** main station loop

QCode=-1L & LastStn=0
;QCode = -1L & on_ioerror, badIO

;badIO: if (QCode NE -1L) then print, "  > Bad input. Try again."
while (QCode NE 99) do begin
  if (QManualAuto EQ 0) then begin
    read, prompt="Code (99=exit):", QCode		; get station code
    QSetYears=1
  endif else if (QManualAuto EQ -1) then begin
    XBatch=XBatch+1
    if (XBatch EQ BatchN) then begin
      QCode=99
    endif else begin
      QCode=88 & QSetYears=0
      DumpFile=Batch(XBatch)
      DumpBeg=strpos(DumpFile,'/',/reverse_search)
      DumpEnd=strpos(DumpFile,'.',DumpBeg)
      DumpStn=long(strmid(DumpFile,(DumpBeg+1),(DumpEnd-DumpBeg+1)))
      if (DumpStn NE LastStn) then begin
;        print, "stn = ", DumpStn
        QSetYears=1 & LastStn=DumpStn
      endif
    endelse
  endif else begin
    spawn, 'tail < /cru/tyn1/f014/scratch/f90toidl.pipe', F90info
    QCode = long(strtrim(F90info(0),2))
    print, ""
    print, "QCode = ", QCode
    
    if (QCode EQ 88) then begin
      spawn, 'tail < /cru/tyn1/f014/scratch/f90toidl.pipe', F90info	; get file path
      DumpFile = strtrim(F90info(0),2)
      print, "file  = ", strtrim(DumpFile,2)
      spawn, 'tail < /cru/tyn1/f014/scratch/f90toidl.pipe', F90info
      QSetYears = long(strtrim(F90info(0),2))
      print, "range = ", QSetYears
    endif else if (QCode NE 99) then begin
      spawn, 'tail < /cru/tyn1/f014/scratch/f90toidl.pipe', F90info	; get break year
      SuspectBreak = long(strtrim(F90info(0),2))
      print, "year  = ", SuspectBreak
      QSetYears = 1
    endif
  endelse
  
  if (QCode NE 99) then begin
    if (QCode NE 88) then begin
      XStation = -1 & QFound = -1				
      while (QFound EQ -1 AND XStation LT StationN-1) do begin		; locate station
        XStation = XStation + 1
        if (Code(XStation) EQ QCode) then QFound = XStation
      endwhile
      
      if (QFound EQ -1) then begin
      	print, "  > Station not found. Try again."
      endif else begin
;        XYear0 = YearAD0(QFound)-YearAD(0) & XYear1 = YearAD1(QFound)-YearAD(0)
        LastQFound  = QFound
        
        if (keyword_set(ManYear0)) then begin
          PlotYearAD0=ManYear0 & XYear0 = ManYear0-YearAD(0)
        endif else begin
          PlotYearAD0=YearAD0(QFound) & XYear0 = YearAD0(QFound)-YearAD(0)
        endelse
        
        if (keyword_set(ManYear1)) then begin
          PlotYearAD1=ManYear1 & XYear1 = ManYear1-YearAD(0)
        endif else begin
          PlotYearAD1=YearAD1(QFound) & XYear1 = YearAD1(QFound)-YearAD(0)
        endelse
        
        if (QManualAuto NE -1) then print, "Adopting range:  ", PlotYearAD0,PlotYearAD1
      endelse
    endif else begin							; load station data
      openr, lunDump, DumpFile, /get_lun
      readf, lunDump, DumpWMO, DumpName, DumpCty, DumpYearN, DumpText, format="(i8,x,a20,x,a13,x,i4,x,a8)"
      readf, lunDump, DumpTitle, format="(a80)"
      DumpData=fltarr(12,DumpYearN) & DumpData(*,*)=DumpMissVal & DumpYearAD=intarr(DumpYearN)
      
      for XDumpYear = 0, DumpYearN-1 do begin
        readf, lunDump, HYear,Int0,Int1,Int2,Int3,Int4,Int5,Int6,Int7,Int8,Int9,Int10,Int11, format="(i4,12i5)"
        DumpYearAD(XDumpYear)=HYear
        DumpData(0,XDumpYear)=float(Int0)/Divisor & DumpData(1,XDumpYear)=float(Int1)/Divisor & DumpData(2,XDumpYear)=float(Int2)/Divisor 
        DumpData(3,XDumpYear)=float(Int3)/Divisor & DumpData(4,XDumpYear)=float(Int4)/Divisor & DumpData(5,XDumpYear)=float(Int5)/Divisor 
        DumpData(6,XDumpYear)=float(Int6)/Divisor & DumpData(7,XDumpYear)=float(Int7)/Divisor & DumpData(8,XDumpYear)=float(Int8)/Divisor 
        DumpData(9,XDumpYear)=float(Int9)/Divisor & DumpData(10,XDumpYear)=float(Int10)/Divisor & DumpData(11,XDumpYear)=float(Int11)/Divisor 
      endfor
      free_lun, lunDump
      
      if (Divisor NE 10.0 AND (DumpText EQ "Homogeno" OR strmid(DumpText,4,4) EQ "Diff")) then begin
        print, "Correcting for non-standard multiplier"
        DumpData = DumpData * 10.0
        Missings = where (DumpData EQ (DumpMissVal*10.0))
        DumpData(Missings) = DumpMissVal
      endif
      
      if (MemYearAD0 EQ -999 OR QSetYears GE 0) then begin
       XYear0 =  100000000
       for XDumpYear = 0, DumpYearN-1 do begin
        for XMonth = 0, 11 do begin
          if (XYear0 GT XDumpYear AND DumpData(XMonth,XDumpYear) NE DumpMissVal) then XYear0 = XDumpYear
        endfor
       endfor
      
       XYear1 = -100000000
       for XDumpYear = DumpYearN-1, 0, -1 do begin
        for XMonth = 0, 11 do begin
          if (XYear1 LT XDumpYear AND DumpData(XMonth,XDumpYear) NE DumpMissVal) then XYear1 = XDumpYear
        endfor
       endfor
       
       if (XYear0 NE 100000000) then begin
         PlotYearAD0=DumpYearAD(XYear0) & PlotYearAD1=DumpYearAD(XYear1)
         if (QManualAuto NE -1) then print, "Adopting implicit range:  ", DumpYearAD(XYear0),DumpYearAD(XYear1)
       endif else begin
         XYear1=XYear0
       endelse
      endif else begin		; QSetYears=-1
       XYear0=0 & XYear1=DumpYearN-1
       for XDumpYear = 0, DumpYearN-1 do begin
         if (DumpYearAD(XDumpYear) EQ MemYearAD0) then XYear0=XDumpYear
         if (DumpYearAD(XDumpYear) EQ MemYearAD1) then XYear1=XDumpYear
       endfor
       PlotYearAD0=MemYearAD0 & PlotYearAD1=MemYearAD1
       if (QManualAuto NE -1) then print, "Adopting memorised range: ", MemYearAD0,MemYearAD1
      endelse
      
      Missings = where (DumpData EQ DumpMissVal, Count)
      if (Count GT 0) then DumpData [Missings] = !Values.F_NaN
    endelse
  endif
  
  if ((QCode EQ 88 OR (QCode NE 99 AND QFound NE -1))) then begin
   if (XYear1 GT XYear0) then begin
    if (QCode EQ 88) then begin
      Reformed = reform(DumpData(*,XYear0:XYear1))
      
      if            (strtrim(DumpText,2) EQ 'Homogeno') then begin
        DataMin = 0-StDevRange & DataMax = StDevRange
      endif else if (strmid(DumpText,4,4) EQ 'Diff') then begin
        DataMin = 0-(StDevRange*2) & DataMax = (StDevRange*2) 
      endif else if (strtrim(DumpText,2) EQ 'Rev-Orig') then begin
      	DataMin=OrigDataMin & DataMax=OrigDataMax
      endif else if (strmid(DumpText,4,4) EQ 'Orig' AND $
      				strmid(LastDumpText,4,4) EQ 'Comb') then begin
      	DataMin=OrigDataMin & DataMax=OrigDataMax
      endif else if (strmid(DumpText,0,4) EQ 'Anom') then begin
        DataMax = max(abs(Reformed), /nan) & DataMin = 0-DataMax
      endif else begin
        if (TmpPre EQ 0) then begin
          DataMin = min(Reformed,max=DataMax, /nan)
        endif else begin
          Moments = moment(Reformed,/nan)
          DataMin = Moments(0)-(2.0*sqrt(Moments(1)))
          DataMax = Moments(0)+(2.0*sqrt(Moments(1)))
          if (DataMin LT 0) then DataMin=0
        endelse
        OrigDataMin=DataMin & OrigDataMax=DataMax
      endelse
      
      Years = DumpYearAD(XYear0:XYear1)
      
      PlotTitle = strtrim(DumpName,2) + ' (' + strtrim(DumpCty,2) + ', ' + $
    				strtrim(string(DumpWMO),2) + ')'
      SavePath = '/cru/tyn1/f014/scratch/plotcruts.' + strtrim(string(DumpText),2) $
      			+ '.' + strtrim(string(DumpWMO),2) + '.eps'
    endif else begin
      Reformed = reform(Data(QFound,*,XYear0:XYear1))

      if (TmpPre EQ 0) then begin
          DataMin = min(Reformed,max=DataMax, /nan)
      endif else begin
          Moments = moment(Reformed,/nan)
          DataMin = Moments(0)-(2.0*sqrt(Moments(1)))
          DataMax = Moments(0)+(2.0*sqrt(Moments(1)))
          if (DataMin LT 0) then DataMin=0
      endelse

      Years = YearAD(XYear0:XYear1)

      PlotTitle = strtrim(NameStn(QFound),2) + ' (' + strtrim(NameCty(QFound),2) + ', ' + $
    				strtrim(string(QCode),2) + ')'
      DumpText='Raw-Orig'
      SavePath = '/cru/tyn1/f014/scratch/plotcruts.Original.' + strtrim(string(QCode),2) + '.eps'
    endelse
    
    if (QSetYears EQ 1) then begin		; set memory
      MemYearAD0 = PlotYearAD0 & MemYearAD1 = PlotYearAD1
      if (QManualAuto NE -1) then print, "Memorising range:         ", MemYearAD0,MemYearAD1
    endif
    
    YearADFlt = float(Years)/10.0
    YearADInt = round(YearADFlt)
    YTickVal  = where (float(YearADInt) EQ YearADFlt, YTickCount)    
    YearADInt = YearADInt * 10
    
    if (ytickcount eq 0) then begin	; no decade data range, do start/end instead
    	YearADInt = Years
    	YearADIntDim = size (YearADInt)
    	ytickcount = 2
    	ytickval = intarr(ytickcount)
    	ytickval(*) = [0,YearADIntDim[1]-1]
    endif
    
    if (not keyword_set(ForceScale) and (strtrim(DumpText,2) EQ 'Homogeno' OR strmid(DumpText,4,4) EQ 'Diff' OR $
    		strtrim(DumpText,2) EQ 'Anomalys')) then loadanyct, AnoScale
    
    if (keyword_set(ForceMin)) then DataMin = ForceMin
    if (keyword_set(ForceMax)) then DataMax = ForceMax
    if (DataMin EQ DataMax) then begin
      DataMin=DataMin-1 ; DataMax=DataMax+1
    endif
    
    GetContours, [DataMin,DataMax], [1,253], Contours, ContourColors, 21
    ContoursSparse = fltarr(5)
    for XSparse = 0, 4 do begin
      ContoursSparse(XSparse) = DataMin + (XSparse * ((DataMax-DataMin) / 4))
    endfor
    
    TrimToScale, Reformed, TrimSlice, [DataMin,DataMax], /SinglePrec
    
    Set_Plot, 'ps', /copy
    device, filename=SavePath, bits_per_pixel=8, xsize=XPixels, ysize=YPixels,/Color,/encapsulated
    
    if (keyword_set(NoMetaPlot)) then begin
      plot, [PlotYearAD0,PlotYearAD1], xrange=[0.5,11.5], yrange=[PlotYearAD0,PlotYearAD1], position=PlotCorners, $
    	ticklen=1.0,color=Black,/noerase,/nodata,charsize=0.8, $
    	xstyle=1, xticks=11, xthick=1, xtickv=months, xtickname=xlabel, font=5,  $
    	ystyle=1, yticks=YTickCount-1, ythick=1, ytickv=YearADInt[YTickVal]
    endif else begin
      plot, [PlotYearAD0,PlotYearAD1], xrange=[0.5,11.5], yrange=[PlotYearAD0,PlotYearAD1], position=PlotCorners, $
    	ticklen=1.0,color=Black,/noerase,/nodata,subtitle=PlotTitle,title=DumpText,charsize=0.8, $
    	xstyle=1, xticks=11, xthick=1, xtickv=months, xtickname=xlabel, font=5,  $
    	ystyle=1, yticks=YTickCount-1, ythick=1, ytickv=YearADInt[YTickVal]
    endelse
    
    contour, TrimSlice, Months, Years, levels=Contours, $
    	c_colors=ContourColors, /cell_fill, /overplot,/noerase

    if (keyword_set(NoMetaPlot)) then begin
      plot, [PlotYearAD0,PlotYearAD1], xrange=[0.5,11.5], yrange=[PlotYearAD0,PlotYearAD1], position=PlotCorners, $
    	ticklen=1.0,color=Black,/noerase,/nodata,charsize=0.8, $
    	xstyle=1, xticks=11, xthick=1, xtickv=months, xtickname=xlabel, font=5, $
    	ystyle=1, yticks=YTickCount-1, ythick=1, ytickv=YearADInt[YTickVal]
    endif else begin
      plot, [PlotYearAD0,PlotYearAD1], xrange=[0.5,11.5], yrange=[PlotYearAD0,PlotYearAD1], position=PlotCorners, $
    	ticklen=1.0,color=Black,/noerase,/nodata,subtitle=PlotTitle,title=DumpText,charsize=0.8, $
    	xstyle=1, xticks=11, xthick=1, xtickv=months, xtickname=xlabel, font=5, $
    	ystyle=1, yticks=YTickCount-1, ythick=1, ytickv=YearADInt[YTickVal]
    endelse
    
    if (SuspectBreak NE -999) then plots, [0.5,11.5], [SuspectBreak+0.5,SuspectBreak+0.5], $
    	linestyle=1, color=Black
    DrawConScale,ScaleCorners,Contours,ContourColors,$
    	ContoursSparse,"(f7.2)",0.8,1,"",XPixels,YPixels,0
    
    device, /close
    set_plot, 'X'	

    if (not keyword_set(ForceScale) and (strtrim(DumpText,2) EQ 'Homogeno' OR strmid(DumpText,4,4) EQ 'Diff' OR $
    		strtrim(DumpText,2) EQ 'Anomalys')) then loadanyct, AbsScale
    
    LastDumpText=DumpText & DumpText="standard"
    
    GVcommand = 'ghostview ' + strtrim(SavePath,2) + ' &'
    if (QManualAuto NE -1) then spawn, GVcommand, StOutText
   endif else begin
    if (QManualAuto NE -1) then print, "Insufficiently long period for plot."
   endelse
  endif
endwhile

; *************************************** end

print, ""
if (QManualAuto EQ 1) then spawn, 'rm /cru/tyn1/f014/scratch/f90toidl.pipe'
spawn, 'rm /cru/tyn1/f014/scratch/plotcruts*'
spawn, 'rm /cru/tyn1/f014/scratch/??????.*.txt'

end