;
; calculate vector field at a spherical surface of r=r
;
pro sb_ss06,apar,nmax,rss,r,iph,jth,g,h,bbr,bbt,bbp

    if N_params() LT 1 then begin
       print,'pro sb_ss,apar,nmax,rss,r,iph,jth,g,h,bbr,bbt,bbp'
       return
    endif

    sshf06,apar,nmax,rss,r,Hr,dHr,d2Hr
    zgrid06,iph,jth,phd,thd,lad,cth,sth
    csmph06,phd,nmax,cmph,smph
    pdp06,cth,nmax,P,dP
    gfn06,g,h,cmph,smph,Gf,dGf
    Gfi = fltarr(nmax+1,nmax+1)
    dGfi = Gfi
    pj = Gfi
    dpj = pj
    bbr = fltarr(iph,jth)
    bbt = bbr
    bbp = bbr
    FOR i = 0, iph-1 do begin
        Gfi(*,*) = Gf(*,*,i)
        dGfi(*,*) = dGf(*,*,i)
        FOR j = 0, jth-1 do begin
            sthj = sth(j)
            pj(*,*) = P(*,*,j)
            dpj(*,*) = dP(*,*,j)
            brtp06,r,sthj,Hr,dHr,pj,dpj,Gfi,dGfi,Br,Bt,Bp
            bbr(i,j) = Br
            bbt(i,j) = Bt
            bbp(i,j) = Bp
        ENDFOR
    ENDFOR

end