;
; WRITTEN:  15JUN2007   Xuepu Zhao for DRC 
; MODIFIED: 24MAY2009   Xuepu Zhao
;
pro drcsc_updsc20090810b,crn,cl0,itl,updsc0,updscr,updscl,updscrncld0,updscrncldr,updscrncldl,drcsc,longi,xticscb,xticsct,xticsfb,xticsft,PEVT=pevt,DISP=disp,CRND=crnd,CLOG=clog

    if n_params( ) lt 1 then begin
       print,'ldrcsc_updsc2009080l0b,crn,cl0,itl,updsc0,updscr,'
       print,'                       drcsc,longi,...,/pevt'
       return
    end
; Get phd, slar2 and slar4 for calculating feature's differential rotation 
    sz=SIZE(updsc0) & xsz=sz(1) & ysz=sz(2) & hxsz=xsz/4
    ppd=xsz/360.
    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
; Specify the reference time fdoy0 (in days) for 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
; 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, the time difference dt with erspect to fdoy0, and
;  longitude shift and field redistribution
    longi=fltarr(3*xsz,ysz) & field=longi
    for ii=hxsz,xsz-1 do begin
        scrncli=updscrncldl(ii)
        time_crncl0_ctimes,scrncli,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
        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
            longi(ii-hxsz,jj)=phd(ii)-360-dlogij1
            field(ii-hxsz,jj)=updsc0(ii,jj)
        endfor
    endfor
    i1=ii
    for ii=0,xsz-1 do begin
        scrncli=updscrncld0(ii)
        time_crncl0_ctimes,scrncli,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
        dlogCR1=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        for jj=0,ysz-1 do begin
            dlogij10=(aa + bb*slar2(jj) + cc*slar4(jj))*dt
            dlogij1=dlogij10 - dlogCR1
            longi(ii+i1,jj)=phd(ii)-dlogij1
            field(ii+i1,jj)=updsc0(ii,jj)
        endfor
    endfor
    i2=ii+i1
; Find the ending point for MSCR and get the contribution from MSCR
    dtt =fltarr(xsz)
    for ii=0,xsz-1 do begin
        scrncli=updscrncldr(ii)
        time_crncl0_ctimes,scrncli,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)=dt
        dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        dlogij20=(aa + bb*slar2(ysz-1) + cc*slar4(ysz-1))*dt
        dlogij2=dlogij20 - dlogCR 
        if ABS(dlogij20)  ge ABS(dlogCR1) then goto,lbr
    endfor
  lbr: iire=ii 
    for ii=0,iire-1 do begin
        dt=dtt(ii)
        dlogCR=(aa + bb*slar2(ind1) + cc*slar4(ind1))*dt
        for jj=0,ysz-1 do begin
            dlogij0=(aa + bb*slar2(jj) + cc*slar4(jj))*dt
            dlogij=dlogij0 - dlogCR
            longi(ii+i2,jj)=360+phd(ii)-dlogij
            field(ii+i2,jj)=updscR(ii,jj)
        endfor
    endfor
; Get differential-rotation-corrected MSC, drcsc
    drcsc=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
            phi=phd(i)
        ;   indi=where(longi_midj ge i-0.5 AND longi_midj lt i+0.5,ci)
          ; indi=where(longi_midj ge i-1.0 AND longi_midj lt i+1.0,ci)
            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
               drcsc(i,j)=br
            endif else begin
               drcsc(i,j)=!values.f_nan
            endelse
        endfor
        indii=where(FINITE(drcsc(*,j)) eq 1,cii)
     ;  indin=where(FINITE(drcsc(*,j)) ne 1,cin)
; if cii lt 361 then print,'j,cii:',j,cii
        if cii lt 360 then begin
           ti=indgen(360)
           xi=indii
           yi=drcsc(indii,j)
           zi=SPLINE(xi,yi,ti)
           drcsc(*,j)=zi
        endif
    endfor

end