
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! XX                                                                  XX
! XX  Complements a la tsygalib pour le calcul de lignes de forces    XX
! XX  et le calcul de magnetopause.                                   XX
! XX  Utilisee pour comagneto et tous les animagneto et anikp         XX
! XX  P. Robert, 1995 - derniere modif Avril 98                       XX
! XX                    revision Octobre 2003 pour bug routine        XX
! XX                    revision Janvier 2008 ajou modules 3D         XX
! XX                    revision Juin    2010 conversion f90          XX
! XX                     et portabilite ifort, gfortran linux/windows XX
! XX                                                                  XX
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine califerjou(ex,ikp,parmod,lines,nblines,nbppl,rep)
     
      common /dimlin/ Nlimax,Nptmax
      
      character (len=*) :: ex,rep
      real (kind=4), dimension(10)       :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl


      real (kind=4), dimension(3,2000)   :: line1       
      real (kind=4), dimension(1000)     :: xl1, yl1, zl1
      real (kind=4), dimension(1000)     :: xl2, yl2, zl2
      real (kind=4), dimension(30)       :: xv 
      
      save line1
      save xl1, yl1, zl1
      save xl2, yl2, zl2
      save xv

!     ******************************************************************
!     calcul des lignes de forces fermees cote jour
!     on part de l'axe X du SM (perpendiculaire au dipole) dans le GSM
!     ******************************************************************

      print*
      print*, 'calcul des lignes de forces fermees cote jour'
      print*, '---------------------------------------------'
      print*
      print*, 'on part de l''axe X du SM (perpendiculaire au dipole) ',&
              'dans le GSM'
      print*

      nblines=0
      ixlast=50

      call  caxvjou(ex,ikp,xv)

      do 20 ix=1,30

      if(ix.gt.ixlast) go to 22

      xism= xv(ix)
   24 continue
      call tsmagsm(xism,0.,0.,xi,yi,zi)
! *** astuce si gse:
      if(rep.eq.'gse') then 
                       call tgsegsm(xi,yi,zi, xii,yii,zii)
                       xi=xii
                       yi=yii
                       zi=zii
      endif
                            
! /// print*, 'ix, xi,yi,zi=',ix,xi,yi,zi

      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)

      if(ite1.eq.0.or.ite2.eq.0) then
! /// print*, 'ligne ouverte cote jour - calcul fin de la limite'
! *** ixlast est calcule pour sortir de la boucle en ix 
                                 ixlast=ix
                                 xism=xism-0.01
                                 go to 24
                                 endif

! *** la ligne calculee est ordonnee puis sauvee dans le tab. lines
!     (soit celle de la boucle, soit la derniere fermee)

      nblines=nblines+1

      if(nblines.gt.Nlimax) then
                            print*,'***warning, nblines max atteint=', &
                                   Nlimax
                            nblines=Nlimax       
                            go to 22
                            endif
                            
! *** calcul de la ligne ordonnee du sud vers le nord

      call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)

! *** remplissage du tableau de toutes les lignes

      call calitout(line1,lines,nblines,nl1,nl2,nbppl)

   20 continue
   22 continue

!   * xism correspond a la derniere ligne fermee
!   * on calcule 5 lignes en dessous, et on choisi comme limite celle
!     qui donne un cusp nord le plus etroit (zmax le plus grand)

      call calmaxtab(zl2,nl2,zmaxp)

! /// print*, 'xism=',xism
! /// print*, 'zmaxp=',zmaxp
! /// print*, ite1,ite2

      do 30 ili5=2,5

      xism=xism-0.01

      call tsmagsm(xism,0.,0.,xi,yi,zi)
      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)

      if(ite1.eq.0.or.ite2.eq.0) go to 30
      call calmaxtab(zl2,nl2,zmax)

! /// print*, 'xism=',xism
! /// print*, 'zmax=',zmax
! /// print*, ite1,ite2

! *** on choisi la meilleure, celle qui donne un zmax le plus grand
!     pour avoir un cusp le plus fin.
!     Pour celle-ci, on met en odre les points du sud vers le nord,
!     le long de B, et on charge le tableau des lignes 
!     lines(2,nbp,num de la ligne)

      if(zmax.gt.zmaxp) then
              call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)
              call calitout(line1,lines,nblines,nl1,nl2,nbppl)
                        zmaxp=zmax
                        endif
   30 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine califernui(ex,ikp,parmod,lines,nblines,nbppl,rep)

      common /dimlin/ Nlimax,Nptmax
      
      character(len=*) :: ex,rep
      real(kind=4), dimension(10)        :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

      real (kind=4), dimension(3,2000)   :: line1       
      real (kind=4), dimension(1000)     :: xl1, yl1, zl1
      real (kind=4), dimension(1000)     :: xl2, yl2, zl2
      
      save line1
      save xl1, yl1, zl1
      save xl2, yl2, zl2

!     ******************************************************************
!     calcul des lignes de forces fermees cote nuit
!     on part de l'axe X du SM (perpendiculaire au dipole) dans le GSM
!     ******************************************************************

      print*
      print*, 'calcul des lignes de forces fermees cote nuit'
      print*, '---------------------------------------------'
      print*
      print*, 'on part de l''axe X du SM (perpendiculaire au dipole) ', &
              'dans le GSM'
      print*

      nblines=0
      ixlast=40

! *   attention, pour le remplissage des liges cote nuit, ixlast
!     doit etre inferieur ou egal a la limite de la figure

      do 10 ix=2,40,2

      if(ix.gt.ixlast) go to 12

      xism=-float(ix)
!  24 continue
      call tsmagsm(xism,0.,0.,xi1,yi1,zi1)
! *** astuce si gse:
      if(rep.eq.'gse') then 
                       call tgsegsm(xi,yi,zi, xii,yii,zii)
                       xi=xii
                       yi=yii
                       zi=zii
      endif

! *   astuce pour prendre en compte le deplacement en Z de la queue

      xi=(xi1+xism)/2.
      yi=0.
      zi=(zi1+0.  )/2.
      
! *   apres 10 RT, on part du Z correspondant a l'etirement maximum

      if(ix.gt.10) then
                   xi=xq -float(ix-8)/2.
                   zi=zq
                   if(abs(xi).gt.float(ixlast)) go to 12
                   endif
      
! ///      print*, 'ix, xi,yi,zi=',ix,xi,yi,zi

      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)

!      if(ite1.eq.0.or.ite2.eq.0) then
! ///        print*, 'ligne ouverte cote nuit - calcul fin de la limite'
!                                 ixlast=ix
!                                 xism=xism+0.05
!                                 go to 24
!                                 endif

      nblines=nblines+1

      if(nblines.gt.Nlimax) then
                            print*,'***warning, nblines max atteint=', &
                                   Nlimax
                            nblines=Nlimax       
                            go to 12
                            endif
                            
! *** calcul de la ligne ordonnee du sud vers le nord

      call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)

! *** calcul du point le plus eloigne ou x change de sens

      call calimin(line1,nl1+nl2, xq,yq,zq)

! *** remplissage du tableau de toutes les lignes

      call calitout(line1,lines,nblines,nl1,nl2,nbppl)

   10 continue
   12 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calidipnui(ex,ikp,parmod,lines,nblines,nbppl)

      common /dimlin/ Nlimax,Nptmax
      
      character (len=*) :: ex
      real (kind=4), dimension(10)       :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

      real (kind=4), dimension(3,2000)   :: line1       
      real (kind=4), dimension(1000)     :: xl1, yl1, zl1
      real (kind=4), dimension(1000)     :: xl2, yl2, zl2
      
      save line1
      save xl1, yl1, zl1
      save xl2, yl2, zl2

!     ******************************************************************
!     calcul des lignes de forces fermees cote nuit
!     on part de l'axe X du SM (perpendiculaire au dipole) dans le GSM
!     ******************************************************************

      print*
      print*, 'calcul des lignes de forces fermees cote nuit'
      print*, '---------------------------------------------'
      print*

      nblines=0
      ixlast=50

      do 10 ix=2,50,2

      if(ix.gt.ixlast) go to 12

      xism=-float(ix)
   24 continue
      call tsmagsm(xism,0.,0.,xi,yi,zi)
! ///      print*, 'ix, xi,yi,zi=',ix,xi,yi,zi

      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)

      if(ite1.eq.0.or.ite2.eq.0) then
! ///        print*, 'ligne ouverte cote nuit - calcul fin de la limite'
                                 ixlast=ix
                                 xism=xism+0.05
                                 go to 24
                                 endif

      nblines=nblines+1

      if(nblines.gt.Nlimax) then
                            print*,'***warning, nblines max atteint=', &
                                   Nlimax
                            nblines=Nlimax       
                            go to 12
                            endif
                            
! *** calcul de la ligne ordonnee du sud vers le nord

      call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)

! *** remplissage du tableau de toutes les lignes

      call calitout(line1,lines,nblines,nl1,nl2,nbppl)

   10 continue
   12 continue


      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine caliouvnor(ex,ikp,parmod,lines,nblines,nbppl,rep)

      common /dimlin/ Nlimax,Nptmax
      
      character (len=*) :: ex,rep
      real (kind=4), dimension(10)        :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

!     ******************************************************************
!     calcul des lignes ouvertes au Nord
!     ******************************************************************

      call caliouvertes(ex,ikp,parmod,lines,nblines,nbppl, 1,rep)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine caliouvsud(ex,ikp,parmod,lines,nblines,nbppl,rep)

      common /dimlin/ Nlimax,Nptmax
      
      character (len=*) :: ex,rep
      real (kind=4), dimension(10)        :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

!     ******************************************************************
!     calcul des lignes ouvertes au Sud
!     ******************************************************************

      call caliouvertes(ex,ikp,parmod,lines,nblines,nbppl,-1,rep)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine caliouvertes(ex,ikp,parmod,lines,nblines,nbppl,ipm,rep)

      common /dimlin/ Nlimax,Nptmax
      
      character (len=*) :: ex,rep
      real (kind=4), dimension(10)       :: parmod

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

      real (kind=4), dimension(3,2000)   :: line1       
      real (kind=4), dimension(1000)     :: xl1, yl1, zl1
      real (kind=4), dimension(1000)     :: xl2, yl2, zl2
      real (kind=4), dimension(30)       :: zv 
      
      save line1
      save xl1, yl1, zl1
      save xl2, yl2, zl2
      save zv
      
      character*4 sudnor(3)
      save sudnor
      data sudnor /'sud ',' ?? ','nord'/

!     ******************************************************************
!     calcul des lignes de forces ouvertes cote nord (ipm=1) ou sud (-1)

!     on part a x=-40 et Z=0 dans le GSM
!     puis on se deplace le long de +/- Z donc perpendiculairement
!     a la direction du soleil.
!     ******************************************************************

      if(ipm.ne.-1.and.ipm.ne.1) stop 'caliouvertes: ipm must be +-1)'

      print*
      print*, 'calcul des lignes de forces ouvertes cote ',sudnor(ipm+2)
      print*, '----------------------------------------------'
      print*
      print*, 'on part a x=-40 et Z=0 dans le GSM', &
             'puis on se deplace le long de Z cote nord et -Z cote sud,'
      print*, 'donc perpendiculairement a la direction du soleil'
      print*

      xi=-40.
      yi=  0.
      izlast=50
      nblines=0

! *** calcul du pas variable le long de Z

      call tsmagsm(xi,0.,0.,xi1,yi1,zi1)
      call cazvjou(ex,ikp,zv)


      do 40 iz=1,30

      if(iz.gt.izlast) go to 42

      zi= float(ipm)*zv(iz) +zi1/6.
! *** astuce si gse:
      if(rep.eq.'gse') then 
                       call tgsegsm(xi,yi,zi, xii,yii,zii)
                       xi=xii
                       yi=yii
                       zi=zii
      endif
   44 continue
! ///      print*, 'iz, xi,yi,zi=',xi,yi,zi

      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)

      if(ite1.eq.1.and.ite2.eq.1) then
! ///                             print*, 'ligne fermee, rejetee'
                                  go to 40
                                  endif

! *** on elimine les lignes trop courtes

      if(nl1.le.2.or.nl2.le.2) go to 40

      if(ite1.eq.0.and.ite2.eq.0) then
! ///  print*, 'ligne sans pied, rejetee - calcul de la dern. avec pied'
! *    si la ligne fait un allez-retour dans la queue, on la garde,
!      et on ne considere pas que c'est la limite; par contre, si
!      un des 2 pieds ouverts est cote jour, on la jette, et on cherche
!      la limite de la precedente.
                       if(xl1(nl1).gt.0. .or. xl2(nl2).gt.0.) then
                                                   izlast=iz
                                                   zi=zi-float(ipm)*0.01
                                                   go to 44
                                                              endif
                                  endif

      nblines=nblines+1
      
      if(nblines.gt.Nlimax) then
                            print*,'***warning, nblines max atteint=', &
                                   Nlimax
                            nblines=Nlimax       
                            go to 42
                            endif

! *** calcul de la ligne ordonnee du sud vers le nord

      call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)

! *** remplissage du tableau de toutes les lignes

      call calitout(line1,lines,nblines,nl1,nl2,nbppl)

   40 continue
   42 continue

! *** amelioration pour le calcul du cusp

!   * zi correspond a la derniere ligne fermee
!   * on calcule 5 lignes en dessous, et on choisi comme limite celle
!     qui donne un cusp le plus etroit (xmax le plus grand)

      call calmaxtab(xl1,nl1,xmaxp)

! /// print*, 'zi=',zi
! /// print*, 'xmaxp=',xmaxp
! /// print*, ite1,ite2

      do 30 ili5=2,5

      zi=zi-float(ipm)*0.01

      call calinca(ex,ikp,parmod, 1.,xi,yi,zi,xl1,yl1,zl1,nl1)
      call calinca(ex,ikp,parmod,-1.,xi,yi,zi,xl2,yl2,zl2,nl2)

      call calinte(xl1(nl1),yl1(nl1),zl1(nl1),ite1)
      call calinte(xl2(nl2),yl2(nl2),zl2(nl2),ite2)
      call calmaxtab(xl1,nl1,xmax)

! /// print*, 'zi=',zi
! /// print*, 'xmax=',xmax
! /// print*, ite1,ite2

! *** on elimine les lignes trop courtes

      if(nl1.le.2.or.nl2.le.2) go to 30

! *** on choisi la meilleure

      if(xmax.gt.xmaxp) then
              call caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)
              call calitout(line1,lines,nblines,nl1,nl2,nbppl)
                        xmaxp=xmax
                        endif
   30 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine caliordo(line1,  xl1,yl1,zl1,nl1,  xl2,yl2,zl2,nl2)

      real (kind=4), dimension(3,2000)   :: line1       
      real (kind=4), dimension(1000)     :: xl1, yl1, zl1
      real (kind=4), dimension(1000)     :: xl2, yl2, zl2
      
!     ******************************************************************
!     concatene et ordonne les 2 lignes en une seule
!     ******************************************************************

      do 10 i=1,nl2
      ii=nl2-i+1
      line1(1,i)=xl2(ii)
      line1(2,i)=yl2(ii)
      line1(3,i)=zl2(ii)
   10 continue

      do 20 j=1,nl1
      jj=nl2+j
      line1(1,jj)=xl1(j)
      line1(2,jj)=yl1(j)
      line1(3,jj)=zl1(j)
   20 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calimin(line1,nbp, xq,yq,zq)

      real (kind=4), dimension(3,2000)   :: line1    
      
!     ******************************************************************
! *** repere l'endroit ou dx change de signe (bout de la queue)
!     ******************************************************************

      dxpre=0.
      
      do 10 i=2,nbp
      dx= line1(1,i) -line1(1,i-1)
      if(dx*dxpre .lt. 0.) then
                           xq=line1(1,i)
                           yq=line1(2,i)
                           zq=line1(3,i)
                           return
                           endif
   10 continue
   
      xq=line1(1,nbp/2)
      yq=line1(2,nbp/2)
      zq=line1(3,nbp/2)
      
      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calitout(line1,lines,nblines,nl1,nl2,nbppl)

      common /dimlin/ Nlimax,Nptmax
      
      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl

      real(kind=4), dimension(3,2000)    :: line1 

!     ******************************************************************
!     charge la ligne courante line1 dans le big tableau lines
!     ******************************************************************

      nbppl(nblines)=nl1+nl2
      jmax=nbppl(nblines)
      
      if(jmax.gt.Nptmax) then
                         print*, '*** warning, jmax=',jmax, &
                                 ' mis a ',Nptmax
                         jmax=Nptmax
                         endif

      do 10 i=1,3
      do 20 j=1,Jmax
      lines(i,j,nblines)=line1(i,j)
   20 continue
   10 continue

      print 100, 'ligne ',nblines,'  nb point:',nbppl(nblines), &
              '  xi,zi=',line1(1,1),line1(3,1),                 &
              '  xf,zf=',line1(1,nbppl(nblines)),line1(3,nbppl(nblines))

  100 format(a,i2,a,i4,a,2f10.4,a,2f10.4)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calimagpau(linjou,linnor,linsud, &
                            nlijou,nlinor,nlisud, &
                            npljou,nplnor,nplsud, &
                            magpau,npau)

      common /dimlin/ Nlimax,Nptmax

      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou, linnor,linsud

      integer,       dimension(Nlimax)          :: npljou, nplnor,nplsud

      real (kind=4), dimension(3,3600) :: magpau
      real (kind=4), dimension(3,3600) :: work
      
      save work

!     ******************************************************************
!     calcul de la ligne representant la magnetopause
!     ******************************************************************

      ns=nplsud(nlisud)
      
      if(ns.gt.2*Nptmax) then
                         print*, '***warning, calimagpau, ns=',ns
                         print*, '    mis a ', 2*Nptmax
                         ns=2*Nptmax
                         endif

      do 10 i=1,3
      do 20 j=1,ns
      work(i,j)=linsud(i,j,nlisud)
   20 continue

      nj=npljou(nlijou)

      if(nj.gt.2*Nptmax) then
                         print*, '***warning, calimagpau, nj=',nj
                         print*, '    mis a ', 2*Nptmax
                         nj=2*Nptmax
                         endif
      do 30 j=1,nj
      j2=ns+j
      jj=nj-j+1
      work(i,j2)=linjou(i,jj,nlijou)
   30 continue

      nn=nplnor(nlinor)

      if(nn.gt.2*Nptmax) then
                         print*, '***warning, calimagpau, nn=',nn
                         print*, '    mis a ', 2*Nptmax
                         nn=2*Nptmax
                         endif
      do 40 j=1,nn
      jj=ns+nj+j
      work(i,jj)=linnor(i,j,nlinor)
   40 continue

   10 continue

      npau=ns+nj+nn

      if(npau.gt.3600) then
                        print*, '***warning, calimagpau, npau=',npau
                        print*, '    mis a 3600'
                        npau=3600
                        endif
                        
! *** inversion pour depart depuis le nord

      do 50 i=1,npau
      ii=npau-i+1
      magpau(1,i)=work(1,ii)
      magpau(2,i)=work(2,ii)
      magpau(3,i)=work(3,ii)
   50 continue

! *** verif que pas de point a l'interieur de la Terre
!     et on ne fait pas descendre la magnetopause a moins de 1/100 RT

      do 60 i=1,npau
      r=sqrt(magpau(1,i)**2 + magpau(2,i)**2 + magpau(3,i)**2)
      if(r.lt.1.01) then
                  magpau(1,i)=magpau(1,i)*1.01/r 
                  magpau(2,i)=magpau(2,i)*1.01/r 
                  magpau(3,i)=magpau(3,i)*1.01/r 
                  endif
   60 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calinte(xf,yf,zf,ite)

!     *****************************************************************
!     calcule si la ligne revient sur terre ou non
!     *****************************************************************


      rf=sqrt(xf**2 + yf**2 + zf**2)
      if(rf.le.1.0001) then
                       ite=1
                       else
                       ite=0
                       endif

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine caxvjou(ex,ikp,xv)

      character(len=*) :: ex

      real(kind=4), dimension(30) :: xv
      real(kind=4), dimension(30) :: xs
      save xs

      data xs/2.0,  4.0,  6.0,  7.5,  8.9,  9.7, 10.5, 11.3, 12.1, 12.9,&
             13.5, 14.0, 14.5, 15.0, 15.5, 16.5, 17.0, 17.5, 18.0, 20.0,&
             22.0, 24.0, 26.0, 28.0, 30.0, 33.0, 35.0, 40.0, 45.0, 50.0/

!     ******************************************************************
!     calcul un pas variable en X pour le depart des lignes cote jour
!     ******************************************************************

      do 10 ix=1,30
      xv(ix)=xs(ix)
   10 continue

      if(ex.ne.'exdish') return

! *** valeur particulieres pour le melange dipole + exshor
!     Dans ce cas, ikp represente le % de melange
!     B= Bigrf + (ikp/100) Bshor
! 
!     On prend donc l'intermediaire entre l'espacement dipole pur
!     et l'espacement exshor

      do 20 ix=1,30
      xvdip=float(ix*2)
      xvsho=xs(ix)
      xv(ix)= (xvdip*(100.-float(ikp)) + xvsho*float(ikp))/100.
   20 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cazvjou(ex,ikp,zv)

      character*(*) ex
      real zv(30)

!     ******************************************************************
!     calcul un pas variable en Z pour le depart des lignes ouvertes
!     ******************************************************************

! *** repartition des Zi sur 30 lignes de forces au maximum
!     le pas se resserre petit a petit

      do 10 iz=1,30
! ///      zv(iz)= 40.*sqrt(sin(3.14158*(iz-1)/60.))
      if(iz.le.7) then
                 zv(iz)=2.*float(iz)
                 else
                 zv(iz)= sqrt(31.4**2 -(35.-float(iz))**2)
                 endif
   10 continue

      if(ex.ne.'exdish') return

! *** valeur particulieres pour le melange dipole + exshor
!     Dans ce cas, ikp represente le % de melange
!     B= Bigrf + (ikp/100) Bshor
! 
!     On prend donc l'intermediaire entre l'espacement dipole pour
!     et l'espacement exshor

      do 20 iz=1,30
      zvdip=float(iz-1)*2.
      zvsho=zv(iz)
      zv(iz)= (zvdip*(100.-float(ikp)) + zvsho*float(ikp))/100.
   20 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cobowsho(bowsho,nbow)

      real (kind=4), dimension(3,3600) :: bowsho

!     ******************************************************************
!     calcul du bowshock  en GSE ou GSM (idem car symetrie autour de x)
!     ******************************************************************

      x1=-40.
      x2= 13.5
      dx=  0.1
      nx= int((x2-x1+0.0001)/dx) +1
      nbow=2*nx

      if(nbow.gt.3600) then
                       print*, 'nbow=',nbow
                       print*, 'GT.3600 *** ERROR***'
                       stop 'error cobowsho'
                       endif

      do 10 ix=1,nx

      x=x1+float(ix-1)*dx
      if(ix.eq.1) x=-75.

      call co_zbowsho(x,z,ierr)

      if(ierr.ne.0) go to 10

      bowsho(1,ix)= x
      bowsho(2,ix)= 0.
      bowsho(3,ix)= z

   10 continue

! *** remplissage pour la partie des z negatifs

      do 20 ix=nx+1,nbow

      bowsho(1,ix)= bowsho(1,nbow-ix+1)
      bowsho(2,ix)= bowsho(2,nbow-ix+1)
      bowsho(3,ix)=-bowsho(3,nbow-ix+1)

   20 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine coflow(gdx,gdz,xflow,zflow,n)

      real (kind=4), dimension(n) :: xflow,zflow
      
!     ******************************************************************
!     calcul du flow  en GSE ou GSM (idem car symetrie autour de x)
!     ******************************************************************

      x1=-40.
      x2= 40.
      x0= 13.5
      dx1=(x2-x1)/float(n)

! *   pas variable: apres le choc, on espace un peu plus les pts

      xb=x1
      dx=dx1
      zp=0.

      do 10 ix=1,n

      xb= xb +dx

! *   on prend le flow comme le bowschock avec une homothetie
!     de centre 0 et de rapport (x0-gdx)/x0

      call co_zbowsho(xb,zb,ierr)

      if(ierr.eq.0.and.abs(zb).gt.abs(gdz)) then
                    xflow(ix)= xb*(x0-gdx)/x0
                    zflow(ix)= sign(zb*(x0-gdx)/x0,gdz)
                    zp=zflow(ix)
                    else
                    xflow(ix)=xb-gdx
                    zflow(ix)=zp
                    dx=dx1*1.5
                    endif

! *   on rajoute une homotethie en x de centre x0-gdx

! /   if(x.lt.x0-gdx) then
! /               xflow(ix)= x -(x0-x)*0.3/abs(gdz)
! /               else
! /               xflow(ix)= x
! /               endif

   10 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine co_zbowsho(x,z,ierr)

!     ******************************************************************
!     calcul du bowshock  en GSE ou GSM (idem car symetrie autour de x)
!     d'apres CLUSTER: Mission payload and supporting activities,
!     ESA SP-1159, page 264
!     ******************************************************************

      c=  56.7
      d= 177.7
      s=  13.5

      p1= (d+s-x)**2
      p2= d*d
      p3=p1/p2-1.

      r2=p3*c*c + x*x
      ierr=0

      if(r2.lt.0.) then
                   ierr=1
                   return
                   endif
! *   z2=r2-x*x
      z2=p3*c*c

      if(z2.lt.0.) then
                   ierr=2
                   return
                   endif

      z=sqrt(z2)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine calmaxtab(zl,nl,zmax)

      real (kind=4), dimension(nl) ::zl

!     ******************************************************************
!     calcul le max tu tableau zl
!     ******************************************************************

      zmax=0.

      do 10 i=1,nl
      zmax=max(zmax,zl(i))
   10 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine verif_lines(lines,nblines,nbppl,factor)

      common /dimlin/ Nlimax,Nptmax

      real (kind=4), dimension(3,Nptmax,Nlimax) :: lines
      integer,       dimension(Nlimax)          :: nbppl


!     *****************************************************************
!     verification des lignes: chaque deltaR ne doit pas etre 
!     plus grand que le deltaR precedent*factor 
!     *****************************************************************

      drp=1.e30

      do 10 nli=1,nblines
      do 20 np =2,nbppl(nli)

      dx2=(lines(1,np,nli)-lines(1,np-1,nli))**2
      dy2=(lines(2,np,nli)-lines(2,np-1,nli))**2
      dz2=(lines(3,np,nli)-lines(3,np-1,nli))**2
      dr=sqrt(dx2 +dy2 +dz2)

      if(dr.gt.drp*factor) then
                           nbppl(nli)=np-1
                           go to 10
                           endif
      drp=dr
   20 continue
   10 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

