;
; WRITTEN:	31Jan2001
;
pro brss_csss_scbf,scbf,obs,ghfcsss,brfrss,IPH=iph,JTH=jth,$
                    NMAX=nmax,APAR=apar,RCP=rcp,RSS=rss
                     
    if N_params() lt 1 then begin
       print,'SYNTAX: brss_csss_scbf,scbf,obs,ghfcsss,brfrss,JPH=,ITH=,NMAX=,APAR='
       print,'                        RCP=,RSS='
       print,"INPUT: scbf (filename for photospheric field), obs ('wso' or 'mdi'),"
       print," ghfcsss, brfrss (filenames for calculated spherical coefficients"
       print," and the radial field at source surface"
       return
    endif

 ;  Default values for parameters

    if not keyword_set(iph) then iph=72
    if not keyword_set(jth) then jth=30
    if not keyword_set(nmax) then nmax=9
    if not keyword_set(rcp) then rcp=2.25
    if not keyword_set(rss) then rss=15.0
    if not keyword_set(apar) then apar=0.2

 ;  Get spherical coefficients at the photosphere based on synoptic maps

    if obs NE 'wso' then gh_mbrsc,scbf,nmax,g,h,/ifl $
                    else gh_wbrsc,scbf,nmax,g,h,/ifl,/cali 
    g=100*g & h=100*h   ; in mT

 ;  Grid spacing

    zgrid,iph,jth,phd,thd,lad,cth,sth

    pdp,cth,Nmax,P,dP

    csmph,phd,nmax,cmph,smph

 ;  Get field distribution on the cusp surface of rcp

    sb_hc,Nmax,apar,rcp,jth-1,iph-1,sth,cmph,smph,g,h,P,dP,brcp,btcp,bpcp

 ;  Get spherical coefficients at the cusp surface based on brcp,btcp,bpcp 

    ghv_csss,apar,Nmax,rcp,rss,jth-1,iph-1,sth,cmph,smph,P,dP,brcp,btcp,bpcp,$
                 ghfcsss,gcsss,hcsss

 ;  Calculate the radial field at the source surface

    sb_csss,Nmax,apar,rcp,rss,rss,phd,thd,g,h,gcsss,hcsss,br

    writefits,brfrss,br

end