
  program comp_field_lines_XZ

  use dim_grid

  character(len=255) :: grid_path,grid_name,flfile
  character(len=1)   :: BJ,pm,ckp

  DIMENSION XXP(1000),YYP(1000),ZZP(1000)
  DIMENSION XXM(1000),YYM(1000),ZZM(1000)
  DIMENSION XXT(2000),YYT(2000),ZZT(2000)
  

! -------------------------------------------------------------------
  common /Rweigh/ Rint
! -------------------------------------------------------------------

  print*,'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
  print*
  print*, 'grid 3D  name? ex: grid_B/grid_Bxyz_160160160_p00_AKp.dat'
  read (*,'(a)') grid_path
  print*, grid_path

  i=index(grid_path,'/',back=.true.)
  if(i>0) then
      grid_name=grid_path(i+1:)
          else
      grid_name=grid_path
  endif
  
  print*, 'grid_name=',trim(grid_name)
  read(grid_name,'(5x,a1,14x,a1,i2,1x,a3)') BJ,pm,idta,ckp
  dta=float(idta)
  if(pm == 'm') dta=-dta

  print*, '------------------------------------------------------------------'
  print*,'grid_name= ',trim(grid_name),'  dta=',dta,' kp=',ckp
  print*, '------------------------------------------------------------------'

  i=len(trim(grid_name))

  flfile='Field_Lines/FL_'//grid_name(6:6)//'XZ'//grid_name(10:i)
  
! ---------------
! read grid data
! ---------------
  
  call read_3D_grid(1,grid_path)
  
! ============================================================== 
  
  call alloc_grid1D()

  print*, 'compute field lines in XZ plane'
  print*, '-------------------------------'


  iyear=2005
  imon=3
  iday=5
  ih=20
  im=36
  is=0

  call cdoyear(iyear,imon,iday,idoty)
  call recalc(iyear,idoty,ih,im,is)

  R0=1.
  RLIM=30. ! enleve car fou la merde
  Rint=1. ! distance max d'interpolation

  print*, 'fill ',trim(flfile)

  open(2,file=flfile)
  
  nltot=0

  print*, 'calcul avec tolerance sur Y'  
  print*
  
  do j=-1,1
     YI=float(j)/2.  
     print*,'Yi=',YI
     
! boucle sur les distances radiales     
  do ir= 2, 20,2
     print*,'    ir=',ir
     
! boucle sur un cercle dans le plan XZ     
  do it=0,360,2

        XI=float(ir)*cos(float(it)*3.14159/180.)
        ZI=float(ir)*sin(float(it)*3.14159/180.)
        
        call give_model_B(XI,YI,ZI,BX,BY,BZ,iout)
        if(iout == 1) cycle
        
        CALL TRACE_PR(XI,YI,ZI, 1.,R0,XXP,YYP,ZZP,LP)      
        CALL TRACE_PR(XI,YI,ZI,-1.,R0,XXM,YYM,ZZM,LM)

        if(LP == 0 .and. LM == 0) cycle

! on reuni les deux lignes en une seule
        
        do l=1,LM
           XXT(LM -l +1)=XXM(l)
           YYT(LM -l +1)=YYM(l)
           ZZT(LM -l +1)=ZZM(l)
        end do
                
        do l=1,LP
           XXT(l+LM)=XXP(l)
           YYT(l+LM)=YYP(l)
           ZZT(l+LM)=ZZP(l)
        end do
        
        LT=LM+LP
!          print*,'XI,k,ZI,iout,LP,LM,LT=',XI,k,ZI,iout,LP,LM,LT   

! on rajoute la valeur du champ en chaque point de la ligne
          
        write(2,100) LT
        do l=1,LT
           call give_model_B(XXT(l),YYT(l),ZZT(l),BX,BY,BZ,iout)
           if(iout == 1) then
                         BX=0.
                         BY=0.
                         BZ=0.
           endif
           write(2,200) XXT(l),YYT(l),ZZT(l),BX,BY,BZ
        end do
        nltot=nltot+1
        
  end do
  end do
  end do
  
  100 format(i6)
  200 format(3F10.4, 3f10.1)
  
  close(2)

  print*
  print*, 'File  ',trim(flfile),' created, nb computed lines: ',nltot

  stop 'comp_field_lines_XZ : normal termination'
  end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
