
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! XX                                                                  XX
! XX  modules pour le trace d'une magnetosphere complete dans le plan XX
! XX  meridien, en gsm ou gse, Soleil a gauche ou a droite.           XX
! XX  La Terre est orientee convenablement a l'heure du calcul.       XX
! XX                                                                  XX
! XX  Utilise la maglinelib pour le calcul des lignes de force        XX
! XX  et la Rogralib pour la visu et la production du PostScript.     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 re_co_magneto(iyear,imon,iday,ih,im,is,ex,ikp, &
                       parmod,rep,sunpo)

      common /sunsid/ isun

      character(len=*) :: ex, sunpo, rep
      real(kind=4), dimension(10) :: parmod

!     ------------------------------------------------------------------
!     lecture des donnees pour le programme co_magneto
!     ------------------------------------------------------------------

      print*

      call redate(iyear,imon,iday)
      call retime(ih,im,is)

      print*, 'iyear,imon,iday=',iyear,imon,iday
      print*, 'ih,im,is=',ih,im,is

      call rap_exikp(ex,ikp)
      call rap_parmod(parmod,10)
      call rap_repgs(rep)
      call rap_sunpo(sunpo)

! *** calcul du isun

      if(sunpo(1:1).eq.'l') then
                            isun=-1
                            else
                            isun= 1
                            endif

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine re_plot_magneto(ifill,isolwi, &
                             irvs1,irvs2,icbs1,icbs2,irbs1,irbs2, &
                             icmp1,icmp2,irmp1,irmp2,ircl1,ircl2, &
                             irps1,irps2,cocli,cooli,codip,corot,coarr)

      character(len=*) :: cocli,cooli,codip,corot,coarr  
      character(len=3) :: eof
      save eof
      
      integer :: ieof=0

!     ------------------------------------------------------------------
!     lecture des options de remplissage ou de plot, et
!     lecture des couleurs pour les lignes de force et le remplissage
!     ------------------------------------------------------------------

      print*, 'remplissage des zones ? (1=oui, 0=non)'
      read *, ifill
      print*, ifill

      print*, 'plot du vent solaire ?'
      read *, isolwi
      print*, isolwi

! *** on ne lit rien si l'eof des couleurs a deja ete rencontre
!     dans un precedent call

      if(ieof.eq.1) return

! *** lecture des couleurs pour la magnetosphere

      read(*,100,end=10) eof
      if(eof.eq.'eof') go to 20

      print*, 'lecture des couleurs...'

      read(*,*)  irvs1,irvs2
      read(*,*)  icbs1,icbs2
      read(*,*)  irbs1,irbs2
      read(*,*)  icmp1,icmp2
      read(*,*)  irmp1,irmp2
      read(*,*)  ircl1,ircl2
      read(*,*)  irps1,irps2

      read(*,'(a1)')  cocli
      read(*,'(a1)')  cooli
      read(*,'(a1)')  codip
      read(*,'(a1)')  corot
      read(*,'(a1)')  coarr

      print*, '    lectcoulmag / couleurs lues:'

      print *, '    irvs1,irvs2=',  irvs1,irvs2
      print *, '    icbs1,icbs2=',  icbs1,icbs2
      print *, '    irbs1,irbs2=',  irbs1,irbs2
      print *, '    icmp1,icmp2=',  icmp1,icmp2
      print *, '    irmp1,irmp2=',  irmp1,irmp2
      print *, '    ircl1,ircl2=',  ircl1,ircl2
      print *, '    irps1,irps2=',  irps1,irps2

      print 100, '     cocli=',  cocli
      print 100, '     cooli=',  cooli
      print 100, '     codip=',  codip
      print 100, '     corot=',  corot
      print 100, '     coarr=',  coarr

      return

   10 continue

      print*, '*** error lectcoulmag: EOF on input file'
      print*, '          colors values are expected, premature ending'
      stop    '*** error lectcoulmag: EOF on input file'

   20 continue

      print*, '    lectcoulmag : eof=',eof,' les couleurs on etees lues'
      print*, '    ou precedemment definies.'
      print*, '    les couleurs sont donc maintenant fixes'
      print*
      ieof=1

  100 format(2a)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine rap_sunpo(sunpo)

      character(len=*) :: sunpo

      common /sunsid/ isun

!     ------------------------------------------------------------------
!     lecture de la direction du Soleil
!     ------------------------------------------------------------------

      print*, 'Sun position? (left or right)'
      read(5,100) sunpo
      print  100, sunpo

      if(sunpo(1:1).ne.'l'.and. &
         sunpo(1:1).ne.'r'   )  then
                                print*, 'Sun position unknown'
                                print*, 'default taken=left'
                                sunpo='left'
                                endif

  100 format(a5)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine rap_dtnbp(deltat,nbp)

!     ------------------------------------------------------------------
!     lecture du deltatT  et du nb de points pour les animations
!     ------------------------------------------------------------------

      print*, 'deltat T (minutes) et nb. de points ?'
      read *,  deltat, nbp
      print*,  deltat, nbp

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine rap_repgs(rep)

      character(len=*) :: rep

!     ------------------------------------------------------------------
!     lecture du repere (gse ou gsm)
!     ------------------------------------------------------------------

      print*, 'repere? (gse or gsm)'
      read  100, rep
      print 100, rep

      if(rep.eq.'gse'.or.rep.eq.'gsm') then
                                       return
                                       else
                                 stop 'que gse ou gsm dans rap_repgs'
                                       endif

  100 format(a)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine co_magneto(iyear,imon,iday,ih,im,is,ex,ikp,parmod,rep,&
                            linjou,linnui,linnor,linsud, &
                            nlijou,nlinui,nlinor,nlisud, &
                            npljou,nplnui,nplnor,nplsud, &
                            magpau,npau,bowsho,nbow,     &
                            dxrep,dyrep,dzrep, rxrep,ryrep,rzrep)

      common /dimlin/ Nlimax,Nptmax

      character(len=*) :: ex, rep

      real(kind=4), dimension(10) :: parmod
      
      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou,linnui, &
                                                   linnor,linsud
      integer,       dimension(Nlimax)    :: npljou,nplnui,nplnor,nplsud

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

!     ------------------------------------------------------------------
!
!     Calcul d'une magnetosphere d'apres le modele de Tsyganenko
!     pour le trace d'une carte meridienne des lignes de force
!     le repere peut etre gse ou gsm
!
!     Calcul selon ecrements 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 rogralib.
!
!     P. Robert, CETP, 1996
!     revu Decembre 2001
!     conversion F90 Juin 2010, P.Ro. LPP
!     ------------------------------------------------------------------

      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 califerjou(ex,ikp,parmod,linjou,nlijou,npljou,rep)
      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

! *** calcul du bowshock dans le gsm

      call cobowsho(bowsho,nbow)

      print*, 'bow shock   : nbp=',nbow      

! ***  lignes pretes pretes pour le trace

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plot_magneto(linjou,linnui,linnor,linsud, &
                              nlijou,nlinui,nlinor,nlisud, &
                              npljou,nplnui,nplnor,nplsud, &
                              magpau,npau,bowsho,nbow,     &
                              imawor,idiwor,               &
                             dxrep,dyrep,dzrep, rxrep,ryrep,rzrep,  &
                             ofx,ofy,sfx,sfy,x1,x2,y1,y2,oarx,oary, &
                             stipox,stipoy,stisil,stisis,grasix,grasiy,&
                             ifill,isolwi, &
                             irvs1,irvs2,icbs1,icbs2,irbs1,irbs2, &
                             icmp1,icmp2,irmp1,irmp2,ircl1,ircl2, &
                             irps1,irps2,cocli,cooli,codip,corot,coarr)

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

      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou,linnui, &
                                                   linnor,linsud
      integer,       dimension(Nlimax)    :: npljou,nplnui,nplnor,nplsud

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

      integer, dimension(idiwor,idiwor) :: imawor

      character(len=*) :: stipox,stipoy
      character(len=*) :: cocli,cooli, codip,corot, coarr

! *   dimension des fleches

      data sarbl,sarbw,sarhl,sarhw /5., 0.35, 1.5, 1.5/
      data iwidBS, iwidMP, iwidLI /7, 7, 2/

!     ------------------------------------------------------------------
!
!     Trace d'une magnetosphere complete calculee par co_magneto.
!     trace d'une carte meridienne des lignes de force
!     le soleil est a gauche (sunpo='left') ou a droite (sunpo='right'),
!     on rajoute la terre et le bow shock
!     le repere peut etre gse ou gsm
!     les couleurs de remplissage sont choisies pour mettre en valeur
!     les zones d'interet
!     le fichier graphique est suppose ouvert.
!
!     P. Robert, CETP, 1996
!     revu Decembre 2001
!     conversion F90 Juin 2010, P.Ro. LPP
!     ------------------------------------------------------------------


      print*, '--- plot de la magnetosphere complete:'
      print*, '    ---------------------------------'

      print*, '--- inversion eventuelle pour avoir le soleil a gauche'

      x1s=  x1
      x2s=  x2
      oarxs=oarx

      if(isun.eq.-1) then
                    
                     x1 = -x2
                     x2 = -x1s 
                     oarx= -oarx
                     aar =  0.
                     print*, '    plot_magneto: Soleil a gauche'
                     else
                     aar =180.
                     print*, '    plot_magneto: Soleil a droite'
                     endif

      print*, '--- definition de la figure'

      call dfontyp_('h')
      call dfigsiz_(sfx,sfy)
      call dfigori_(ofx,ofy)
      call dfigzat_(0.6)

      print*, '--- remplissage du cadre'

      if(ifill.eq.1) then
                     call dlinihs_(irvs1,irvs2)
                     call pfigfra_
                     call pfilzon_
                     endif

      print*, '--- plot des axes'

      call dfontyp_('h')
      call dlinwid_(3)
      call dlincol_('n')

      call dfiglimx(x1,x2)
      call dfiglimy(y1,y2)

      call dstipos_(stipox,stipoy)
      call dstisiz_(stisil,stisis)
      call dgrasiz_(grasix,grasiy)

      call pfiggrax(0.,10.,1.,'(i3)')
      call pfiggray(0.,10.,1.,'(i3)')

      print*, '--- plot de l''image de la terre + un contour circulaire'

      call pfigima_(0.,0.,0,2.,2.,0.,imawor,idiwor,idiwor)

      print*, '--- remplissage du bow shock et plot'

      if(ifill.eq.1) then
                     call dlinwid_(1)
                     call dlinihs_(irbs1,irbs2)
                     call fillmagpau(bowsho,nbow,isun)
                     endif

      call dlinwid_(iwidBS)
      call dlinihs_(icbs1,icbs2)
      call plotmagpau(bowsho,nbow,x1,x2,y1,y2)

      print*, '--- remplissage de la magnetopause'

      if(ifill.eq.1) then
                     call dlinwid_(1)
                     call dlinihs_(irmp1,irmp2)
                     call fillmagpau(magpau,npau,isun)

      print*, '--- remplissage des lignes fermees jour et nuit'

                     call dlinwid_(1)
                     call dlinihs_(ircl1,ircl2)
                     call fillline(linjou,nlijou,npljou,-isun)
                     call fillline(linnui,nlinui,nplnui, isun)

      print*, '--- remplissage de la plasmasphere'

                     call dlinihs_(irps1,irps2)
                     call fillline(linjou,2,npljou,-isun)
                     call fillline(linnui,2,nplnui, isun)
                     endif

      print*, '--- plot des lignes de force ouvertes et fermees'

      call dlinwid_(iwidLI)
      call dlincol_(cocli)
      call plotline(linjou,nlijou,npljou)
      call plotline(linnui,nlinui,nplnui)

      call dlincol_(cooli)
      call plotline(linnor,nlinor,nplnor)
      call plotline(linsud,nlisud,nplsud)

      print*, '--- plot de la magnetopause par dessus la derniere ligne'

      call dlinwid_(iwidMP)
      call dlinihs_(icmp1,icmp2)
      call plotmagpau(magpau,npau,x1,x2,y1,y2)

      print*, '--- plot du contour circulaire de la terre'

      call dlinwid_(1)
      call dlincol_('b')
      call pfigcir_(0.,0.,1.,0.,360.,2.)

      print*, '--- plot du dipole'

      call dlincol_(codip)
      call plotdip(dxrep,dyrep,dzrep,3.0)

      print*, '--- plot de l''axe de rotation de la Terre'

      call dlincol_(corot)
      call plotrot(rxrep,ryrep,rzrep,6.)

      print*, '--- plot des fleches pour la direction du soleil'

      if(oarx.lt.x1.or.oarx.gt.x2) go to 20
      if(oary.lt.y1.or.oary.gt.y2) go to 20

      call gfigsca_(scax,scay)

      call dlincol_(coarr)
      call pfigarr_(oarx, oary,-1,sarbl,sarbw,sarhl,sarhw ,aar)
      call pfilzon_
      call pfigarr_(oarx,-oary,-1,sarbl,sarbw,sarhl,sarhw ,aar)
      call pfilzon_

      call dlincol_('n')
      call dlinwid_(2)
      call pfigarr_(oarx, oary,-1,sarbl,sarbw,sarhl,sarhw ,aar)
      call pfigarr_(oarx,-oary,-1,sarbl,sarbw,sarhl,sarhw ,aar)

   20 continue

      print*, '--- plot du vent solaire'

      if(isolwi.eq.1) call plot_various_flow(y1,y2)

      print*, '--- reimpression des graduations'

      call dfontyp_('h')
      call dlinwid_(3)
      call dlincol_('n')
      call pfiggrax(0.,10.,1.,'(i3)')
      call pfiggray(0.,10.,1.,'(i3)')

      print*, '--- on restore les valeurs initiales des bornes'

      x1=x1s
      x2=x2s
      oarx=oarxs

      print*, '--- plot de la magnetosphere fait.'
      print*, '    ------------------------------'

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plotline(linjou,nlijou,npljou)

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

      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou
      integer,       dimension(Nlimax)          :: npljou

      real (kind=4), dimension(2000) :: x,z ! 2000 max dans Tsyganenko
      save x,z

!     ------------------------------------------------------------------
!     trace d'une famille de ligne de champ
!     ------------------------------------------------------------------

      do 10 i=1,nlijou
      if(i.lt.1 .or. i.gt.Nlimax) stop ' pb i in plotline'
      do 20 j=1,npljou(i)
      if(j.lt.1 .or. j.gt.Nptmax) stop ' pb j in plotline'
      x(j)= float(isun)*linjou(1,j,i)
      z(j)=             linjou(3,j,i)
   20 continue

      call pfigcur_(x,z,max(1,npljou(i)))
   10 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine fillline(linjou,nlijou,npljou,isens)

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

      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou
      integer,       dimension(Nlimax)          :: npljou

      real (kind=4), dimension(2000) :: x,z
      save x,z

!     ------------------------------------------------------------------
!     remplissage d'une famille de ligne de champ
!     ------------------------------------------------------------------

      i=nlijou
      nbp=npljou(i)

      do 10 j=1,nbp
      x(j)= float(isun)*linjou(1,j,i)
      z(j)=             linjou(3,j,i)
   10 continue

      call ccarpol_(x(  1),z(  1),r1,tet1)
      call ccarpol_(x(nbp),z(nbp),r2,tet2)

! *** info sur la zone remplie definie par les lignes de champ

      print*, '    Zone remplie ,teta1,teta2=',tet1,tet2
      print*, '    soit x1,z1=      ', x(1)  ,z(1)
      print*, '         x2,z2=      ', x(nbp),z(nbp)

      if(isens.eq.1) then
                     ntet=int((tet2-tet1)/10.) +1
                     else
                     ntet=int((tet1-tet2+360.)/10.) +1
                     endif

      print*, '         ntet=',ntet

! *** on complete la zone remplie par le contour de la Terre

      do 20 j=1,ntet
      tet= (-float(isens*(j-1))*10. + tet2)*3.14159/180.
      x(nbp+j)= cos(tet)
      z(nbp+j)= sin(tet)
      call ccarpol_(x(nbp+j),z(nbp+j),rj,tetj)
   20 continue

      call pfigcur_(x,z,nbp+ntet)
      call pfilzon_

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plotmagpau(magpau,npau,x1,x2,z1,z2)

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

      common /sunsid/ isun

!     ------------------------------------------------------------------
!     trace de la magnetopause, ou d'un contour non ferme quelconque
!     ------------------------------------------------------------------

      print*, '    plot d''un contour non ferme, n=', npau

      n1=1
      n2=npau

      do 5 j=1,npau
      x(j)= float(isun)*magpau(1,j)
      z(j)=             magpau(3,j)
    5 continue

      do 10 j=1,npau
      if((x(j).gt.x1).and.(x(j).lt.x2).and. &
         (z(j).gt.z1).and.(z(j).lt.z2)) then
                                        n1=j
                                        go to 20
                                        endif
   10 continue
   20 continue

      do 30 j=1,npau
      if(j.le.n1) go to 30
      if((x(j).lt.x1).or.(x(j).gt.x2).or. &
         (z(j).lt.z1).or.(z(j).gt.z2)) then
                                       n2=j-1
                                       go to 40
                                       endif
   30 continue
   40 continue
      print*, '    n1, n2 pris pour le trace=',n1,n2


      if(n1.lt.1) n1=1
      if(n2.lt.1) n2=1
      if(n2.gt.npau) n2=npau

      call pfigcurl(x,z,n1,n2,1)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine fillmagpau(magpau,npau,isens)

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

      common /sunsid/ isun

!     ------------------------------------------------------------------
!     remplissage de la magnetopause, ou d'un contour non ferme 
!     ------------------------------------------------------------------

      print*, '    remplissage d''un contour non ferme'

      do 10 j=1,npau
      x(j)= float(isun)*magpau(1,j)
      z(j)=             magpau(3,j)
   10 continue

      x(1)=x(1)-float(isun)*10.
      x(npau)=x(npau)-float(isun)*10.

      x(npau+1)= x(npau)
      z(npau+1)= 0.

      x(npau+2)= 1.
      z(npau+2)= 0.

      do 20 j=1,37
      tet= float(isens*(j-1))*10.*3.14/180.
      x(npau+2+j)= cos(tet)
      z(npau+2+j)= sin(tet)
   20 continue

      x(npau+40)=x(npau)
      z(npau+40)=0.

      x(npau+41)=x(1)
      z(npau+41)=z(1)

      call pfigcur_(x,z,npau+41)
      call pfilzon_

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

! ///      subroutine fillexbowsho(bowsho,nbow,x1,x2,y1,y2)
! ///c
! ///      real bowsho(3,3600)
! ///c
! ///c------------------------------------------------------------------
! ///c     remplissage de l'exterieur du bow shock
! ///c------------------------------------------------------------------
! ///c
! ///      call plotmagpau(bowsho,nbow,x1,x2,y1,y2)
! ///c
! ///      call pfigpmo_(x1,y1)
! ///      call pfigpmo_(x1,y2)
! ///      call pfilzon_
! ///c
! ///      return
! ///      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plotdip(x,y,z,sca)

      character*2 col

      common /sunsid/ isun

!     ------------------------------------------------------------------
!     trace du dipole terrestre avec une fleche
!     ------------------------------------------------------------------

! *** direction du soleil

      xp=x*float(isun)

! *** axe du dipole

      call dlinwid_(3)
      call glincol_(col)

      xf=xp*sca
      zf= z*sca

      call ccarpol_(xp,z,r,tet)
      call cpolcar_(1.,tet,xs,zs)

      rr=sqrt((zs-zf)**2 +(xs-xf)**2)
      salp=(zs-zf)/rr
      calp=(xs-xf)/rr
      aar=atan2(salp,calp)*180./3.14159

      awid=sca/7.5
      hwid=awid*2.
      hlon=hwid

      if(y*float(isun).lt.0.) then
            call dlincol_(col)
            call pfigarr_( xp, z ,-1, sca,awid,0.,0.,aar+180.)
            call pfilzon_
            call dlincol_('n')
            call pfigarr_( xp, z ,-1, sca,awid,0.,0.,aar+180.)

            call dlincol_(col)
            call pfigarr_(-xs,-zs,-1, sca,awid,hwid,hlon,aar)
            call pfilzon_
            call dlincol_('n')
            call pfigarr_(-xs,-zs,-1, sca,awid,hwid,hlon,aar)
                              else
            call dlincol_(col)
            call pfigarr_( xs, zs,-1, sca,awid,0.,0.,aar+180.)
            call pfilzon_
            call dlincol_('n')
            call pfigarr_( xs, zs,-1, sca,awid,0.,0.,aar+180.)

            call dlincol_(col)
            call pfigarr_(-xp,-z ,-1, sca,awid,hwid,hlon,aar)
            call pfilzon_
            call dlincol_('n')
            call pfigarr_(-xp,-z ,-1, sca,awid,hwid,hlon,aar)
                              endif

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plotrot(x,y,z,sca)

      common /sunsid/ isun

!     ------------------------------------------------------------------
!     trace de l'axe de rotation de la Terre
!     ------------------------------------------------------------------

      xp=x*float(isun)

      call dlinwid_(5)

      call ccarpol_(xp,z,r,tet)
      call cpolcar_(1.,tet,xs,zs)

      xf=xp*sca
      zf= z*sca

! *** axe de rotation

      if(y*float(isun).lt.0.) then
                              call pfiglin_( xp, z , xf, zf)
                              call pfiglin_(-xs,-zs,-xf,-zf)
                              else
                              call pfiglin_( xs, zs, xf, zf)
                              call pfiglin_(-xp,-z ,-xf,-zf)
                              endif

! *** symbole de rotation

      call dlinwid_(3)

      xf=xp*(sca-0.5)
      zf= z*(sca-0.5)

      call pfigros_(xf,zf,0.5,tet)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine convgse(linjou,linnui,linnor,linsud, &
                         nlijou,nlinui,nlinor,nlisud, &
                         npljou,nplnui,nplnor,nplsud, &
                         magpau,npau,bowsho,nbow)

      common /dimlin/ Nlimax,Nptmax

      real (kind=4), dimension(3,Nptmax,Nlimax) :: linjou,linnui, &
                                                   linnor,linsud
      integer,       dimension(Nlimax)          :: npljou,nplnui, &
                                                   nplnor,nplsud

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

!     ------------------------------------------------------------------
!     passage en gse de toutes les lignes de force, magnetopause etc.
!     ------------------------------------------------------------------

      print*
      print*, 'convgse :'
            
      print*, '    passage des lignes jour dans le gse'

      do 10 i=1,nlijou
      do 10 j=1,npljou(nlijou)

      x=linjou(1,j,i)
      y=linjou(2,j,i)
      z=linjou(3,j,i)

      call tgsmgse(x,y,z,linjou(1,j,i),linjou(2,j,i),linjou(3,j,i))

   10 continue

      print*, '    passage des lignes nuit dans le gse'

      do 20 i=1,nlinui
      do 20 j=1,nplnui(nlinui)

      x=linnui(1,j,i)
      y=linnui(2,j,i)
      z=linnui(3,j,i)

      call tgsmgse(x,y,z,linnui(1,j,i),linnui(2,j,i),linnui(3,j,i))

   20 continue

      print*, '    passage des lignes nord dans le gse'

      do 30 i=1,nlinor
      do 30 j=1,nplnor(nlinor)

      x=linnor(1,j,i)
      y=linnor(2,j,i)
      z=linnor(3,j,i)

      call tgsmgse(x,y,z,linnor(1,j,i),linnor(2,j,i),linnor(3,j,i))

   30 continue

      print*, '    passage des lignes sud dans le gse'

      do 40 i=1,nlisud
      do 40 j=1,nplsud(nlisud)

      x=linsud(1,j,i)
      y=linsud(2,j,i)
      z=linsud(3,j,i)

      call tgsmgse(x,y,z,linsud(1,j,i),linsud(2,j,i),linsud(3,j,i))

   40 continue

      print*, '    passage de la magnetopause dans le gse'

      do 50 i=1,npau

      x=magpau(1,i)
      y=magpau(2,i)
      z=magpau(3,i)

      call tgsmgse(x,y,z,magpau(1,i),magpau(2,i),magpau(3,i))

   50 continue

      print*, '    passage du bowshock dans le gse:'

!     En principe inutile 
!     car le bowshock est idem dans le gse (symetrie autour de x)
!     garde par compatibilite

      do 60 i=1,nbow

      x=bowsho(1,i)
      y=bowsho(2,i)
      z=bowsho(3,i)

      call tgsmgse(x,y,z,bowsho(1,i),bowsho(2,i),bowsho(3,i))

   60 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine co_plot_flow(gdx,gdz,zmin,zmax,ts)
 
      common /sunsid/ isun  
 
!     -----------------------------------------------------------------
!     calcule et plot le vent solaire pour une seule ligne
!     P. Robert, CETP, juin 2003
!     -----------------------------------------------------------------
      
      character*1 col
      real (kind=4), dimension(150) :: xflow,zflow
      save xflow,zflow
      
      data n /150/
      save n
 
      call coflow(gdx,gdz,xflow,zflow,n)
      call glincol_(col)
      call dlincol_('y')
 
      do 10 i=1,n
      if(isun.eq.-1) xflow(i)= -xflow(i)
      if(zflow(i).ge.zmax.or.zflow(i).le.zmin) go to 10
 
      call pfigsym_(xflow(i),zflow(i),0,ts,ts,'circ')
      call pfilzon_
   10 continue
 
      call dlincol_(col)
 
      return
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine plot_various_flow(y1,y2)

 
! *** plot du vent solaire
 
      ts=0.12
 
      call co_plot_flow(2., -2.,y1,y2,ts)
      call co_plot_flow(2., -7.,y1,y2,ts)
      call co_plot_flow(2.,-13.,y1,y2,ts)
 
      call co_plot_flow(2.,  2.,y1,y2,ts)
      call co_plot_flow(2.,  7.,y1,y2,ts)
      call co_plot_flow(2., 13.,y1,y2,ts)
 
      return
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine coplot_linesat(ex,ikp,parmod,rep)

      real(kind=4), dimension(10) :: parmod
      real(kind=4), dimension(4)  :: sposx,sposy,sposz

      common /sat_spos/ sposx,sposy,sposz, spgx,spgy,spgz
      common /sunsid/ isun

      character*(*) rep,ex
      character*1 col(4)

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

      data col /'n','r','g','b'/
      save col

      if(rep(1:3).eq.'GSE') rep='gse'
      if(rep(1:3).eq.'GSM') rep='gsm'

      do 10 isat=1,4

      xs=sposx(isat)
      ys=sposy(isat)
      zs=sposz(isat)

      if(rep(1:3).eq.'gsm') then
                   xi=xs
                   yi=ys
                   zi=zs
                            else
                   if(rep(1:3).ne.'gse') then
                                 print*, '*** coplot_linesat: rep=',rep
                                 stop 'coplot_linesat: gsm/gse only'
                                         endif
                   call tgsegsm(xs,ys,zs,xi,yi,zi)
                            endif

      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)

! *   on elimine les lignes qui n'ont aucun pies sur terre

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

      if(ite1.eq.0 .and. ite2.eq.0) then
                                    nl1=1
                                    nl2=1
                                    endif

      if(rep(1:3).eq.'gse') then
                            do i=1,nl1
                               x=xl1(i)
                               y=yl1(i)
                               z=zl1(i)
                               call tgsmgse(x,y,z,xl1(i),yl1(i),zl1(i))
                            enddo

                            do i=1,nl2
                               x=xl2(i)
                               y=yl2(i)
                               z=zl2(i)
                               call tgsmgse(x,y,z,xl2(i),yl2(i),zl2(i))
                            enddo
                            endif

      do i=1,nl1
         xl1(i)=xl1(i)*float(isun)
      enddo

      do i=1,nl2
         xl2(i)=xl2(i)*float(isun)
      enddo

      call dlincol_(col(isat))
      call pfigcur_(xl1,zl1,nl1)
      call pfigcur_(xl2,zl2,nl2)

      print*,'col=',col(isat)
      print*, 'xl1:', (xl1(i),i=20,25)
      print*, 'zl1:', (zl1(i),i=20,25)
   10 continue

      call dlincol_('n')
 
      return
      end 
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine sort_arr_xyzc(x,y,z,c,nbp,isens)
!
      character*(*) c(nbp)
      dimension x(nbp),y(nbp),z(nbp)
!
      real memo1,memo2,memo3
      character*255 memoc
!
!
!     ---------------------------------------------------------------0--
!
!     tri du tableau reel x et rearrangement de 3 tableaux associes
!     y et z sont de type reel, c de type character
!
!     nbp: nb de points des tableaux
!     isens= 1 tri ascendant, isens=-1 tri descendant
!
!     P. Robert, CRPE, 1984, revu juin 2003
!
!     ---------------------------------------------------------------0--
!
!
      ideb=1
      ifin=nbp
!
      if(ideb.gt.0) go to 7
    7 if(isens.ne.1.and.isens.ne.-1) go to 6
      ifinm1=ifin - 1
!
      do 1 io=ideb,ifinm1
      iop1=io + 1
!
      do 2 ki=iop1,ifin
      if(isens) 4,6,3
    3 if(x(io).le.x(ki)) go to 2
!
    5 memo1=x(io)
      x(io)=x(ki)
      x(ki)=memo1
!
      memo2=y(io)
      y(io)=y(ki)
      y(ki)=memo2
!
      memo3=z(io)
      z(io)=z(ki)
      z(ki)=memo3
!
      memoc=c(io)
      c(io)=c(ki)
      c(ki)=memoc
!
      go to 2
!
    4 if(x(io).le.x(ki)) go to 5
    2 continue
    1 continue
!
      return
!
    6 print 100, isens
  100 format('*** sort_arr_xyzc: sens de tri ',i3,' doit etre +/-1')

      stop
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine dlinihs_(ihue22,isat11)

!     ---------------------------------------------------------------+--
! *   Object : define_line_ihsb color  1<ihue22<22, 1<isat11<11
! *   Class  : define modules of Rogralib Software
! *   Author : P. Robert, CRPE, 1992
!     ---------------------------------------------------------------+--

! *** pour compatibilite anterieure a l'introduction des mires

      isat=isat11
      if(isat11.lt. 1) isat=1
      if(isat11.gt.11) isat=11

      sat=1.- float(isat-1)/10.
      call dlinihsb(ihue22,sat,1.)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine dlinihsb(ihue22,sat,bri)

!     ---------------------------------------------------------------+--
! *   Object : define_line_ihsb color  1<ihue22<22, O.<sat<1. O.<bri<1.
! *   Class  : define modules of Rogralib Software
! *   Author : P. Robert, CRPE, 1992
!     ---------------------------------------------------------------+--

      character*64 pal_type
      save pal_type

      real hue22(256)
      real sat22(256)
      real bri22(256)

      save hue22

      data icall /0/
      save icall

      icall =icall+1

! *** valeurs possibles:

      if(ihue22.lt. 1) ihue22=1
      if(ihue22.gt.22) ihue22=22

      pal_type='classic022'
      nlev=256

      if(icall.eq.1) then
                     call dcolmap_(pal_type)
                     call gcolmapv(pal_type,nlev,hue22,sat22,bri22)
                     endif

      hue=hue22(ihue22)

      call dlinhsb_(hue,sat,bri)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine upper_case(string,nc)
!
      character*(*) string
!
!     met la string en majuscule si elle ne l'est pas
!
      do 10 i=1,nc
      ic=ichar(string(i:i))
      if(ic.ge.97.and.ic.le.122) string(i:i)=char(ic-32)
   10 continue

      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine lower_case(string,nc)
!
      character*(*) string
!
!     met la string en minuscule si elle ne l'est pas
!
      do 10 i=1,nc
      ic=ichar(string(i:i))
      if(ic.ge.65.and.ic.le.90) string(i:i)=char(ic+32)
   10 continue

      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX


