; program to construct cloud correlation coefficients (with sunshine %)
;	assumes that the relationship in each location is same throughout year
; method approximately follows New et al 2000
; this program is required because Mark New has lost both
;	the correlation data file, and construction files
; written by Tim Mitchell 10.01.03

pro CloudCorrSPCAnn9196

MissVal=-999.0 & Ratio=1.0

NYear=6
MonthNum = strarr(12)
MonthNum = ['01','02','03','04','05','06','07','08','09','10','11','12']
CLD = fltarr (NYear,12,72,36) & SPC = fltarr (NYear,12,72,36)
Aye = fltarr (72,36) & Bee = fltarr (72,36)
Aye(*,*) = MissVal & Bee(*,*) = MissVal

print,"loading..."
for XYear = 0, NYear-1 do begin		; load the data
  YearAD = 1991 + XYear
  
  for XMonth = 0, 11 do begin
    fileCLD = "/cru/tyn1/f709762/scratch/cld9196." + strtrim(MonthNum(XMonth),2) + "." + $
    			strtrim(string(YearAD),2) + ".glo"
    loadglo,23,CLDmon,fileCLD,fileinfo,/silent
    CLD (XYear,XMonth,*,*) = CLDmon (*,*)
    
    fileSPC = "/cru/tyn1/f709762/scratch/spc9196." + strtrim(MonthNum(XMonth),2) + "." + $
    			strtrim(string(YearAD),2) + ".glo"
    loadglo,23,SPCmon,fileSPC,fileinfo,/silent
    SPC (XYear,XMonth,*,*) = SPCmon (*,*)
  endfor
endfor

print,"regressing..."
for XLon = 0, 71 do begin		; now do the correlations
 for XLat = 0, 35 do begin
    SPCvec = reform(SPC(*,*,XLon,XLat),NYear*12)
    CLDvec = reform(CLD(*,*,XLon,XLat),NYear*12)
    robust_reg_tdm,SPCvec,CLDvec,Kay1,Kay2,YFit,-1,miss=0.0,Ratio=Ratio
    if finite(Kay1) and finite(Kay2) and finite(Ratio) then begin
      if Ratio LT 0.3 then begin
        Aye(XLon,XLat) = Kay1 & Bee(XLon,XLat) = Kay2
      endif
    endif
 endfor
endfor

print,"saving 5.0 ..."
for XMonth = 0, 11 do begin
  savefileA="/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/a.9196." + strtrim(MonthNum(XMonth),2) + ".txt"
  savefileB="/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/b.9196." + strtrim(MonthNum(XMonth),2) + ".txt"
  openw, lunSaveA, savefileA, /get_lun
  openw, lunSaveB, savefileB, /get_lun
  for XLon = 0, 71 do begin
    for XLat = 0, 35 do begin
      if (Aye(XLon,XLat) NE MissVal AND Bee(XLon,XLat) NE MissVal) then begin
      	if (abs(Aye(XLon,XLat)) LE 99999 AND abs(Bee(XLon,XLat)) LE 99999) then begin
      	  lat = ((float(xlat)+0.5) * 5.0) -  90.0
      	  lon = ((float(xlon)+0.5) * 5.0) - 180.0
      	  printf,lunSaveA,lat,lon,1.0,Aye(XLon,XLat),1.0,format="(2f8.2,f8.1,f13.5,i7)"
      	  printf,lunSaveB,lat,lon,1.0,Bee(XLon,XLat),1.0,format="(2f8.2,f8.1,f13.5,i7)"
        endif
      endif
    endfor
  endfor
  free_lun, lunSaveA
  free_lun, lunSaveB
endfor

print,"interpolating to 2.5 ..."
quick_interp_tdm2,9196,9196,"/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/a.25.",1000,/info, $
   binfac=1000.0,gs=2.5,pts_prefix="/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/a.",/dumpbin,/dumpglo
quick_interp_tdm2,9196,9196,"/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/b.25.",1000,/info, $
   binfac=1000.0,gs=2.5,pts_prefix="/cru/tyn1/f709762/cru_ts_2.0/_constants/_9196/spc2cld/_ann/b.",/dumpbin,/dumpglo

end