pro calrunavg2_wc,n,cal,caltemp,mandb,mandboutd,mandboutds,table ;Modified by BWBerger 7/12/99 to handle a new format for cal - ;it used to have the date as yymmdd in the first column but was changed ;to three columns yyyy mm dd. All mastercal files were updated so the previous ;format is never used. ;***CALRUNAVG2 returns a list of daily slopes and intercepts given the ;***slopes and intercepts contained in the raw calibration files. The list is ;***based on a weighted running average of the raw data. ;*** ;***Input n - number of half days to average over (need 5 good values ;*** in the period to calculate an average value) ;*** cal - raw calibration values: mastercal_I_L.dat files ;*** I=c for CO2, =q for H2O; L= level ;*** ;*** ;***Output caltemp - an array based on cal but has all dates. ;*** Badvals (-999.) are inserted for all dates with no ;*** data. ;*** mandb - an array of [index,m,b]. The points are based ;*** on an average over n points in CALTEMP and the value ;*** placed at the center of the interval (unless there ;*** are insufficient points to average - then a badval ;*** is placed there). ;*** If n is a constant, mandb would be n points shorter ;*** than caltemp since the first averaged value ;*** corresponds to the index=(n+1)/2 of caltemp. To make ;*** mandb the same length as caltemp the ends are ;*** determined by an average taken over a window width that ;*** decreases by 2 pts for each step taken toward the end. ;*** The ends themselves are set to the average of the ;*** ending 3 points (which is the same as the second to ;*** last point). ;*** mandbout - an array exactly like mandb, but it has been repaired ;*** for holes by extrapolating between the endpoints. ;*** The endpoints are generally based on the average of ;*** 3 points at each end of the gap. If only one or two ;*** points exist, the first point next to the gap is used ;*** for the extrapolation. If a gap is at the very ;*** beginning or very end of mandb, usually 3 points to ;*** the right or left, respectively, are averaged and ;*** the entire gap is given the avg. value. If only one or ;*** two points are available, the first point to the ;*** right and left, respectively, are used. ;*** Format: [index, m, b],where m and b correspond ;*** to indexes which are for each half-day ;*** mandboutd - an array based on mandbout except that m and b are ;*** listed for the day (not the half day). ;*** Format: [Julian day,m,b], where m and b ;*** correspond to mid-day of the Julian day ;*** mandboutds - an array based on mandboutd except that the m ;*** and b values are modified to conform with earlier ;*** definitions still used in WLEF data processing codes. ;*** Format: [Julian day,slope=1/m,intercept=-b/m], where ;*** slope and intercept correspond to mid-day of Julian day ;*** table - a table of the gaps in mandb ;*** format:[blocknumber,leftindex,rightindex] ;*** (if no gaps, table=-1) ;*** ;***Written by BWBerger 8/24/98 ; modified by weiguo wang 11/3/99 ;get the first and last day from the input array CAL jdatestart=jdate2(cal(0,0),cal(1,0),cal(2,0)) jdatestop=jdate2(cal(0,row(cal)-1),cal(1,row(cal)-1),cal(2,row(cal)-1))+cal(3,row(cal)-1) print,'Starting julian day:',jdatestart print,'Stopping julian day (up to):',jdatestop ndy=jdatestop-jdatestart ;the number of days in cal ;CAL contains dates and ndays for fits. We create a new matrix CALTEMP ;with indexes for every 12 hours. Then we can insert the CAL data into ;CALTEMP at its proper time location. After CALTEMP is filled it can be ;filled in for gaps and then the calibrations for each mid-day can be ;put in the MANDBOUT matrix. ;First initialize CALTEMP ;All values other than the first column are -999. nct=(jdatestop-jdatestart+1)*2-2 ;the max index of caltemp caltemp=fltarr(5,nct+1)-999. caltemp(0,*)=lindgen(row(caltemp)) ;put index in 1st col of caltemp if row(caltemp) lt n then begin print,'Averaging window width, n, too large. Must be less than:',$ row(caltemp) goto, finish endif ;load each row of caltemp where there are corresponding values from cal for j=0,len(cal)-1 do begin i=(jdate2(cal(0,j),cal(1,j),cal(2,j))-jdatestart+1)*2-2+cal(3,j) caltemp(1:4,i)=[cal(4:5,j),cal(8:9,j)] endfor ;Now create MANDB matrix of index,slope and intercept ;It is caltemp with weighted averages for each half day index ;The window width of the average is n mandb=fltarr(3,row(caltemp)) mandb(0,*)=lindgen(row(mandb)) ;***Load up center of mandb (i.e., where window width n is constant) for i=0,len(caltemp)-n do begin mtemp=fltarr(2,n) ;columns for ma,sigmasq btemp=fltarr(2,n) ;columns for ba,sigbasq mtemp=[caltemp(1,i:i+n-1),caltemp(3,i:i+n-1)] btemp=[caltemp(2,i:i+n-1),caltemp(4,i:i+n-1)] rm=where(mtemp(0,*) ne -999.,cm) if cm ge 5 then begin ;i.e., if there are at least 5 good m values in mtemp mandb(1,((n-1)/2)+i)=total(mtemp(0,rm)/mtemp(1,rm))/total(1/mtemp(1,rm)) endif else begin mandb(1,((n-1)/2)+i)=-999. endelse rb=where(btemp(0,*) ne -999.,cb) if cb ge 5 then begin ;i.e., if at least 5 good b values in mtemp mandb(2,((n-1)/2)+i)=total(btemp(0,rb)/btemp(1,rb))/total(1/btemp(1,rb)) endif else begin mandb(2,((n-1)/2)+i)=-999. endelse endfor ;***Load up ends of mandb (i.e., where window width n is not constant) for nn=3,n-2,2 do begin mtemp1=fltarr(nn) btemp1=fltarr(nn) mtemp2=fltarr(nn) btemp2=fltarr(nn) mtemp1=[caltemp(1,0:nn-1),caltemp(3,0:nn-1)] btemp1=[caltemp(2,0:nn-1),caltemp(4,0:nn-1)] mtemp2=[caltemp(1,len(mandb)-nn:len(mandb)-1),caltemp(3,len(mandb)-nn:len(mandb)-1)] btemp2=[caltemp(2,len(mandb)-nn:len(mandb)-1),caltemp(4,len(mandb)-nn:len(mandb)-1)] rm1=where(mtemp1(0,*) ne -999.,cm1) if cm1 ge 2 then begin ;i.e., if there are at least 2 good m values in mtemp1 mandb(1,((nn-1)/2))=total(mtemp1(0,rm1)/mtemp1(1,rm1))/total(1/mtemp1(1,rm1)) endif else begin mandb(1,((nn-1)/2))=-999. endelse rm2=where(mtemp2(0,*) ne -999.,cm2) if cm2 ge 2 then begin ;i.e., if there are at least 2 good m values in mtemp2 mandb(1,(len(mandb)-1)-((nn-1)/2))=total(mtemp2(0,rm2)/mtemp2(1,rm2))/total(1/mtemp2(1,rm2)) endif else begin mandb(1,(len(mandb)-1)-((nn-1)/2))=-999. endelse rb1=where(btemp1(0,*) ne -999.,cb1) if cb1 ge 2 then begin ;i.e., if at least 2 good b values in mtemp1 mandb(2,((nn-1)/2))=total(btemp1(0,rb1)/btemp1(1,rb1))/total(1/btemp1(1,rb1)) endif else begin mandb(2,((nn-1)/2))=-999. endelse rb2=where(btemp2(0,*) ne -999.,cb2) if cb2 ge 2 then begin ;i.e., if at least 2 good b values in mtemp2 mandb(2,(len(mandb)-1)-((nn-1)/2))=total(btemp2(0,rb2)/btemp2(1,rb2))/total(1/btemp2(1,rb2)) endif else begin mandb(2,(len(mandb)-1)-((nn-1)/2))=-999. endelse endfor ;endpoints themselves are set to next point in that should be based ;on a three point average. mandb(1,0)=mandb(1,1) mandb(1,len(mandb)-1)=mandb(1,len(mandb)-2) mandb(2,0)=mandb(2,1) mandb(2,len(mandb)-1)=mandb(2,len(mandb)-2) ;*************************************************************** ;Now all the averaged values are in place. Next, all the gaps ;that could not be bridged with the running average window (provided ;that there were five good values in the window) will be filled. ;*************************************************************** ;***Define an r vector containing all the indicies of bad points in mandb r=where(mandb(1,*) eq -999.,c) if (c ne 0) and (c ne row(mandb)) then begin ;***Create a table of all the bad blocks of points in r if len(r) eq 1 then begin print,'Only one bad pt in series.' table=[0,r(0),r(0)] endif else begin table=fltarr(3,100) i=1 done=0 blockcount=1 blocksize=1 while (done ne 1) do begin if r(i) eq r(i-1)+1 then begin blocksize=blocksize+1 if (i eq len(r)-1) then begin done=1 ;add last blk to table table(*,blockcount-1)=[blockcount,r(i-blocksize+1),r(i)] table=table(*,0:blockcount-1) endif else begin i=i+1 endelse endif else begin table(*,blockcount-1)=[blockcount,r(i-blocksize),r(i-1)] if i eq len(r)-1 then begin table(*,blockcount)=[blockcount+1,r(i),r(i)] table=table(*,0:blockcount) done=1 endif else begin blockcount=blockcount+1 blocksize=1 i=i+1 endelse endelse endwhile endelse ;***Check first block for low and high end points to extrapolate ;***across the block repairtype=0 ;check low end of first block if table(1,0) ge 3 then begin ;at least 3 low end pts to avg for extrap sl=2 dl=3 print,'First gap block - at least 3 low end pts' endif if (table(1,0) eq 1) or (table(1,0) eq 2) then begin ;only 1 or 2 low end pts sl=1 dl=1 print,'First gap block - only 1 or 2 low end pts' endif if table(1,0) eq 0 then begin ;no low end pts sl=0 dl=0 repairtype=1 print,'First gap block - no low end pts' endif ;check high end of first block if row(table) gt 1 then begin ;there's more than 1 block print,'There is more than one block.' if table(1,1)-table(2,0) ge 4 then begin ;enough high end pts sh=2 dh=3 print,'First gap block - at least 3 high end pts' endif if (table(1,1)-table(2,0) eq 3) or (table(1,1)-table(2,0) eq 2) then begin sh=1 ; only 1 or 2 pts dh=1 print,'First gap block - only 1 or 2 high end pts' endif endif else begin ;First block is also the last block if table(2,0) le row(mandb)-1-3 then begin ;enough high end pts sh=2 dh=3 print,'First & Last gap block - at least 3 high end pts' endif if (table(2,0) eq row(mandb)-1-2) or $ (table(2,0) eq row(mandb)-1-1) then begin sh=1 ;only 1 or 2 pts dh=1 print,'First & Last gap block - only 1 or 2 high end pts' endif if table(2,0) eq row(mandb)-1 then begin ;no high end pts sh=0 dh=0 repairtype=2 print,'First & Last gap block - no high end pts' endif endelse ;repair mandb extrapo,mandb,table(1,0),table(2,0),sl,sh,dl,dh,repairtype,mandbout if row(table) eq 1 then begin goto, finish endif ;***If there was more than one block, must now check last block. repairtype=0 ;check low end of last block if table(1,row(table)-1)-table(2,row(table)-1-1) ge 4 then begin sl=2 ; enough low end pts dl=3 print,'Last gap block - at least 3 low end pts' endif if (table(1,row(table)-1)-table(2,row(table)-1-1) eq 2) $ or (table(1,row(table)-1)-table(2,row(table)-1-1) eq 3) then begin sl=1 ;only 1 or 2 low end pts dl=1 print,'Last gap block - only 1 or 2 low end pts' endif ;check high end of last block if table(2,row(table)-1) le row(mandb)-1-3 then begin ;enough high end pts sh=2 dh=3 print,'Last gap block - at least 3 high end pts' endif if (table(2,row(table)-1) eq row(mandb)-1-2)$ or (table(2,row(table)-1) eq row(mandb)-1-1) then begin sh=1 ;only 1 or 2 high end pts dh=1 print,'Last gap block - only 1 or 2 high end pts' endif if table(2,row(table)-1) eq row(mandb)-1 then begin ;no high end pts sh=0 dh=0 repairtype=2 print,'Last gap block - no high end pts' endif ;repair last block extrapo,mandbout,table(1,len(table)-1),table(2,len(table)-1),sl,sh,dl,dh,$ repairtype,mandbout if row(table) eq 2 then begin ;only 2 blocks to fix print,'Only two gap blocks' goto, finish endif ;***If more than two blocks, check middle blocks repairtype=0 for i=1,row(table)-2 do begin ;check low end of block if table(1,i)-table(2,i-1) ge 4 then begin ; enough low end pts sl=2 dl=3 print,'Middle gap block - at least 3 low end pts' endif if (table(1,i)-table(2,i-1) eq 3) or (table(1,i)-table(2,i-1) eq 2) then begin sl=1 ;only 1 or 2 pts dl=1 print,'Middle gap block - only 1 or 2 low end pts' endif ;check high end of block if table(1,i+1)-table(2,i) ge 4 then begin ;enough high end pts sh=2 dh=3 print,'Middle gap block - at least 3 high end pts' endif if (table(1,i+1)-table(2,i) eq 3) or (table(1,i+1)-table(2,i) eq 2) then begin sh=1 ; only 1 or 2 pts dh=1 print,'Middle gap block - only 1 or 2 high end pts' endif ;repair middle block extrapo,mandbout,table(1,i),table(2,i),sl,sh,dl,dh,repairtype,mandbout endfor endif else begin mandbout=mandb if c eq 0 then begin table=-1 print,'No gaps.' endif else begin table=-2 print,'No good points in cal series.' endelse endelse finish: ;Extract the values for the mid-day of each day and create output table ;The first column is the jdate, the second m, and the third b mandboutd=fltarr(3,ndy) for i=0,ndy-1 do begin mandboutd(*,i)=[jdatestart+i,mandbout(1,i*2+1),mandbout(2,i*2+1)] endfor mandboutds=fltarr(3,ndy) for i=0,ndy-1 do begin mandboutds(*,i)=[mandboutd(0,i),1/mandboutd(1,i),-mandboutd(2,i)/mandboutd(1,i)] endfor return end