;
; WRITTEN:  15JUN2007   Xuepu Zhao for DRC 
; MODIFIED: 24MAY2009   Xuepu Zhao
;
pro drcmsc_updmsc20090810,crn,cl0,itl,updmsc0,updmscl,updmscr,updscrnphd0,updscrnphdl,updscrnphdr,drcmsc,longi,xticscb,xticsct,xticsfb,xticsft,PEVT=pevt,DISP=disp,CRND=crnd,CLOG=clog

    if n_params( ) lt 1 then begin
       print,'drcmsc_updmsc20090810,crn,cl0,itl,updmsc0,updmscl,updmscr,'
       print,' updscrnphd0,updscrnphdl,updscrnphdr,drcmsc,...,/pevt,..'
       print,' updated arraies contain the information on left or center'
       return
    end
; Specify the reference time fdoy0 (in days) corresponding to crn and cl0
    scrn0=STRTRIM(crn,2)
    CASE 1 OF
     cl0 lt 10: scl0='00'+STRTRIM(cl0,2)
     cl0 ge 10 AND cl0 lt 100: scl0='0'+STRTRIM(cl0,2)
     else: scl0=STRTRIM(cl0,2)
    ENDCASE
    scrncl0=scrn0+':'+scl0    
    time_crncl0_ctimes,scrncl0,symd_hms0
    year_sec_date,symd_hms0,yyyy0,mm0,dd0,hh0,mi0,ss0
    dmy2doy,yyyy0,mm0,dd0,doy0
    hh0=(hh0+mi0/60.+ss0/3600.)/24.
    fdoy0=doy0+hh0
print,'fdoy0:',fdoy0
; Get updmsct,updscrnphdt 
    sz=SIZE(updmsc0) & xsz=sz(1) & ysz=sz(2) 
    ppd=xsz/360.
    updmsct=fltarr(3*xsz,ysz)
    updscrnphdt=strarr(3*xsz)
    updmsct(0:(xsz-1),*)=updmscl
    updmsct(xsz:(2*xsz-1),*)=updmsc0
    updmsct(2*xsz:(3*xsz-1),*)=updmscr
    updscrnphdt(0:(xsz-1))=updscrnphdl
    updscrnphdt(xsz:(2*xsz-1))=updscrnphd0
    updscrnphdt(2*xsz:(3*xsz-1))=updscrnphdr
print,'updscrnphdr(0:5):',updscrnphdr(0:5)
print,'updscrnphdr(178:183):',updscrnphdr(178:183)
print,'updscrnphdr(355:359):',updscrnphdr(355:359)
; Get phd, slar2 and slar4 for calculating feature's differential rotation 
    zgrid,xsz,ysz,phd,thd,lad,cth,sth
    slar2=fltarr(ysz) & slar4=slar2
    for jj=0,ysz-1 do begin
        slarj=SIN(lad(jj)*!DTOR)
        slar2(jj)=slarj*slarj
        slar4(jj)=slarj^4.
    endfor
; Based on Pevtsov (2000) or Snograss (1983) formulae to get aa,bb,cc and
; the latitude of sunspot where synodic rate=27.2753
    if keyword_set (pevt) then begin
       lat4carringtonrate,xsz,ysz,aa,bb,cc,lat2,ind,/pevt
    endif else begin
       lat4carringtonrate,xsz,ysz,aa,bb,cc,lat2,ind,/snog
    end
    ind0=ind(0) & ind1=ind(1)
; Get fdoyi and the time difference dt with erspect to fdoy0
    if itl eq 180 then is=xsz/2 else is=3*xsz/4
    ie=5*xsz/2-1
    it=ie-is+1
  ; dtt=fltarr(it)
    longi=fltarr(it,ysz) & field=longi
    for ii=is,ie do begin
        updscrncli=updscrnphdt(ii)
        time_crncl0_ctimes,updscrncli,symd_hmsi
        year_sec_date,symd_hmsi,yyyyi,mmi,ddi,hhi,mii,ssi
        dmy2doy,yyyyi,mmi,ddi,doyi
        hhi=(hhi+mii/60.+ssi/3600.)/24.
        fdoyi=doyi+hhi
        CASE 1 OF
           yyyy0 lt yyyyi: begin
               dr=yyyy0 MOD 4
               if dr eq 0 then dt=(366+fdoyi)-fdoy0 else dt=(365+fdoyi)-fdoy0
                           end
           yyyy0 gt yyyyi: begin
               dr=yyyyi MOD 4
               if dr eq 0 then dt=fdoyi-(366+fdoy0) else dt=fdoyi-(365+fdoy0)
                          end
           else: dt=fdoyi-fdoy0
        ENDCASE
  ;     dtt(ii-is)=dt
; Get the latitude-dependent longitude shift relative to CR rotation
        dlogCR1=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
    ;   dlogCR=(360/27.2753)*dt=13.1988*dt
        for jj=0,ysz-1 do begin
            dlogij10=(aa + bb*slar2(jj) + cc*slar4(jj))*dt
            dlogij1=dlogij10 - dlogCR1
            CASE 1 OF
              ii lt xsz: longi(ii-is,jj)=phd(ii)-360-dlogij1
              ii ge 2*xsz: longi(ii-is,jj)=phd(ii-2*xsz)+360-dlogij1
              else: longi(ii-is,jj)=phd(ii-xsz)-dlogij1
            ENDCASE
            field(ii-is,jj)=updmsct(ii,jj)
        endfor
    endfor
; Get differential-rotation-corrected MSC, drcmsc
    drcmsc=fltarr(xsz,ysz)
    for j=0,ysz-1 do begin
        longi_midj0=longi(*,j)
        field_midj0=field(*,j)
        ind=WHERE(field_midj0 ne 0,cc)
        longi_midj=longi_midj0(ind)
        field_midj=field_midj0(ind)
        for i=0,xsz-1 do begin
        ;   indi=where(longi_midj ge i-0.5 AND longi_midj lt i+0.5,ci)
            phi=phd(i)
            indi=where((longi_midj ge phi-1.0) AND (longi_midj lt phi+1.0),ci)
            if ci gt 0 then begin
               br=TOTAL(field_midj(indi))/ci
               drcmsc(i,j)=br
            endif else begin
               drcmsc(i,j)=!values.f_nan
            endelse
        endfor
        indii=where(FINITE(drcmsc(*,j)) eq 1,cii)
 ;if cii lt xsz  then print,'j,cii:',j,cii
        if (cii lt xsz) and (cii gt 2) then begin
           ti=indgen(xsz)
           xi=indii
           yi=drcmsc(indii,j)
           zi=SPLINE(xi,yi,ti)
           drcmsc(*,j)=zi
        endif
    endfor

end