;
; PURPOSE: Calculate inward field lines using PFSS model
; WRITTEN:  Aug. 1998  Xuepu Zhao
;
pro iflfss_gaha06,ga,ha,iph,jth,iflf,Rss=rss,RBOTTOM=Rbottom,DSTEP=dstep

    if N_params() lt 1 then begin
       print,'pro iflfss_gaha06,ga,ha,iph,jth,iflf,rss=,rbottom=,dstep='
       return
    endif
    if not keyword_set(rss) then rss=2.50
    if not keyword_set(rbottom) then rbottom=1.0
    if not keyword_set(dstep) then dstep=0.02
    zgrid06,iph,jth,phd,thd,lad,cth,sth
    get_lun,unit
    openw,unit,iflf
    for i=0,iph-1 do begin
        phri=phd(i)*!DTOR
        for j=0,jth-1 do begin
            thrj=thd(j)*!DTOR
            if keyword_set(rss) then begin
               if keyword_set(rbottom) then begin
                  if keyword_set(dstep) then begin
                     iflijss_gaha06,ga,ha,phri,thrj,np,fte,xv,yv,zv,rv,tvv,$
                                pv,rss=rss,rbottom=rbottom,dstep=dstep
                  endif else begin
                     iflijss_gaha06,ga,ha,phri,thrj,np,fte,xv,yv,zv,rv,tvv,$
                                pv,rss=rss,rbottom=rbottom
                  endelse
               endif else begin
                  iflijss_gaha06,ga,ha,phri,thrj,np,fte,xv,yv,zv,rv,tvv,$
                                pv,rss=rss
               endelse
            endif else begin
               iflijss_gaha06,ga,ha,phri,thrj,np,fte,xv,yv,zv,rv,tvv,pv
            endelse
            printf,unit,np,fte
            anp=abs(np)
    ;       cnp=n_elements(xv)
    ;       if ABS(cnp - anp) gt 0.5 then print,'cnp,anp,phd(i),thd(j): ',cnp,anp,phd(i),thd(j)
            for k=0,anp-1 do printf,unit,xv(k),yv(k),zv(k),rv(k),tvv(k),pv(k)
        endfor
    endfor
    free_lun,unit

end