;
; WRITTEN:  12APR2006	Xuepu
; MODIFIED: 15JUN2007   Xuepu for DRC 
; MODIFIED: 02July2009  Xuepu for replacing time_crncl0.pro by time_crncl0_ctimes.pro
pro drcmsc_updmsc20090706,crn,cl0,updmsc0,updmscR,updmscL,drcmsc,longi,xticscb,xticsct,xticsfb,xticsft,SNOG=snog,PEVT=pevt,DISP=disp,SDT=sdt,CRND=crnd,CLOG=clog

    if n_params( ) lt 1 then begin
       print,'drcmsc_updmsc20090706,crn,cl0,updmsc0,updmscR,updmscL,'
       print,'       drcmsc,longi,/snog(/pevt),/left,/disp,crnd=,clog='
       return
    end
    mx=100
  ; xz=360 & yz=180 & hxz=xz/4
    sz=SIZE(updmsc0) & xsz=sz(1) & ysz=sz(2) & hxsz=xsz/4
    ppd=xsz/360.
; Get sphd, slar2 and slar4 for calculating feature's differential rotation
    zgrid,xsz,ysz,phd,thd,lad,cth,sth
    sphd=strarr(xsz)
    for i=0,xsz-1 do begin
        CASE 1 OF
         phd(i) lt 10: sphd(i)='00'+strtrim(phd(i),2)
         phd(i) ge 10 and phd(i) lt 100: sphd(i)='0'+strtrim(phd(i),2)
         else: sphd(i)=strtrim(phd(i),2)
        ENDCASE
    endfor
    sphd=STRMID(sphd,0,7)
; print,'sphd:',sphd
    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
; Get the reference time doy0 & hh0 for scrncl0 and getting dt
    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,scrncl0,symd_hms0
    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 Snofrass (1983) formulae and the CR rate 
; (13.1988 deg/day) to get the latitude of sunspot used to obtain 27.2753
    if keyword_set (pevt) then begin
    ;  lat4carringtonrate,lat2,ind,lim,it,jt,aa,bb,cc,/pevt
       lat4carringtonrate,360,180,aa,bb,cc,lat2,ind,/pevt
    endif
    if keyword_set (snog) then begin     
    ;  lat4carringtonrate,lat2,ind,lim,it,jt,aa,bb,cc,/snog
       lat4carringtonrate,360,180,aa,bb,cc,lat2,ind,/snog
    end
    ind0=ind(0) & ind1=ind(1)
    sinl0=SIN(lad(ind0)*!DTOR) & sinl1=SIN(lad(ind1)*!DTOR)
; Get doyi,hhi and the time difference dt between scrnL and scrn0, and
; Find the starting point and get the longitude diferrence for MSCL  
    longi=fltarr(1.5*xz,yz)
    field=longi
    scrnL=strtrim(crn+1,2)
    is=0 & xs=270 
 ;  for ii=iils,xz-1 do begin
    for ii=360-hxz,xz-1 do begin
        scrncli=scrnL+':'+sphd(ii)
      ; time_crncl0,scrncli,symd_hmsi
        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.
        dt=(doyi-doy0)+(hhi-hh0)                 ; dt > 0 in days for MSCL
        for jj=0,yz-1 do begin
        ;   dlogCR=(aa + bb*slar2(ind0) + cc*slar4(ind0))*dt
            dlogCR=13.1988*dt
            dlogij=(aa + bb*slar2(jj) + cc*slar4(jj))*dt - dlogCR
            longi(is,jj)=phd(ii)-dlogij-360
            field(is,jj)=updmscL(ii,jj)
        endfor
        is=is+1
    endfor
    xszl=is
print,'xszl:',xszl
; get the longitude difference from MSC0
    for ii=0,xz-1 do begin
        scrncli=scrn0+':'+sphd(ii)
      ; time_crncl0,scrncli,symd_hmsi
        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.
        dt=(doyi-doy0)+(hhi-hh0)     ; 0 > dt > 0 in days
        for jj=0,yz-1 do begin
        ;   dlogCR=(aa + bb*slar2(ind(0)) + cc*slar4(ind(0)))*dt
            dlogCR=13.1988*dt
            dlogij=(aa + bb*slar2(jj) + cc*slar4(jj))*dt - dlogCR
            longi(is,jj)=phd(ii)-dlogij
            field(is,jj)=updmsc0(ii,jj)
        endfor
        is=is+1
    endfor
    xsz0=is
print,'xsz0:',xsz0
; Find the ending point for MSCR and get the contribution from MSCR 
    scrnR=strtrim(crn-1,2)
;   dttR=fltarr(360)
;   for ii=0,xz-1 do begin
;       scrncli=scrnR+':'+sphd(ii)
;       time_crncl0,scrncli,symd_hmsi
;       yyyyi=FIX(strmid(symd_hmsi,0,4))
;       dmy2doy,yyyyi,FIX(strmid(symd_hmsi,5,2)),$
;               FIX(strmid(symd_hmsi,8,2)),doyi
;       hhi=FIX(strmid(symd_hmsi,11,2))/24.
;       dt=(doyi-doy0)+(hhi-hh0)     ; dt < 0 in days
;       dttR(ii)=dt
;       dlogCR=(aa + bb*slar2(ind(1)) + cc*slar4(ind(1)))*dt
;       dlogij=(aa + bb*slar2(yz-1) + cc*slar4(yz-1))*dt - dlogCR
;       longi_i179=360+phd(ii)-dlogij
;       if longi_i179 gt 360 then goto,lbr
;   endfor
; lbr: iire=ii 
;   for ii=0,iire do begin
    for ii=0,hxz-1 do begin
        scrncli=scrnR+':'+sphd(ii)
      ; time_crncl0,scrncli,symd_hmsi
        time_crncl0_ctimes,scrncli,symd_hmsi
        yyyyi=FIX(strmid(symd_hmsi,0,4))
        dmy2doy,yyyyi,FIX(strmid(symd_hmsi,5,2)),$
                FIX(strmid(symd_hmsi,8,2)),doyi
        hhi=FIX(strmid(symd_hmsi,11,2))/24.
        dt=(doyi-doy0)+(hhi-hh0)     ; dt < 0 in days
        for jj=0,yz-1 do begin
        ;   dlogCR=(aa + bb*slar2(ind(0)) + cc*slar4(ind(0)))*dt
            dlogCR=13.1988*dt
            dlogij=(aa + bb*slar2(jj) + cc*slar4(jj))*dt - dlogCR
            longi(is,jj)=360+phd(ii)-dlogij
            field(is,jj)=updmscR(ii,jj)
        endfor
        is=is+1
    endfor
    xszr=is
    xe=360+phd(ii) & ie=is
; Get differential-rotation-corrected MSC, drcmsc
    drcmsc=fltarr(xz+1,yz)
    for j=0,yz-1 do begin
        longi_midj=longi(*,j)
        field_midj=field(*,j)
        for i=0,xz do begin
        ;   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)
            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)
     ;  indin=where(FINITE(drcmsc(*,j)) ne 1,cin)
if cii lt 361 then print,'j,cii:',j,cii
        if cii lt 361 then begin
           ti=indgen(361)
           xi=indii
           yi=drcmsc(indii,j)
           zi=SPLINE(xi,yi,ti)
           drcmsc(*,j)=zi
        endif
    endfor
    drcmsc=drcmsc(1:xz,*)

end