
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

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

      character*6 ex
      character*3 rep
      character*5 sunpo
      character*2 stipox,stipoy,ckp
      character*80 orbidir, bid
      character*64 info(40) 

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

      real (kind=4), dimension(3,1200,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,600/ ! 600*2=1200 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.',           &
      ' ',                                                             &
      ' On rajoute l''orbite de Cluster et les lignes de force ',      &
      ' passant par chaque satellite a l''instant considere.',         &
      ' Utilise  les bibliotheques suivantes:',                        &
      '          comagnetolib.f',                                      &
      '          maglinelib.f',                                        &
      '          geopacklib.f90',                                      &
      '          tsygalib.f90',                                        &
      '          tsauxlib.f90',                                        &
      '          rocotlib.f90',                                        &
      '          rogralib_V9p6.f',                                     &
      '          worldlib.f',                                          &
      '          cluster_orbitlib.f90',                                &
      ' ',                                                             &
      ' 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 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)

      print*
      print*, 'lecture du fichier monde et calcul de l''image projetee'
      print*, '******************************************************'

      call rewormap(1,icowor,3)
      call coworpic(icowor,imawor,idiwor,rep,sunpo)

      print*
      print*, 'lecture des donnees pour le programme plot_magneto'
      print*, '**************************************************'

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


      print*, 'ouverture du fichier graphique/initialisations/header'
      print*, '*****************************************************'

      call dopegra_(2,'coplot_magsphe_cluorbit.ps')
      call dpagfor_('l')

      call dfontyp_('h')
      call ppaghea_('P. Robert, LPP / CNRS  - ')

      print*
      print*, 'plot de la magnetosphere '
      print*, '*************************'

      call 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)

      print*
      print*, 'ajou d''une legende au dessus de la figure'
      print*, '******************************************'

      x=ofx
      y=ofy+sfy+0.3

      call plot_legende(iyear,imon,iday,ih,im,is,ex,rep,ikp,ckp,x,y,tc)

      print*
      print*, 'calcul et trace de la position de 4 sat Cluster'
      print*, '***********************************************'

! *** lecture du directory des fichiers d'orbite


      read(*,100) bid
      print*, 'directory des fichiers d''orbite ?'
      print*, 'ex: /NFS/nas-cluster1/CLUSTER/ORBIT/DATA'
      read(*,100) orbidir
      print  100, orbidir
  100 format(a)

      call coplot_cluster(orbidir,iyear,imon,iday,ih,im,is,rep,ifill)

      print*
      print*, 'calcul et plot des 4 lignes de force passant par les sat'
      print*, '********************************************************'

      print*, 'plot of field lines crossing each spacecraft ? (y/n)'
      read(*,*) kline
      print * , kline

      if(kline.eq.1) then
              call co_clu_pos(orbidir,iyear,imon,iday,ih,im,is,rep)
              call coplot_linesat(ex,ikp,parmod,rep)
                     endif

      print*
      print*, 'fermeture du fichier graphique'
      print*
      
      call dclogra_
      
      print*
      print*, '**********************************************'
      print*, 'coplot_magsphe_cluorbit : normal termination ;'
      print*, 'see file coplot_magsphe_cluorbit.ps'
      print*, '**********************************************'

 
      stop 'coplot_magsphe_cluorbit: NORMAL TERMINATION'
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine co_clu_pos(orbidir,iyear,imon,iday,ih,im,is,rep)
 
! *** calcul l'orbite moyenne et la position de 4 sat Cluster
 
      integer ifc(4)

      real(kind=4), dimension(3,4) :: satpos, satvit
      real(kind=4), dimension(4)   :: satrev

      real(kind=4), dimension(4)  :: sposx,sposy,sposz
      
      common /sat_spos/ sposx,sposy,sposz, spgx,spgy,spgz
 
      character*(*) orbidir
      character*(*) rep
 
      data ifc /11,12,13,14/
 
! *** ouverture des fichiers d'orbite
 
      call openorbi(orbidir,iyear,imon,iday,ifc)
 
! *** calcul des positions de cluster dans le GEI pour l'instant donne
 
      call cluposvit(iyear,imon,iday,ih,im,is,ifc,satpos,satvit,satrev)
 
! *   passage dans le gse ou le gsm
 
      call ctimpar(iyear,imon,iday,ih,im,is)
 
      if(rep(1:3).eq.'GSE'.or.rep(1:3).eq.'gse') call gei_to_gse(satpos)
      if(rep(1:3).eq.'GSM'.or.rep(1:3).eq.'gsm') call gei_to_gsm(satpos)
 
! *** passage en RT
 
      do j=1,4
      sposx(j)=satpos(1,j)/6370.
      sposy(j)=satpos(2,j)/6370.
      sposz(j)=satpos(3,j)/6370.
      enddo
 
! *** position du centre de gravite
 
      spgx= (sposx(1) + sposx(2) +sposx(3) +sposx(4))/4.
      spgy= (sposy(1) + sposy(2) +sposy(3) +sposy(4))/4.
      spgz= (sposz(1) + sposz(2) +sposz(3) +sposz(4))/4.

      call closeorbi(ifc)
 
      return
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine plot_clu_pos(isens)

      real(kind=4), dimension(4)  :: sposx,sposy,sposz
      
      common /sat_spos/ sposx,sposy,sposz, spgx,spgy,spgzspgz

      character*1 coul(4)
 
! *** trace des 4 satellites dasn le plan x-z
 
      sx=0.2
      sy=0.2
 
! *   ordre du trace
 
      coul(1)='n'
      coul(2)='r'
      coul(3)='g'
      coul(4)='b'
 
      print*
      do j=1,4
         sposx(j)= sposx(j)*float(isens)
         print*, 'sat ',j,coul(j),' xyz= ',sposx(j),sposy(j),sposz(j)
      enddo
 
      call sort_arr_xyzc(sposy,sposx,sposz,coul,4,isens)
 
      print*
      print*, 'ordre du trace: ',coul
 
! *   trace des cercles pleins et des contours en noir
 
      call dlinwid_(1)
 
      call dlincol_(coul(1))
      call pfigsym_(sposx(1),sposz(1),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(sposx(1),sposz(1),0,sx,sy,'dode')
 
      call dlincol_(coul(2))
      call pfigsym_(sposx(2),sposz(2),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(sposx(2),sposz(2),0,sx,sy,'dode')
 
      call dlincol_(coul(3))
      call pfigsym_(sposx(3),sposz(3),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(sposx(3),sposz(3),0,sx,sy,'dode')
 
 
      call dlincol_(coul(4))
      call pfigsym_(sposx(4),sposz(4),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(sposx(4),sposz(4),0,sx,sy,'dode')
 
      return
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine coplot_cluster(orbidir,iyear,imon,iday,ih,im,is, &
                                rep,ifill)

      character*(*) orbidir,rep

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

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

!                     -----------------------

! *** calcul et trace de la position des 4 sat Cluster +orbite

      deltat=2.0
      durorb=2.5
      norb=int(durorb*24./deltat)
      print*, 'deltat,durorb,norb=',deltat,durorb,norb

      call cdoyear(iyear,imon,iday,idoty)
      idotyc=idoty
      iyearc=iyear

      call cjd2000(iyear,imon,iday,jd00)

      norbt=0
      dechmag=float(ih) +float(im)/60. +float(is)/3600.

! *   boucle sur le delta T pour le trace des satellites

      do 10 itime=1,norb

! *   sous boucle pour le trace de l'orbite a haute resolution
!      entre 2 satellites

      do 20 istim=1,20

      dech=float(itime-1)*deltat +float(istim-1)*deltat/20.
      print*, 'itime, istim, heure decimale=',itime,istim,dech

! *   reset si > 24h
      nbjou=int(dech/24.)
      dech=dech-24.*float(nbjou)
      jdi= jd00 +nbjou

! *   calcul de la date et de l'heure courante

      call cdatj00(jdi,kyy,kmm,kdd)
      call ctimhou(dech,kh,km,ks)

! *   calcul des positions des 4 sat (charge le common /sat_spos/)

      call co_clu_pos(orbidir,kyy,kmm,kdd,kh,km,ks,rep)

! *   repere pour l'heure ou est calculee la magnetosphere

      if(ifill.eq.1) then
                     call dlincol_('w')
                     else
                     call dlincol_('n')
                     endif

      call dlinwid_(6)
      if(jdi.eq.jd00 .and. abs(dech-dechmag).lt.deltat/20.) then
                  call pfigsym_(spgx*float(isun),spgz,0,0.4,0.4,'dode')
                                                            endif

! *   plot de l'orbite a haute resolution temporelle

      if(norbt.ge.1) then
                     call pfiglin_(pgxp,pgzp, spgx*float(isun),spgz)
                     endif

! *   plot des 4 sat a basse resolution temporelle

      if(istim.eq.1) call plot_clu_pos(isun)
      pgxp=spgx*float(isun)
      pgzp=spgz
      norbt=norbt +1

   20 continue
   10 continue
 
      return
      end 
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
      subroutine plot_legende(iyear,imon,iday,ih,im,is,ex,rep,ikp,ckp, &
                              x,y,tc)

!     ajou d''une legende au dessus de la figure

      character*(*) ex,rep,ckp
      character*6 mod
      character*4 sys
      character*9 cmonth,ctime
      character*84 com

      call cmoncar_(imon,cmonth)
      call ctimcha_(ih,im,is,ctime)
      call cchalen_(cmonth,nm)

      mod=ex
      call upper_case(mod,1)
      call strlen(mod,nc)
      sys=rep
      call upper_case(sys,3)

      print*, 'date/time=',iyear,cmonth(1:nm),iday,ctime(1:5)

     if(ckp.ne.'??') then
                     write(com,200) iyear,cmonth(1:nm),iday,ctime(1:5),&
                                    mod(1:nc), ikp,' : Kp=',ckp,', ',sys
                     else
                     write(com,200) iyear,cmonth(1:nm),iday,ctime(1:5),&
                                    mod(1:nc), ikp,', ',sys
                     endif

  200 format(i4,1x,a,i3,', ',a,' From Tsyganenko model ',a,', ikp=',i1,&
             4a,' system')

      call cchalen_(com,nc)

      print*, com(1:nc)

      call ppagcha_(x,y,-1,tc,tc,0.,com)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

