
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      program co_magpau_3D

      real(kind=4), dimension(10) :: parmod

      character*6 ex
      character*3 rep
      character*5 sunpo
      character*2 stipox,stipoy,ckp
      character*64 info(36) 

!      character(len=1) :: cocli,cooli,codip,corot,coarr

!      real (kind=4), dimension(3,2000,20) :: linjou,linnui,linnor,linsud
!      integer,       dimension(20)        :: npljou,nplnui,nplnor,nplsud

!      real (kind=4), dimension(3,3600) :: magpau,bowsho ! 3600=3*1200

!      integer, dimension(180,360) :: icowor
!      integer, dimension(101,101) :: imawor

      common /sunsid/ isun
      common /dimlin/ Nlimax,Nptmax

      data Nlimax, Nptmax /20,1000/ ! 1000*2=2000 points (// et \//)
      data idiwor /101/
      data info /                                                      &
      ' ',                                                             &
      '--------------------------------------------------------------',&
      ' Trace d''une magnetosphere 2D dans le plan meridien',          &
      ' d''apres le modele de Tsyganenko.',                            &
      ' Le soleil est a gauche ou a droite.',                          &
      ' On rajoute la terre (bien orientee) et le bow shock,',         &
      ' le repere peut etre gse ou gsm.',                              &
      ' Les couleurs de remplissage sont choisies pour mettre en ',    &
      ' valeurles zones d''interet',                                   &
      ' ',                                                             &
      ' Calcul selon des ecrements suivant l axe X le long de SM',     &
      ' pour les lignes fermees',                                      &
      ' (et donc le long d''un axe X perpendiculaire au dipole Z et',  &
      ' avec la direction du soleil dans le plan XZ)',                 &
      ' et le long d''un axe // a Z GSM pour les lignes ouvertes',     &
      ' (et donc perpendiculaire a la direction du soleil)',           &
      ' ',                                                             &
      ' Rappel: SM et GSM ont Y en commun.',                           &
      ' Le depart des lignes de force est dans le plan XZ.',           &
      ' Utilise  les bibliotheques suivantes:',                        &
      '          comagnetolib.f',                                      &
      '          maglinelib.f',                                        &
      '          geopacklib.f90',                                      &
      '          tsygalib.f90',                                        &
      '          tsauxlib.f90',                                        &
      '          rocotlib.f90',                                        &
      '          rogralib_V9p6.f',                                     &
      '          worldlib.f',                                          &
      ' ',                                                             &
      ' P. Robert, CETP, 1996',                                        &
      '    revu Decembre 2001',                                        &
      '    revu Janvier  2010 migration SUN -> linux & PC compatible ',&
      '    revu Juin     2010 pour compatibilite avec les Roprocs',    &
      '                       et les differents modeles de Tsyganenko',&
      ' ',                                                             &
      '--------------------------------------------------------------' /
     
! *** parametres de la figure
 
      data ofx,ofy / 1.1, 1.8/
      data sfx,sfy /26. ,16. /
 
!      data x1,x2   /-40., 25./
!      data y1,y2   /-20., 20./
 
! *   option grande echelle (attention garder le rapport 65/40)
 
      data x1,x2   /-80., 50./
      data y1,y2   /-40., 40./
 
      data oarx,oary / 18., 10. /
      data tc        /  0.30    /
 
      data stipox,stipoy /'oi','oi'/
      data stisil,stisis /0.30,0.15/
      data grasix,grasiy /0.40,0.40/
!                          _____________________
                          
      
      print '(a)', info
 
      print*
      print*, 'lecture des donnees pour le programme co_magneto'
      print*, '************************************************'
 
      call re_co_magneto(iyear,imon,iday,ih,im,is,ex,ikpl, &
                         parmod,rep,sunpo)
 
      print*
      print*, 'calcul du ikp pour la date et l''heure consideree'
      print*, '*************************************************'
 
      call cokp(1,iyear,imon,iday,ih,rkp)
      call convkptsy(ex,rkp,ikp)
      call convrkp(rkp,ckp)
 
! *   si pas de Kp dans la base, on prend la valeur lue
 
      if(ikp.lt.0) then
                   print*, 'base des kp innaccessible'
                   ikp=ikpl
                   print*, 'kp pris=kp lu =',ikp
                   else
                   print*, 'rkp,ckp,ikp=',rkp,ckp,ikp
                   endif

      print*
      print*, 'calcul et trace de la magnetosphere complete'
      print*, '********************************************'
 

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

! *** calcul du dipole dans le gsm ou le gse

      call gdipdir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo)
      
      if(rep.eq.'gsm') then
                       call tgeogsm(dxgeo,dygeo,dzgeo,dxrep,dyrep,dzrep)
                       else
                       call tgeogse(dxgeo,dygeo,dzgeo,dxrep,dyrep,dzrep)
                       endif

! *** calcul de l'axe de rotation de la Terre dans le gsm ou le gse

      if(rep.eq.'gsm') then
                       call tgeogsm(0.,0.,1.,rxrep,ryrep,rzrep)
                       else
                       call tgeogse(0.,0.,1.,rxrep,ryrep,rzrep)
                       endif

! *** calcules des familles de lignes de forces dans le gsm

      call califerjou3D(ex,ikp,parmod)
!      call califernui(ex,ikp,parmod,linnui,nlinui,nplnui,rep)
!      call caliouvnor(ex,ikp,parmod,linnor,nlinor,nplnor,rep)
!      call caliouvsud(ex,ikp,parmod,linsud,nlisud,nplsud,rep)

      print*
      print*, 'co_plot_magneto:'
      print*, '--------------- '
      print*
!      print*, 'nb de lignes jour: ',nlijou
!      print*, 'nb de lignes nuit: ',nlinui
!      print*, 'nb de lignes nord: ',nlinor
!      print*, 'nb de lignes sud : ',nlisud

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

      factor=4.

! /// print*
! /// print*, 'verification de l''espacement des points dans les lignes'
! /// print*, 'limitation a dr/drpre <',factor
! /// print*

! /// call verif_lines(linjou,nlijou,nbljou,factor)
! /// call verif_lines(linnui,nlinui,nblnui,factor)
! /// call verif_lines(linnor,nlinor,nblnor,factor)
! /// call verif_lines(linsud,nlisud,nblsud,factor)

! *** calcul de la magnetopause dans le gsm

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

!      print*
!      print*, 'magnetopause: nbp=',npau



!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      stop 'co_magpau_3D: NORMAL TERMINATION'
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine califerjou3D(ex,ikp,parmod)
      
      character (len=*) :: ex
      real (kind=4), dimension(10)       :: parmod


      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

      print*
      print*, 'calcul de la magnetopause 3D'
      print*, '****************************'

!     fichier resultat

      open(1,file='magpau.resu')

! position du cusp dans le plan xz du gsm (coplot_magsphe.out)

!cote jour
!ligne  6  nb point:  95  xi,zi=    0.3706    0.9273  xf,zf=    0.0199   -0.9963
!ligne  7  nb point: 112  xi,zi=    0.3342    0.9356  xf,zf=    0.0050   -0.9946
!cote nord
!ligne 11  nb point:  92  xi,zi=    0.2441    0.9650  xf,zf=  -77.8018   22.1864
!ligne 12  nb point: 110  xi,zi=    0.3222    0.9445  xf,zf=  -77.8012   22.7808
!cote sud 
!ligne 13  nb point: 114  xi,zi=  -77.8036  -21.3201  xf,zf=   -0.0432   -0.9912

! au nord:

      xij=0.3342
      zij=0.9356

      xin=0.3222
      zin=0.9445

      xij=0.3706
      zij=0.9273

      xin=0.2441
      zin=0.9650

      rcusp=sqrt((xin-xij)**2 + (zin-zij)**2)/2.
      xcusp=(xij+xin)/2.
      zcusp=(zij+zin)/2.

      pisd=acos(-1.)/180.

! boucle en phi

      iphi1=0
      iphi2=360
      iphi3=10

      nphi=(iphi2-iphi1)/iphi3 +1

      write(1,200) nphi

      do iphi=iphi1,iphi2,iphi3
      phi=float(iphi)
      write(*,*) 'phi=',phi

      xp=rcusp*cos(phi*pisd)
      yp=rcusp*sin(phi*pisd)
      zp= 0.

      cosa=0.5*(xin-xij)/rcusp
      sina=0.5*(zin-zij)/rcusp

      xi= xcusp + xp*cosa +zp*sina
      yi=  yp
      zi= zcusp - xp*sina + zp*cosa

                            
! /// 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)
                 
! *** calcul de la ligne ordonnee du sud vers le nord

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

! *   ecriture de la ligne dans le fichier resu

      nli=nl1+nl2
      write(1,200) nli

      do i=1,nli
        write(1,100) line1(1,i),line1(2,i),line1(3,i)
      enddo

      enddo ! sur iphi

      close(1)

  100 format(3E14.5)
  200 format(i5)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

