;
; calculate vector field at a spherical surface of r=r
; MODIFIED:	01APR2004	Xuepu	added /dbl
;
pro zsb_ss,g,h,r,sth,P,dP,cmph,smph,bbr,bbt,bbp,RSS=rss,$
           NMAX=nmax,N0=n0,DBL=dbl
 
    if N_params() LT 1 then begin
       print,'pro sb_ss,g,h,r,pp,dpp,cmph,smph,bbr,bbt,bbp,'
       print,'          rss=,nmax=,n0=,/dbl'
       return
    endif
    if not keyword_set(rss) then rss=2.5
    if not keyword_set(nmax) then begin
       sz1=SIZE(g)
       nmax=sz1(1)-1
    endif
    if keyword_set(dbl) then begin
       zsshf,nmax,rss,r,Hr,dHr,d2Hr,/dbl
       zgfn,cmph,smph,g,h,Gf,dGf,nmax=nmax,/dbl
    endif else begin
       zsshf,nmax,rss,r,Hr,dHr,d2Hr
       zgfn,cmph,smph,g,h,Gf,dGf,nmax=nmax
    endelse
    sz=SIZE(cmph)
    iph=sz(2)
    sz=SIZE(P)
    jth=sz(3)
    if keyword_set(dbl) then begin
       Gfj = dblarr(nmax+1,nmax+1) 
       bbr=dblarr(iph,jth)
    endif else begin
       Gfj = fltarr(nmax+1,nmax+1)
       bbr = fltarr(iph,jth)
    endelse
    dGfj = Gfj
    pi = Gfj
    dpi = pi
    bbt = bbr
    bbp = bbr
    FOR i = 0, iph-1 do begin
        Gfj(*,*) = Gf(0:nmax,0:nmax,i)
        dGfj(*,*) = dGf(0:nmax,0:nmax,i)
        FOR j = 0, jth-1 do begin
            sthi = sth(j)
            pi(*,*) = P(0:nmax,0:nmax,j)
            dpi(*,*) = dP(0:nmax,0:nmax,j)
            if keyword_set(n0) then $
               zbrtp,r,sthi,Hr,dHr,pi,dpi,Gfj,dGfj,Br,Bt,Bp,/n0 else $
               zbrtp,r,sthi,Hr,dHr,pi,dpi,Gfj,dGfj,Br,Bt,Bp
            bbr(i,j) = Br
            bbt(i,j) = Bt
            bbp(i,j) = Bp
        ENDFOR
    ENDFOR

end