
*JMT(
*	Ai modifie les appels a VCFT2 en appels a NVCFT2
*JMT)
*
*
******************************
*
*****  VERSION NE NECESSITANT QUE DES MODIFICATIONS DANS LE DATA
*
*      A PART LORSQUE RESTART=.TRUE.
*
*             ATTENTION--ATTENTION--ATTENTION
*
*            1024 EST UN CHIFFRE CLE DE VECTORISATION
*
*              ++++++ IL APPARAIT DANS LPLXD !! DON'T FORGET ++++++++
*
******************************
*         VERSION COMPATIBLE AVEC BORDEAUX
*     POUR L'INSTANT ENLEVER : MODIF==> XXXXVP
*
*   VERSION PERMETTANT DE SELECTIONNER L'ANGLE EN DEGRES
*
*  LE MUR REFLECHISSANT A ETE FAIT
*  LE COURANT APPLIQUE A ETE MODIFIER
*  TEMPLOC A {T{ MODIFIER EN CONS{QUENCE
*  LA SUBROUTINE OBLIC A {T{ FAITE
*  LA D{TERMINATION DE L'ANGLE THETA A {T{ FAITE
*
*
*
*
*        BOITE POUR RUN RESISTIF OU NON RESISTIF ==> VOIR FICHIER DATA
*        BOITE POUR RUN RESISTIF OU NON RESISTIF ==> VOIR FICHIER DATA
*        BOITE POUR RUN RESISTIF OU NON RESISTIF ==> VOIR FICHIER DATA
*        BOITE POUR RUN RESISTIF OU NON RESISTIF ==> VOIR FICHIER DATA
*
*       DEPEND DU PARAMETRE RESIST DANS LE DATA
*
*
*      **************  CODE  VECTORISE  ****************
*
*    LE NOMBRE DE CHAQUE ESPECE DE PARTICULES DOIT ETRE DIVISIBLE
*
*    PAR L'ENTIER  *************   2048   *******************
*
*    ON DOIT RENTRER LA VALEUR DE LPLX DANS DENSIT
*
*  RUN ENTIEREMENT VECTORISE PAR JEAN-MARIE TEULER
*
*    *                                            *
*    **********************************************
*
*   ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*   ||                                                            ||
*   ||                      ATTENTION                             ||
*   ||                                                            ||
*   ||  DU FAIT DE CETTE PROCEDURE TOUS LES FICHIERS ONT DU       ||
*   ||                                                            ||
*   ||  ETRE DECLARES DANS LE JCL (A CAUSE DU COMPAT1 NECESSAIRE  ||
*   ||                                                            ||
*   ||  AU CLOSE( ) ) DONC LE VJOBID NE SERT PLUS A RIEN          ||
*   ||                                                            ||
*   ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*
*
*
*  POSSIBILITE D'AGGRANDIR LA BOITE EN COURS DE RUN , LA SEULE LIMITE
*
*  EST QU'IL FAUT TOUJOURS GARDER LA MEME DENSITE PAR CELLULE (PAS
*  ENCORE FAIS)
*
*  LE GRDATA A ETE ETOFFE DE PLUS J'AI RAJOUTER GRDRES POUR LES
*
*  MESSAGES IMPORTANT AU COUR DU RUN ET J'AI ENLEVER LA SORTIE FINALE
*
*  DES RESTART (ELLE ETAIT EN PERRUQUE PAR RAPPORT A LA PROCEDURE
*
*  NOUVELLE "CF" ABSYS )
*
*
************
*
** PERMET DE DESSINER LES COMPOSANTES D'UN RUN OBLIQUE
*
************
*
*  LE PARAMETRE MUS REPRENSENTE LE COEFFICIENT (ENTIER) DETERMINANT
*
*  LE NOMBRE DE PARTICULES DESSINEES (1/MUS)
*
*
**
**  TOUTES LES SORTIES POUR LE POSTPROCESSING CE FAIT EN
**  SIMPLE PRECISION , DE PLUS LA SORTIE DU TEMPS NT EST
**  CONVERTIE EN REEL
**
**
**  LE CALCUL DES DENSITE ET DES COURANTS CE FAIT GRACE A UN APPEL
**
**  ON A RAJOUTER UN DIAGNOSTIQUE PERMETTANT D'ENLEVER LA CONTRIBUTION
**  DU PISTON DANS LE CALCUL DE L'ENERGIE DU CHAMP ELECTROSTATIQUE
**  SUBROUTINE PISTON ET FELECT
**  ON UTILISE UN TEMPLOC 1-D VECTORISE
**  LA BOITE EST MISE A 2048X128 ET LPLX A 1025
**  LE NOMBRE TOTALES DE PARTICULES EST ALORS DE 8192
*
* PROCEDURE PERMETTANT DE FAIRE DES RESTARTS DANS DES FICHIERS
* DIFFERENTS |
* AJONCTION DU DIAGNOSTIQUE TEMPLOC
* QUAND ON CHANGE LE NOMBRE DE BOITES NE PAS OUBLIER CELL1 ET CELL2
* TRACE L'ONDE DE DILATATION POUR LA COMPOSANTE DU CHAMP BOY DOMINANT
* JE SORS EN BLOC LES FICHIERS "PARTICULES" SUR BANDE
* AMELIORATION PAR LA VECTORISATION
*
*   1.JE MET LES BOUCLES DO DANS LE BON SENS (LA OU IL N'Y A AUCUN
*      RISQUE
*   2.JE TRAITE LES IF EN SCALAIRE LORSQU'ILS NE SONT UTILISES
*     QUE PEU DE FOIS
*   3.JE NE BALAYE PAS TOUTES LES PARTICULES EN MEME TEMPS
*     POUR LE CALCUL DES DENSITES ET POUR LE DIAGNOSTIQUE TEMPLC
*  NO2 DOIT ETRE DIVISIBLE PAR 2048
*  |4.VECTORISATION DE LA BOUCLE DES DENSITES
*  |  POUR CE FAIRE ON CREE --128 TABLEAUX DE 2048X32
*  |            (WCUX,WCUY,WCUZ,WG1X)
*  |                        --LE PARAMETRE LWH=128
*  |PAS FAIT CAR DEMANDE TROP DE MEMOIRE
* J'AI MIS TOUTES LES DIMENSIONS ET VARIABLES EN PARAMETER
* POUR CELA IL SUFFIT DE CHERCHER LE MOT CHANGING
*
* J'AI ENLEVER TURESS , SMOOTH ET TOUT LE PROGRAMME CONREC
*
*
*
* JE DISTINGUE DE D2V4 PAR L'UTILISATION DE C06FRE DE NAG POUR
* LES T.F. COMPLEXES -> COMPLEXES (AU LIEU DE DVCFT2)
*
* FAIT :
*   ADAPTATIONS DES ENTR{ES-SORTIES
*   (CALL LINK,BUFFER OUT)
* - ADATATIONS DES APPELS GRAPHIQUES @ L'EXCEPTION DE TDPLOT ET DES
*   TRAC{S DE CONTOURS
* - EMULATION DE SECOND
* - EMULATION DE AMINMX
* - TRONCATURE @ 6 CARACT}RES DES VARIABLES
*   TOUS LES IOXM ET IOYM ONT {T{ TRANSFORM{S EN IXM ET IYM
*   TOUS LES LABEL ONT {T{ TRANSFORM{S EN LABL
*   ATTENTION : ZCITOA -> ZCITO ET PICTUREC -> PICTUR
* - AI RENOMMM{ VCFT2PRP EN VCFT2P
* - AI D{PLI{ LES LIGNES REPLI{ES PAR LE TRANSMIT
* - AI D{PLAC{ DES CONTINUE DE LA COLONNE 6 VERS LA COLONNE 7
*   IDEM POUR DES FORMAT
* - AI CHANG{ LES CARACT}RES " EN '
* - AI MIS LES APPELS @ TDPLOT EN COMMENTAIRES ET AI SORTI TDPLOT
* - IDEM POUR CONSCM
* - AI MIS IMPLICIT DOUBLE PRECISION PARTOUT
* - AI AJOUT{ L'INITIALISATION GRAPHIQUE
* - AI {MUL{ ACONTR
* - AI {MUL{ VCFT2
* - AI SUPPRIM{ LE SOUS-PROGRAMME ARROW (JAMAIS APPEL{)
* - AI INTRODUIT USERID ET JOBID POUR LES NOMS DE FICHIERS
*   (PARAM}TRES @ INITIALISER)
* - AI {MUL{ RANORM ET RANF
* - AI AJOUT{ UN ARGUMENT @ GRAF1 POUR BALAYER SEULEMENT UNE PARTIE DES
*   INDICES. MODIFICATION CORRESPONDANTE DANS MONLI3
* - AI AJOUT{ TROIS ARGUEMENTS @ TURSCM (DEVENU TURESS)
*   INDICE MINIMUM X, INDICE MINIMUM Y, VALEUR DE COUPURE
*
C
C***********    2-1/2 D E.M.CODE **********
C
C++
C++ WARNING *****  WARNING  *****  WARNING
C++
C++  BE CAREFUL  FOR SPECULAR REFLEXION
C++  FOR PLOT PHASE SPACE
C++  BECAUSE WE WANT THE PERPENDICULAR COMPONENT
C++
C++  WHEN YOU CHANGE BYO
C++
C++ WARNING  *****  WARNING  *****  WARNING
C++
C+
C+ PASSAGE EN DOUBLE PR{CISION
C+
C--      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER ABSYS
      INTEGER FACT1,FACT2

        INCLUDE  "PARAM"
*
* ********************
*  POUR MEMOIRE IL EST RAPPELE QUE NTISPA EST LE PRODUIT ENTRE
*  (NBRE DE PARTICULES SAUVER)*(NBRE DE NSS PAR TRANCHE DE CALCUL)
*
*  SI AUCUNE MEMORISATION EST EFFECTUE ALORS MEMOP=.FALSE.
*
*    ==>  NTISPA=4*NGAP*(NEND/NSS)
*
*
*  CHANGING AREA
*
      PARAMETER (NTISPA=5000)
      LOGICAL MEMOP
      DATA    MEMOP/.FALSE./
*
      REAL XELSP(NTISPA),YELSP(NTISPA)
      REAL XIOSP(NTISPA),YIOSP(NTISPA)
      DIMENSION NTSP(NTISPA)
      PARAMETER(FIELD1=FIELD3+2)
      PARAMETER(FACT1=FIELD3/2+1,FACT2=FIELD2/2+1)
      PARAMETER(PLOP=PART/MUS)
C  COMMON STATEMENTS
C
      COMMON/PART/X(PART),Y(PART)
      COMMON/VITE/PX(PART),PY(PART),PZ(PART),PSA(PART)
      COMMON/FORM/SKX(FACT1,FACT2),SKY(FACT1,FACT2),FBX(FACT1,FACT2)
     1  ,FBY(FACT1,FACT2)
      DIMENSION BRX(FIELD1,FIELD2),BRY(FIELD1,FIELD2),BRZ(FIELD1,FIELD2)
      COMMON/XX/T
      COMMON/BOUND/BNCXR,BNCXL,BNCXL2,BNCXR2,ANCYL,ANCYR,BOUNDY,BOUNDX
      COMMON/PERIOD/IR(FIELD3),IL(FIELD3),JR(FIELD2),JL(FIELD2)
C++
C++ TO PLOT ONLY HALF POINTS
C++
      COMMON/PLOT/XG(PLOP),YG(PLOP),PXG(PLOP),PYG(PLOP)
     & ,PZG(PLOP)
      REAL XG,YG,PXG,PYG,PZG
C++
C++ TO PLOT IN THE PERPENDICULAR PLANE IN OBLIQUE PROPAGATION
C++
      COMMON/OBL/PER1(PLOP),PER2(PLOP),PAR(PLOP)
      REAL PER1,PER2,PAR
C++
C++  FOR TEMPLC
C++
      COMMON/ANGLE/SINT,COST
C++
C++ FOR PLOT ONLY FEW PARTICLES
C++
      COMMON/PLOPAR/N8O
C++
C++  FOR SPECULARLY REFLECTION
C++
      COMMON/RESI/RESIST
      LOGICAL RESIST
C++
C++ DIMENSION TO PLOT LOCAL TEMPERATURE
C++
      DIMENSION TPLCE1(PLOH),TPLCE2(PLOH),TPALCE(PLOH)
      DIMENSION TPLCI1(PLOH),TPLCI2(PLOH),TPALCI(PLOH)
C++
C+++    POUR MEMORISER LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART
          DIMENSION NOMPAE(NUMPAE),NOMPAI(NUMPAI),
     &  TRAJEL(NUMPAE,NEMDI,13),TRAJIO(NUMPAI,NEMDI,NUMFIC)
C
C DIMENSION STATEMENTS
C
      DIMENSION CUX(FIELD1,FIELD2),CUY(FIELD1,FIELD2),CUZ(FIELD1,FIELD2)
      DIMENSION G1XION(FIELD1,FIELD2),G1X(FIELD1,FIELD2)
      DIMENSION ELX(FIELD1,FIELD2),ELY(FIELD1,FIELD2)
      DIMENSION ETX(FIELD1,FIELD2),ETY(FIELD1,FIELD2),ETZ(FIELD1,FIELD2)
      DIMENSION ERX(FIELD1,FIELD2),ERY(FIELD1,FIELD2),ERZ(FIELD1,FIELD2)
      DIMENSION ERX1(FIELD1,FIELD2),ERY1(FIELD1,FIELD2)
      DIMENSION BTX(FIELD1,FIELD2),BTY(FIELD1,FIELD2),BTZ(FIELD1,FIELD2)
      DIMENSION ELXEN(FIELD1,FIELD2)
      DIMENSION FEX(FACT1,FACT2),FEY(FACT1,FACT2)
      DIMENSION CKX(FACT1,FACT2),CKY(FACT1,FACT2)
      DIMENSION CKZ(FACT1,FACT2),CKC(FACT1,FACT2)
      DIMENSION REX(FACT1,FACT2),REY(FACT1,FACT2),AK(FACT1,FACT2)
      DIMENSION LNCEL(3)
      DIMENSION EFL(PLOH),EFT(PLOH),BFT(PLOH),EFIEL(PLOH)
      DIMENSION EKINE(PLOH),EKINI(PLOH),EKIN(PLOH),ETOT(PLOH)
      DIMENSION TEMPEX(PLOH),TEMPEY(PLOH),TEMPEZ(PLOH)
      DIMENSION TEMPIX(PLOH),TEMPIY(PLOH),TEMPIZ(PLOH)
      DIMENSION CUEX(PLOH),CUEY(PLOH),CUEZ(PLOH)
      DIMENSION CUIX(PLOH),CUIY(PLOH),CUIZ(PLOH)
      DIMENSION RATIO(PLOH)
C++
C++  DIMENSION STATEMENT FOR ENERGY
C++
      DIMENSION EFLX(PLOH),EFLY(PLOH)
      DIMENSION EFTX(PLOH),EFTY(PLOH),EFTZ(PLOH)
      DIMENSION BFTX(PLOH),BFTY(PLOH),BFTZ(PLOH)
      DIMENSION EFLXEN(PLOH)
C
C AUXILIARY QUANTITIES
C
      COMMON/CONST/NO,NCX,NCY,LVGC
      COMMON/CONST1/RMASS,RTE
      DIMENSION SK(FACT1,FACT2)
C  VARIABLES FOR PLOTTING THE POTENTIAL
      DIMENSION PK(FACT1,FACT2)
      DIMENSION POT(FIELD1,FIELD2)
C  CORRELATION VARIABLES
C     DIMENSION VEL(192),VEBTX(192),VEBTY(192),VEBTZ(192)
C     DIMENSION VG1X(192)
C  AVERAGE ELECTRIC FIELD VARIABLES
      DIMENSION ERXAV(FIELD1,FIELD2),ERYAV(FIELD1,FIELD2)
      DIMENSION ERXXY(FIELD2),ERYXY(FIELD2)
      DIMENSION ERXYX(FIELD3),ERYYX(FIELD3)
C PARTICLE ORBIT VARIABLES
C     DIMENSION STORE(256)
C     DIMENSION XBF(PART)
C  VARIABLES NEEDED FOR TIME AVERAGING
      DIMENSION CUXION(FIELD1,FIELD2)
      DIMENSION G1XAV(FIELD1,FIELD2),G1XIAV(FIELD1,FIELD2)
      DIMENSION CUXAV(FIELD1,FIELD2),CUXIAV(FIELD1,FIELD2)
      DIMENSION DUM(FIELD1,FIELD2)
      DIMENSION CUXX(FIELD3),CUXY(FIELD2)
      DIMENSION G1XX(FIELD3),G1XY(FIELD2)
      DIMENSION G1XIOX(FIELD3),G1XIOY(FIELD2)
      DIMENSION CUXIOX(FIELD3),CUXIOY(FIELD2)
C++
C++  MASKING FUNCTION
C++
      DIMENSION FM(FIELD3)
C++
C++ DIMENSION FOR APPLIQUED CURRENT
C++
      DIMENSION CUZAP(FIELD3,FIELD2)
C  CUZ CURRENT
      DIMENSION CUZAV(FIELD1,FIELD2),CUZIAV(FIELD1,FIELD2)
C--   DIMENSION CUZTAV(FIELD1,FIELD2),CUZION(FIELD1,FIELD2)
      DIMENSION CUZION(FIELD1,FIELD2)
      DIMENSION CUZ1(FIELD1,FIELD2)
      DIMENSION CUZY(FIELD3),CUZIOY(FIELD3),CUZTY(FIELD3)
      DIMENSION CUZX(FIELD3),CUZIOX(FIELD3),CUZTX(FIELD3)
C++
C++ DIMENSION FOR CUY
      DIMENSION CUYION(FIELD1,FIELD2),CUYX(FIELD3),CUYIOX(FIELD3)
C  VARIABLES USED FOR PARTIAL SPATIAL AVERAGING OVER NX AND NY
      DIMENSION G1XYM(FIELD2),G1XIYM(FIELD2)
      DIMENSION CUXYM(FIELD2),CUXIYM(FIELD2)
      DIMENSION CUZYM(FIELD2),CUZIYM(FIELD2),CUZTYM(FIELD2)
      DIMENSION G1XXM(FIELD3),G1XIXM(FIELD3)
      DIMENSION CUXXM(FIELD3),CUXIXM(FIELD3)
      DIMENSION CUZXM(FIELD3),CUZIXM(FIELD3),CUZTXM(FIELD3)
C  TEMPORARY STORE
      DIMENSION CUZTT(FIELD3),CUZTAP(FIELD3)
C TIME HISTORY OF THE APPLIED CURRENT
      DIMENSION CUZAPP(PLOH)
C++
C++  DIMENSION FOR AVERAGE FIELDS ON Y VERSUS X
C++
      DIMENSION ELXYX(FIELD3),ELYYX(FIELD3),POTYX(FIELD3)
      DIMENSION ETXYX(FIELD3),ETYYX(FIELD3),ETZYX(FIELD3)
      DIMENSION BTXYX(FIELD3),BTYYX(FIELD3),BTZYX(FIELD3)

      LOGICAL MOVION,START,SWITCH,HISTO
      CHARACTER*20 LABL2
      CHARACTER*20 LABL3,LABL4,LABL5
      CHARACTER*20 LABL6,LABL7,LABL8
      CHARACTER*20 LABL9,LABL10,LABL11
      CHARACTER*20 LABL12,LABL13,LABL14
*
      CHARACTER*6 JOBID
      COMMON/JOB/JOBID

C
C  SET UP INPUT DATA
C
      NAMELIST/OPTION/NOX,NOY,INCX,INCY,AX,AY,NEND,DT,C
     1  ,START,MOVION,HISTO
      NAMELIST/IONS/RMASS,RTEMP
      NAMELIST/VELO/VTX,VTY,VTZ,VDX,VDY,VDZ
      NAMELIST/MAGF/BO,ANGTHE
      NAMELIST/NAME/JOBID,RESIST
      NAMELIST/APPCUR/CURREN,WM,WN,OMEGA,NPE
      NAMELIST/NEW/NCONT,NCONTB,NLL,NSS,NTIM,NX,NY,LPLX
C+
C+ D{CLARATIONS SUPPL{MENTAIRES
C+    ... ESPACE DE TRAVAIL POUR LES FFT
      real table,work
      common/fft/ table(100+2*(field1+field2)),
     &            work(512*max(field1,field2))

C+    ... POUR LE STOCKAGE DES NOMBRES AL{ATOIRES NORMAUX
      REAL RANORM(PART*3)
C+    ... POUR LA GESTION DU TRAC{
      COMMON/CADNUM/NUMPAG,NUMVUE,NTYPE
      DATA NUMPAG/0/
      DATA NUMVUE/0/
C+    INITIALISATION LA VARIABLE GRAPHIQUE NTYPE
C+    (SI NTYPE = 0, 1 VUE PAR PAGE, SI NTYPE=3, 4 VUES PAR PAGE)
      DATA NTYPE /1/
*JMT(
      REAL TIME1,TIME2,SECOND
      EXTERNAL SECOND
      INTEGER LASF(13)
      DATA LASF/13*1/
*JMT)
C
C     MEASURE THE COMPUTING TIME REQUIRED.
      TIME1 = SECOND()
C+
*     OUVERTURE DE GKS
      CALL OPNGKS
      CALL GSASF(LASF)
      CALL GSTXFP(1,2)

*     4 VUES PAR PAGE
      NTYPE = 0
*     INTIALISATION PAGE ET VUE
      NUMPAG = 0
      NUMVUE = 0

      OPEN(25,FORM='UNFORMATTED')
      OPEN(26,FORM='UNFORMATTED')
      OPEN(22,FORM='UNFORMATTED')
      OPEN(23,FORM='UNFORMATTED')
      OPEN(24,FORM='UNFORMATTED')
C++
C+++    POUR MEMORISER LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART
      IF (NUMPAE  .NE.  1) THEN
      OPEN(50,FORM='FORMATTED')
      OPEN(51,FORM='UNFORMATTED')
      READ(50,*) NOMPAE      
      ENDIF
      IF (NUMPAI .NE. 1) THEN
      OPEN(60,FORM='FORMATTED')
      OPEN(61,FORM='UNFORMATTED')
      READ(60,*) NOMPAI
      ENDIF 
C++
C++  OPEN DISK FILES TO STORE THE FIELDS COMPONENTS
C++
      OPEN(30,FORM='UNFORMATTED')
      OPEN(31,FORM='UNFORMATTED')
      OPEN(32,FORM='UNFORMATTED')
      OPEN(33,FORM='UNFORMATTED')
      OPEN(34,FORM='UNFORMATTED')
      OPEN(35,FORM='UNFORMATTED')
      OPEN(36,FORM='UNFORMATTED')
      OPEN(37,FORM='UNFORMATTED')
      OPEN(39,FORM='UNFORMATTED')

      READ(5,OPTION)
      READ(5,IONS)
      READ(5,VELO)
      READ(5,MAGF)
      READ(5,NAME)
      READ(5,APPCUR)
      READ(5,NEW)
      REWIND 5
C++
C++ PLOT DATA ON PLOT FILE
C++
      CALL GRDATA
C
C SOME INITIAL CONSTANTS
      LABL2=JOBID
      LABL3='TOT CUZ'
      LABL4='X-Y AXIS'
      LABL5='APP CUR'
      LABL6='ELX'
      LABL7='ELY'
      LABL8='ERX'
      LABL9='ERY'
      LABL10='ERZ'
      LABL11='BRX'
      LABL12='BRY'
      LABL13='BRZ'
      LABL14='POT'
C++
C++ FOR TEMPLOC SUBROUTINE
C++
      CELEL=0
      CELIO=0
C
C  PUSHING
      RTE=C*C
      CIN=1./RTE
      VTH=1./C
      VIT=1./VTH
      VC=DT*0.5
C  ELECTRONS AND IONS
      BETE=0.5*DT
      ALFE=DT*VTH*0.5
      BETI=-BETE*RMASS
      ALFI=-ALFE*RMASS
      SQV=SQRT(RMASS/RTEMP)
      VTXI=VTX*SQV
      VTYI=VTY*SQV
      VTZI=VTZ*SQV
C  PARTICLES AND GRIDS
      NO=NOX*NOY
      NO2=2*NO
      NO1=NO+1
      NCX=2**INCX
      ANCX=NCX
C++
C++  TO PLOT ONLY FEW POINTS
      N4O=PLOP
      N8O=N4O/2
      N8O1=N8O+1
C++
C++  AVERAGE DENSITIE
C++
      NCY=2**INCY
      DA=(NOX*NOY)/(NCX*NCY)
      DA10=DA
      NENDO=0
C++
C++ PARAMETER FOR FFT2
C++
      N2=NCX+2
      NCXX=NCX/2
      NCXH1=NCX/2+1
      ANCY=NCY
      NCYY=NCY/2
      NCYH1=NCY/2+1
      NCYH2=NCYY+2
      NCYP2=NCY+2
      NCXY=NCX*NCY
C--   DDELX=ANCX/FLOAT(NOX)
      DDELY=ANCY/FLOAT(NOY)
C--   BOUNDX=ANCX-0.00002000
      BOUNDY=ANCY-0.00002000
      ANCXR=ANCX-0.50001000
      ANCYR=ANCY-0.50001000
C--   ANCXL=-0.49999000
      ANCYL=-0.49999000
C--   XMIN=ANCXL
C--   XMAX=ANCXR
      YMIN=ANCYL
      YMAX=ANCYR
      NOY2=NOY/2
C++
C++  PARTICLES BOUNDARIES IN X
C++
      BNCXL=LPLX-1.49999000
      BNCXR=NCX-1.50001000
      XMAX=BNCXR
      XMIN=BNCXL
      BNCXL2=2.*BNCXL
      BNCXR2=2.*BNCXR
      DDELX=(BNCXR-BNCXL+0.00002)/FLOAT(NOX-1)
CCC        DDELX=(BNCXR-BNCXL)/FLOAT(NOX-1)
C++
C++ PARAMETER FOR THE CHARGE SPACE CORRECTION
C++  
            LXLPOT=LPLX-1
            LXRPOT=NCX

C++
C++  AVERAGE PARAMETER
C++
        MOX=NCX-LPLX
        MOX1=MOX+1
C  FORM FACTORS CONSTANTS
         AEC=FLOAT(NO)
C++
C++  PARAMETER FOR BOUNDED PLASMA
C++
      FGK=(BNCXR-BNCXL+0.00002)*ANCY/AEC
      FCK=FGK*DT
      EKX=6.2831853/ANCX
      EKY=6.2831853/ANCY
      EXX=EKX*AX
      EYY=EKY*AY
      RKX=VIT*DT*EKX
      RKY=VIT*DT*EKY
C  CORRELATION MODES
      NCXSP=12
      NCYSP=8
      NCXYSP=NCXSP*NCYSP
C++
C++ MASKING FUNCTION
C++
      ALPLX=LPLX-1.
      ALPLX=ALPLX*ALPLX
      ALPLX=1./ALPLX
      DO 51 I=1,NCX
      IF (I.LT.LPLX) FM(I)=-ALPLX*(I-LPLX)*(I-LPLX)+1.
      IF (I.GE.LPLX) FM(I)=1.
   51 CONTINUE
C++
C++  EXTERNAL CURRENT PARAMETER
C++
      IF (RESIST) THEN
      CURREN=-ABS(CURREN)
      ELSE
      CURREN=ABS(CURREN)
      ENDIF
      NPE2=NPE/2
      CUREN2=CURREN/2.
      NCUO=LPLX-7
      NCUE=LPLX-2
C++  DEFINE THE NEW PARAMETERS FOR THE TUBE CURRENT
C++
      RAYON=10.
      POSIT=LPLX-RAYON-NCUE
      PI=4.*ATAN(1.D0)
      TWOPI=2.*PI
C++
C++ FIND ANGLE BETWEEN BO AND X,Z
C++
      THET=(PI*ANGTHE)/180.D0
      COST=COS(THET)
      SINT=SIN(THET)
      IF (RESIST) THEN
      BXO=BO*COST
      BZO=BO*SINT
      BYO=0.D0
      ELSE
      BXO=BO*COST
      BYO=BO*SINT
      BZO=0.D0
      ENDIF
      PRINT*,'COST=',COST
      PRINT*,'SINT=',SINT
      PRINT*,'BXO=',BXO
      PRINT*,'BYO=',BYO
      PRINT*,'BZO=',BZO
*
*    INITIALISE L'INCREMENT DES TABLEAUX POUR LE SUIVIT TEMPOREL
*    DE CERTAINES PARTICULES
*
      K1SPA=1
      K2SPA=1
C
C  PREPARE FFT2
C
      LNCEL(1)=INCX-1
      LNCEL(2)=INCY
      LNCEL(3)=0
      LVGC=NCX+2
C++
C+    INITIALISATION DES FFT COMPLEXES 
      IFSET = 0
      CALL NVCFT2(CUX,N2,NCX,NCY,IFSET,TABLE,WORK)
C
C  SET UP PERIODIC COEFFICIENT ARRAYS
C
      DO 951 M=1,NCX
      IR(M)=M+1
  951 IL(M)=M-1
      IR(NCX)=1
      IL(1)=NCX
      DO 952 N=1,NCY
      JR(N)=N+1
  952 JL(N)=N-1
      JR(NCY)=1
      JL(1)=NCY
C
C FORM FACTORS FOR CHARGE DENSITY ,FIELDS AND ENERGY
C
      DO 4119 N=1,NCYH1
      DO 4119 M=1,NCXH1
      FEX(M,N)=0.
      FEY(M,N)=0.
      FBX(M,N)=0.
      FBY(M,N)=0.
      REX(M,N)=0.
      REY(M,N)=0.
      CKX(M,N)=0.
      CKY(M,N)=0.
      CKZ(M,N)=0.
      CKC(M,N)=0.
       PK(M,N)=0.
      AK(M,N)=0.
 4119 CONTINUE
      DO 420 M=1,NCXH1
      DO 420 N=1,NCYH1
      MM=M-1
      NN=N-1
*VOCL STMT,IF(1)
      IF(M.EQ.1.AND.N.EQ.1) MM=NCXX
      WKX=MM*EKX
      WKY=NN*EKY
      SQX=MM*EXX
      SQY=NN*EYY
      SQAA=SQX*SQX+SQY*SQY
*VOCL STMT,IF(1)
      IF(SQAA.GT.20.) GO TO 420
      WKS=WKX*WKX+WKY*WKY
      VARNEW=SQRT(RTE*WKS)
      EXK=EXP(-SQAA)
      FOM=FGK*EXK
      FEX(M,N)=FOM*WKX/WKS
      FEY(M,N)=FOM*WKY/WKS
C NEXT LINE IS FOR POTENTIAL PLOTS
      PK(M,N)=FGK*EXK/WKS
      AK(M,N)=AEC/EXK
*VOCL STMT,IF(1)
      IF(VARNEW.GT.18.0) GO TO 420
      FBX(M,N)=VTH*FEX(M,N)
      FBY(M,N)=VTH*FEY(M,N)
      REX(M,N)=RKX*MM
      REY(M,N)=RKY*NN
      CKZ(M,N)=FCK*EXK
      CKX(M,N)=CKZ(M,N)*WKY*WKY/WKS
      CKY(M,N)=CKZ(M,N)*WKX*WKX/WKS
      CKC(M,N)=CKZ(M,N)*WKX*WKY/WKS
 420  CONTINUE
      FEX(1,1)=0.
      FEY(1,1)=0.
      FBX(1,1)=0.
      FBY(1,1)=0.
      REX(1,1)=0.
      REY(1,1)=0.
      CKX(1,1)=0.
C--   CKY(1,1)=0.
C--   CKZ(1,1)=0.
      CKC(1,1)=0.
      FEX(1,NCYH1)=0.
      FEY(1,NCYH1)=0.
      FBX(1,NCYH1)=0.
      FBY(1,NCYH1)=0.
      REX(1,NCYH1)=0.
      REY(1,NCYH1)=0.
      CKX(1,NCYH1)=0.
      CKY(1,NCYH1)=0.
      CKZ(1,NCYH1)=0.
      CKC(1,NCYH1)=0.
C++
C++  FREE SPACE CORRECTION FOR TANGENTIEL FIELDS
C++
      CKY(1,1)=FCK
      CKZ(1,1)=FCK
      PK(1,1)=0.
      PK(1,NCYH1)=0.
C
C  FORM FACTORS FOR EXTERNAL B FIELD AND VECTOR POTENTIAL
C
      DO 104 N=1,NCYH1
      DO 104 M=1,NCXH1
      SKX(M,N)=0.
      SKY(M,N)=0.
      SK(M,N)=0.
  104 CONTINUE
      DO 106 M=1,NCXH1
      DO 106 N=1,NCYH1
      MM=M-1
      NN=N-1
*VOCL STMT,IF(1)
      IF(M.EQ.1.AND.N.EQ.1) MM=NCXX
      WNX=MM*EKX
      WNY=NN*EKY
      WNS=WNX*WNX+WNY*WNY
      SKX(M,N)=-WNX/WNS
      SKY(M,N)=-WNY/WNS
      SQX=MM*EXX
      SQY=NN*EYY
      SQAA=SQX*SQX+SQY*SQY
*VOCL STMT,IF(1)
      IF(SQAA.GT.10.) GO TO 106
      EXK=EXP(-SQAA)
      SK(M,N)=EXK
  106 CONTINUE
      SKX(1,1)=0.
      SKX(1,NCYH1)=0.
      SKY(1,1)=0.
      SKY(1,NCYH1)=0.
      SK(1,1)=0.
      SK(1,NCYH1)=0.
C
C  INITIALIZATION OF TIME COUNTERS
C
      NC=0
      NT=0
      T=0.
C
C  VALUE OF THE LAST TIME STEP
C
      NSTOP=NEND*DT
C
C  INITIALIZE  ENERGIES ( TIME )
C
C  INITIALIZE FIELD ENERGIES
      DO 4121 N=1,NSTOP
      EFL(N)=0.
      EFT(N)=0.
      BFT(N)=0.
      EFIEL(N)=0.
C  INITIALIZE KINETIC AND TOTAL ENERGIES
      EKINE(N)=0.
      EKINI(N)=0.
      EKIN(N)=0.
      ETOT(N)=0.
C  INITIALIZE TEMPERATURES
      TEMPEX(N)=0.
      TEMPEY(N)=0.
      TEMPEZ(N)=0.
      TEMPIX(N)=0.
      TEMPIY(N)=0.
      TEMPIZ(N)=0.
C  INITIALIZE CURRENTS
      CUEX(N)=0.
      CUEY(N)=0.
      CUEZ(N)=0.
      CUIX(N)=0.
      CUIY(N)=0.
      CUIZ(N)=0.
      RATIO(N)=0.
C  INITIALIZE EXTERNAL CUZ CURRENT
      CUZAPP(N)=0.
C++
C++ INITIALIZE ENERGY
C++
      EFLX(N)=0.
      EFLY(N)=0.
      EFTX(N)=0.
      EFTY(N)=0.
      EFTZ(N)=0.
      BFTX(N)=0.
      BFTY(N)=0.
      BFTZ(N)=0.
 4121 CONTINUE
C
C  INITIALIZE ITERATION ELECTRIC FIELDS ( SPACE )
C
      DO 4118 J=1,NCY
      DO 4118 I=1,NCX
      ETX(I,J)=0.
      ETY(I,J)=0.
      ETZ(I,J)=0.
      ERX(I,J)=0.
      ERY(I,J)=0.
      ERZ(I,J)=0.
      ERX1(I,J)=0.
      ERY1(I,J)=0.
      BTX(I,J)=0.
      BTY(I,J)=0.
      BTZ(I,J)=0.
      BRX(I,J)=0.
      BRY(I,J)=0.
      BRZ(I,J)=0.
      G1XAV(I,J)=0.
      G1XIAV(I,J)=0.
      ERXAV(I,J)=0.
      ERYAV(I,J)=0.
 4118 CONTINUE
C++
C++ FOR THE VECTORISATION ON VP200
C++
C++ I3 AND I4 FOR THE ELECTRONS , OTHERS FOR IONS
C++
      I3=1
      I4=NO/2048
      I5=I4+1
      I6=NO2/2048

C  RESTART PROCEDURE
C
      IF(START) GO TO 1000

      OPEN(12,FORM='UNFORMATTED')
      OPEN(8,FORM='UNFORMATTED')
      OPEN(9,FORM='UNFORMATTED')
      OPEN(15,FORM='UNFORMATTED')
      OPEN(16,FORM='UNFORMATTED')
      OPEN(17,FORM='UNFORMATTED')

      READ(12) NENDO,X,Y
      READ(12) PX,PY,PZ
      READ(12) PSA
      NT=NENDO
      NEND=NEND+NENDO
      READ(8) NENDI,ETX,ETY,ETZ
      READ(8) BTX,BTY,BTZ
      READ(9) NCLAST,EFL,EFT,BFT,EFIEL,EKINE
     1  ,EKINI,EKIN,ETOT
     2  ,EFLX,EFLXEN,EFLY,EFTX,EFTY,EFTZ
     3  ,BFTX,BFTY,BFTZ
     4  ,TEMPEX,TEMPEY,TEMPEZ,TEMPIX
     5  ,TEMPIY,TEMPIZ
     6  ,CUEX,CUEY,CUEZ,CUIX,CUIY,CUIZ
     7  ,TPLCE1,TPLCE2,TPALCE,TPLCI1,TPLCI2,TPALCI
      NC=NCLAST
      GO TO 1001
 1000 CONTINUE

      OPEN(15,FORM='UNFORMATTED')
      OPEN(16,FORM='UNFORMATTED')
      OPEN(17,FORM='UNFORMATTED')
C
C INITIALIZATION OF PARTICLE POSITIONS AND VELOCITIES
C
C  INITIALIZE IONS
      SIX=0.
      SIY=0.
      SIZ=0.
C+
      ISEED = 0
      CALL RANN1(0.,1.,ISEED,RANORM,NOX*NOY*3,IRC)
C
      DO 955 M=1,NOX
      ADPX=(M-1)*DDELX+BNCXL
      DO 955 N=1,NOY
      ADPY=(N-1)*DDELY+ANCYL
      I=N+(M-1)*NOY
      J=I+NO
      X(J)=ADPX
      Y(J)=ADPY
      PX(J)=VTXI*RANORM(3*I-2)
      SIX=SIX+PX(J)
      PY(J)=VTYI*RANORM(3*I-1)
      SIY=SIY+PY(J)
      PZ(J)=VTZI*RANORM(3*I)
 955  SIZ=SIZ+PZ(J)
C  BACKSTAGGER IONS
      SIX=SIX/NO
      SIY=SIY/NO
      SIZ=SIZ/NO
      DO 888 J=NO1,NO2
      PX(J)=PX(J)-SIX
      PY(J)=PY(J)-SIY
      PZ(J)=PZ(J)-SIZ
      PT=PX(J)*PX(J)+PY(J)*PY(J)+PZ(J)*PZ(J)
      PT=SQRT(1.+PT*CIN)
      PSA(J)=PT
      VX=PX(J)/PSA(J)
      VY=PY(J)/PSA(J)
      X(J)=X(J)-VC*VX
      Y(J)=Y(J)-VC*VY
C     XBF(J)=X(J)
C++
C++  REFLECT PARTICLES AT WALL
C++
      IF (X(J).LE.BNCXL) THEN
      X(J)=BNCXL2-X(J)
      IF (RESIST) THEN
      PX(J)=-PX(J)
      PY(J)=-PY(J)
      ELSE
      PX(J)=-PX(J)
      PZ(J)=-PZ(J)
      ENDIF
      ELSE
      IF (X(J).GE.BNCXR) THEN
      X(J)=BNCXR2-X(J)
      IF (RESIST) THEN
      PX(J)=-PX(J)
      PY(J)=-PY(J)
      ELSE
      PX(J)=-PX(J)
      PZ(J)=-PZ(J)
      ENDIF
      ENDIF
      ENDIF
C++
C++ END TEST
C++
  888 CONTINUE
C  INITIALIZE ELECTRONS
      SVX=0.
      SVY=0.
      SVZ=0.
      VDF=0.
C+
      CALL RANN1(0.,1.,ISEED,RANORM,NOX*NOY*3,IRC)
C
      DO 555 M=1,NOX
      DO 555 N=1,NOY
      I=N+(M-1)*NOY
      J=I+NO
      X(I)=X(J)
      Y(I)=Y(J)
      VXR=VTX*RANORM(3*I-2)
      PX(I)=VXR+VDX
      SVX=SVX+VXR
      VYR=VTY*RANORM(3*I-1)
      PY(I)=VYR+VDY
      SVY=SVY+VYR
      VZR=VTZ*RANORM(3*I)
      PZ(I)=VZR+VDZ
      SVZ=SVZ+VZR
 555  CONTINUE
C  BACKSTAGGER ELECTRONS
      SVX=SVX/NO
      SVY=SVY/NO
      SVZ=SVZ/NO
      DO 506 I=1,NO
      PX(I)=PX(I)-SVX
      PY(I)=PY(I)-SVY
      PZ(I)=PZ(I)-SVZ
      PT=PX(I)*PX(I)+PY(I)*PY(I)+PZ(I)*PZ(I)
      PT=SQRT(1.+PT*CIN)
      PSA(I)=PT
      VX=PX(I)/PSA(I)
      VY=PY(I)/PSA(I)
      VDF=VDF+VX
      X(I)=X(I)-VC*VX
      Y(I)=Y(I)-VC*VY
C     XBF(I)=X(I)
C++
C++  REFLECT PARTICLES AT WALL
C++
      IF (X(J).LE.BNCXL) THEN
      X(J)=BNCXL2-X(J)
      IF (RESIST) THEN
      PX(J)=-PX(J)
      PY(J)=-PY(J)
      ELSE
      PX(J)=-PX(J)
      PZ(J)=-PZ(J)
      ENDIF
      ELSE
      IF (X(J).GE.BNCXR) THEN
      X(J)=BNCXR2-X(J)
      IF (RESIST) THEN
      PX(J)=-PX(J)
      PY(J)=-PY(J)
      ELSE
      PX(J)=-PX(J)
      PZ(J)=-PZ(J)
      ENDIF
      ENDIF
      ENDIF
C++
C++  END TEST
C++
 506  CONTINUE
      VDF=VDF/NO
C++
C++  STORE PHASE SPACE INFORMATION ON DISK AT T=0
C++
      WRITE(25) (X(I),I=1,PART)
      WRITE(26) (Y(I),I=1,PART)
      WRITE(22) (PX(I),I=1,PART)
      WRITE(23) (PY(I),I=1,PART)
      WRITE(24) (PZ(I),I=1,PART)
C++
C++  TO PLOT ONLY FEW POINTS
C++
      J=0
      DO 9 I=1,N4O
      J=J+MUS
      XG(I)=X(J)
      YG(I)=Y(J)
      PXG(I)=PX(J)
      PYG(I)=PY(J)
      PZG(I)=PZ(J)
  9   CONTINUE
C++
C++  DRAWING MASKING FUNCTION
C++
      CALL GRAF1(FM,1,NCX,'MASKF',0)
C++
C++ PLOT DISTRIBUTION FUNCTIONS IN OBLIQUE RUN
C++
      CALL OBLIC(COST,SINT)
      CALL RMINMX(PER1,1,N8O,1,ER1MIE,ER1MAE)
      CALL RMINMX(PER2,1,N8O,1,ER2MIE,ER2MAE)
      CALL RMINMX(PAR,1,N8O,1,ARMINE,ARMAXE)
      CALL RMINMX(PER1,N8O1,N4O,1,ER1MII,ER1MAI)
      CALL RMINMX(PER2,N8O1,N4O,1,ER2MII,ER2MAI)
      CALL RMINMX(PAR,N8O1,N4O,1,ARMINI,ARMAXI)
      CALL DIST(PER1,N8O,ER1MAE,ER1MIE,'P1 EL',LABL2)
      CALL DIST(PER2,N8O,ER2MAE,ER2MIE,'P2 EL',LABL2)
      CALL DIST(PAR,N8O,ARMAXE,ARMINE,'PAR EL',LABL2)
      CALL DIST(PER1(N8O1),N8O,ER1MAI,ER1MII,'P1 IO',LABL2)
      CALL DIST(PER2(N8O1),N8O,ER2MAI,ER2MII,'P2 IO',LABL2)
      CALL DIST(PAR(N8O1),N8O,ARMAXI,ARMINI,'PAR IO',LABL2)
 1001 CONTINUE
C
C  INITIALIZE TIMES FOR ENERGIES CALCULATION
C
      NDT=NSS
C++
C+++    POUR MEMORISER LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART
       NKPART=0
C**********************************
C    MAIN TIME LOOP
C**********************************
  100 CONTINUE
      NT=NT+1
      T=NT*DT
      NS=(NT/NDT)*NDT
      NL=(NT/NLL)*NLL
C++
C+++    POUR MEMORISER LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART

	    NTPART=(NT/NTT)*NTT
      NTIME=(NT/NTIM)*NTIM
C
C
C CALCULATION OF CURRENT FROM PARTICLE VELOCITY
C  AND CHARGE DENSITY FROM PARTICLE POSITIONS
C  BY USING SUD METHOD
C
      DO 1 J=1,NCY
      DO 1 I=1,NCX
      G1X(I,J)=0.
      CUX(I,J)=0.
      CUY(I,J)=0.
      CUZ(I,J)=0.
    1 CONTINUE
C  ION CHARGE AND CURRENTS
      IF(MOVION) THEN
      CALL DENSIT(I5,I6,VC,CUX,CUY,CUZ,G1X)
      DO 402 J=1,NCY
      DO 402 I=1,NCX
      CUX(I,J)=-CUX(I,J)
      CUY(I,J)=-CUY(I,J)
      CUZ(I,J)=-CUZ(I,J)
C  STORE THE ION CURRENT
      CUXION(I,J)=CUX(I,J)
      CUYION(I,J)=CUY(I,J)
      CUZION(I,J)=CUZ(I,J)
      G1XION(I,J)=G1X(I,J)
 402  CONTINUE
      ENDIF
C  ELECTRON CHARGE AND CURRENT
      DO 407 J=1,NCY
      DO 407 I=1,NCX
      G1X(I,J)=0.
 407  CONTINUE
      CALL DENSIT(I3,I4,VC,CUX,CUY,CUZ,G1X)
C++
C++  ADD EXTERNAL CURRENT BETWEEN NCUO AND NCUE
C++  THIS IS THE TUBE VERSIOn
      IF (.NOT. RESIST) THEN
      DO 9090 KK=LPLX-2*RAYON-NCUE,LPLX-NCUE
      AA=(KK-POSIT)*(KK-POSIT)
      CIRTEM=SQRT(RAYON*RAYON-AA)+0.5
      K1=CIRTEM
      K2=-1.*K1+FIELD2/2
      K1=K1+FIELD2/2
      IF (NT.GT.NPE2) THEN
      DO 993 J=K2,K1
      CUZ(KK,J)=CUZ(KK,J)+CURREN
      CUZAP(KK,J)=CURREN
  993 CONTINUE
      ELSE
      DO 994 J=K2,K1
      CUZ(KK,J)=CUREN2-CUREN2*SIN((TWOPI*NT/NPE)+PI/2.)
      CUZAP(KK,J)=CUREN2-CUREN2*SIN((TWOPI*NT/NPE)+PI/2.)
  994 CONTINUE
      ENDIF
 9090 CONTINUE
      ELSE
      IF (NT.GT.NPE2) THEN
      DO 992 I=NCUO,NCUE
      DO 992 J=1,NCY
      CUZ(I,J)=CUZ(I,J)+CURREN
      CUZAP(I,J)=CURREN
  992 CONTINUE
      ELSE
      DO 991 I=NCUO,NCUE
      DO 991 J=1,NCY
      CUZ(I,J)=CUREN2-CUREN2*SIN((TWOPI*NT/NPE)+PI/2.)
      CUZAP(I,J)=CUREN2-CUREN2*SIN((TWOPI*NT/NPE)+PI/2.)
  991 CONTINUE
      ENDIF
      ENDIF
      IF(NT.NE.NL) GO TO 625

      DO 603 I=1,NCX
      CUXX(I)=0.
      G1XX(I)=0.
      CUXIOX(I)=0.
      G1XIOX(I)=0.
  603 CONTINUE
*VOCL LOOP,VI(I)
      DO 604 J=1,NCY
      DO 604 I=LPLX,NCX
      G1XX(I)=G1XX(I)+G1X(I,J)
      G1XIOX(I)=G1XIOX(I)+G1XION(I,J)
      CUXX(I)=CUXX(I)+(CUX(I,J)-CUXION(I,J))
      CUXIOX(I)=CUXIOX(I)+CUXION(I,J)
  604 CONTINUE
      DO 605 I=LPLX,NCX
      G1XX(I)=G1XX(I)/NCY
      CUXIOX(I)=CUXIOX(I)/NCY
      CUXX(I)=CUXX(I)/NCY
      G1XIOX(I)=G1XIOX(I)/NCY
 605  CONTINUE
      CALL GRAF1(G1XX,LPLX,NCX,'G1X-X',0)
      CALL GRAF1(G1XIOX,LPLX,NCX,'G1XIO-X',0)
      CALL GRAF1(CUXX,LPLX,NCX,'CUX-X',0)
      CALL GRAF1(CUXIOX,LPLX,NCX,'CUXIO-X',0)

      DO 610 I=1,NCX
      CUZX(I)=0.
      CUZIOX(I)=0.
      CUYX(I)=0.
      CUYIOX(I)=0.
  610 CONTINUE
*VOCL LOOP,VI(I)
      DO 611 J=1,NCY
      DO 611 I=LPLX,NCX
      CUZX(I)=CUZX(I)+(CUZ(I,J)-CUZION(I,J))
      CUZIOX(I)=CUZIOX(I)+CUZION(I,J)
      CUYX(I)=CUYX(I)+(CUY(I,J)-CUYION(I,J))
      CUYIOX(I)=CUYIOX(I)+CUYION(I,J)
  611 CONTINUE
      DO 612 I=LPLX,NCX
      CUZX(I)=CUZX(I)/NCY
      CUZIOX(I)=CUZIOX(I)/NCY
      CUYX(I)=CUYX(I)/NCY
      CUYIOX(I)=CUYIOX(I)/NCY
  612 CONTINUE
      CALL GRAF1(CUYX,LPLX,NCX,'CUY-X',0)
      CALL GRAF1(CUYIOX,LPLX,NCX,'CUYIO-X',0)
      CALL GRAF1(CUZX,LPLX,NCX,'CUZ-X',0)
      CALL GRAF1(CUZIOX,LPLX,NCX,'CUZIO-X',0)

 625  CONTINUE
C
C  CALCULATE TOTAL CHARGE DENSITY
C
*VOCL LOOP,VI(I)
      DO 406 J=1,NCY
      DO 406 I=1,NCX
      G1X(I,J)=G1X(I,J)-G1XION(I,J)
 406  CONTINUE
C TRANSFORM CHARGE AND TOTAL CURRENT COMPONENTS TO K SPACE
C
      IFSET=-1
      CALL NVCFT2(G1X,N2,NCX,NCY,IFSET,TABLE,WORK)
      CALL NVCFT2(CUX,N2,NCX,NCY,IFSET,TABLE,WORK)
      CALL NVCFT2(CUY,N2,NCX,NCY,IFSET,TABLE,WORK)
      CALL NVCFT2(CUZ,N2,NCX,NCY,IFSET,TABLE,WORK)
C
C  CALCULATION OF LONGITUDINAL ELECTRIC FIELDS
C
      DO 451 M=1,NCXX
      LE=M+M
      LO=LE-1
      MM=M
      IF(M.EQ.1) MM=NCXH1
      DO 452 N=1,NCYH1
      FXNEW=FEX(M,N)
      FYNEW=FEY(M,N)
C CALCULATE THE POTENTIAL
      POT(LO,N)=PK(M,N)*G1X(LO,N)
      POT(LE,N)=PK(M,N)*G1X(LE,N)
      ELX(LO,N)=FXNEW*G1X(LE,N)
      ELX(LE,N)=-FXNEW*G1X(LO,N)
      ELY(LO,N)=FYNEW*G1X(LE,N)
  452 ELY(LE,N)=-FYNEW*G1X(LO,N)
      DO 451 N=NCYH2,NCY
      K=NCYP2-N
      FXNEW=FEX(MM,K)
      FYNEW=-FEY(MM,K)
C NEXT LINE IS FOR POTENTIAL CALCULATION.
      POT(LO,N)=PK(MM,K)*G1X(LO,N)
      POT(LE,N)=PK(MM,K)*G1X(LE,N)
      ELX(LO,N)=FXNEW*G1X(LE,N)
      ELX(LE,N)=-FXNEW*G1X(LO,N)
      ELY(LO,N)=FYNEW*G1X(LE,N)
  451 ELY(LE,N)=-FYNEW*G1X(LO,N)
C
C  CALCULATION OF THE SELF-CONSISTENT MAGNETIC FIELD
C  CALCULATION OF THE SELF-CONSISTENT ELECTRIC FIELD
C
      DO 453 M=1,NCXX
      LE=M+M
      LO=LE-1
      MM=M
      IF(M.EQ.1) MM=NCXH1
      DO 454 N=1,NCYH1
      BTX2=BTX(LE,N)
      BTX1=BTX(LO,N)
      BTY2=BTY(LE,N)
      BTY1=BTY(LO,N)
      BTZ2=BTZ(LE,N)
      BTZ1=BTZ(LO,N)
      FF1=REX(M,N)
      FF2=REY(M,N)
      FF3=CKX(M,N)
      FF4=CKY(M,N)
      FF5=CKZ(M,N)
      FF6=CKC(M,N)
      ETX(LO,N)=ETX(LO,N)-FF2*BTZ(LE,N)-FF3*CUX(LO,N)
     1  +FF6*CUY(LO,N)
      ETX(LE,N)=ETX(LE,N)+FF2*BTZ(LO,N)-FF3*CUX(LE,N)
     1  +FF6*CUY(LE,N)
      ETY(LO,N)=ETY(LO,N)+FF1*BTZ(LE,N)-FF4*CUY(LO,N)
     1  +FF6*CUX(LO,N)
      ETY(LE,N)=ETY(LE,N)-FF1*BTZ(LO,N)-FF4*CUY(LE,N)
     1  +FF6*CUX(LE,N)
      ETZ(LO,N)=ETZ(LO,N)-FF1*BTY(LE,N)+FF2*BTX(LE,N)
     1  -FF5*CUZ(LO,N)
      ETZ(LE,N)=ETZ(LE,N)+FF1*BTY(LO,N)-FF2*BTX(LO,N)
     1  -FF5*CUZ(LE,N)
      BTX(LO,N)=BTX(LO,N)+FF2*ETZ(LE,N)
      BTX(LE,N)=BTX(LE,N)-FF2*ETZ(LO,N)
      BTY(LO,N)=BTY(LO,N)-FF1*ETZ(LE,N)
      BTY(LE,N)=BTY(LE,N)+FF1*ETZ(LO,N)
      BTZ(LO,N)=BTZ(LO,N)+FF1*ETY(LE,N)-FF2*ETX(LE,N)
      BTZ(LE,N)=BTZ(LE,N)-FF1*ETY(LO,N)+FF2*ETX(LO,N)
      BRX(LE,N)=0.5*(BTX(LE,N)+BTX2)
      BRX(LO,N)=0.5*(BTX(LO,N)+BTX1)
      BRY(LE,N)=0.5*(BTY(LE,N)+BTY2)
      BRY(LO,N)=0.5*(BTY(LO,N)+BTY1)
      BRZ(LE,N)=0.5*(BTZ(LE,N)+BTZ2)
      BRZ(LO,N)=0.5*(BTZ(LO,N)+BTZ1)
  454 CONTINUE
      DO 453 N=NCYH2,NCY
      K=NCYP2-N
      BTX2=BTX(LE,N)
      BTX1=BTX(LO,N)
      BTY2=BTY(LE,N)
      BTY1=BTY(LO,N)
      BTZ2=BTZ(LE,N)
      BTZ1=BTZ(LO,N)
      FF1=REX(MM,K)
      FF2=-REY(MM,K)
      FF3=CKX(MM,K)
      FF4=CKY(MM,K)
      FF5=CKZ(MM,K)
      FF6=-CKC(MM,K)
      ETX(LO,N)=ETX(LO,N)-FF2*BTZ(LE,N)-FF3*CUX(LO,N)
     1  +FF6*CUY(LO,N)
      ETX(LE,N)=ETX(LE,N)+FF2*BTZ(LO,N)-FF3*CUX(LE,N)
     1  +FF6*CUY(LE,N)
      ETY(LO,N)=ETY(LO,N)+FF1*BTZ(LE,N)-FF4*CUY(LO,N)
     1  +FF6*CUX(LO,N)
      ETY(LE,N)=ETY(LE,N)-FF1*BTZ(LO,N)-FF4*CUY(LE,N)
     1  +FF6*CUX(LE,N)
      ETZ(LO,N)=ETZ(LO,N)-FF1*BTY(LE,N)+FF2*BTX(LE,N)
     1  -FF5*CUZ(LO,N)
      ETZ(LE,N)=ETZ(LE,N)+FF1*BTY(LO,N)-FF2*BTX(LO,N)
     1  -FF5*CUZ(LE,N)
      BTX(LO,N)=BTX(LO,N)+FF2*ETZ(LE,N)
      BTX(LE,N)=BTX(LE,N)-FF2*ETZ(LO,N)
      BTY(LO,N)=BTY(LO,N)-FF1*ETZ(LE,N)
      BTY(LE,N)=BTY(LE,N)+FF1*ETZ(LO,N)
      BTZ(LO,N)=BTZ(LO,N)+FF1*ETY(LE,N)-FF2*ETX(LE,N)
      BTZ(LE,N)=BTZ(LE,N)-FF1*ETY(LO,N)+FF2*ETX(LO,N)
      BRX(LE,N)=0.5*(BTX(LE,N)+BTX2)
      BRX(LO,N)=0.5*(BTX(LO,N)+BTX1)
      BRY(LE,N)=0.5*(BTY(LE,N)+BTY2)
      BRY(LO,N)=0.5*(BTY(LO,N)+BTY1)
      BRZ(LE,N)=0.5*(BTZ(LE,N)+BTZ2)
      BRZ(LO,N)=0.5*(BTZ(LO,N)+BTZ1)
  453 CONTINUE
C++
C++ TRANSFORM ELECTRIC FIELD TO REAL SPACE
C++ FOR MASKING
C++
      CALL NVCFT2(ETX,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ETY,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ETZ,N2,NCX,NCY,1,TABLE,WORK)
*VOCL LOOP,VI(I)
      DO 61 J=1,NCY
      DO 61 I=1,NCX
      ETX(I,J)=FM(I)*ETX(I,J)
      ETY(I,J)=FM(I)*ETY(I,J)
      ETZ(I,J)=FM(I)*ETZ(I,J)
  61  CONTINUE
C++
C++  TRANSFORM TRANSVERSE ELECTRIC FIELD BACK TO K-SPACE
C++
      CALL NVCFT2(ETX,N2,NCX,NCY,-1,TABLE,WORK)
      CALL NVCFT2(ETY,N2,NCX,NCY,-1,TABLE,WORK)
      CALL NVCFT2(ETZ,N2,NCX,NCY,-1,TABLE,WORK)
C
C  CALCULATE AND STORE ALL ENERGIES
C
      IF(NS.NE.NT) GO TO 32
C  FIELDS
      CALL FENERG(ELX,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,EXS)
      CALL FENERG(ELY,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,EYS)
      CALL FENERG(BTX,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,BXS)
      CALL FENERG(BTY,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,BYS)
      CALL FENERG(BTZ,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,BZS)
      EB=BXS+BYS+BZS
      EB=EB+EB
      EL=EXS+EYS
      EL=EL+EL
      CALL FENERG(ETX,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,WXS)
      CALL FENERG(ETY,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,WYS)
      CALL FENERG(ETZ,AK,LVGC,NCX,NCXX,NCXH1,NCYH1,NCYH2,NCY,WZS)
      ET=WXS+WYS+WZS
      ET=ET+ET
C++
C++  ENERGY FIELDS
C++
      EXS=EXS+EXS
      EYS=EYS+EYS
      BXS=BXS+BXS
      BYS=BYS+BYS
      BZS=BZS+BZS
      WXS=WXS+WXS
      WYS=WYS+WYS
      WZS=WZS+WZS
C  KINETIC, TEMPERATURE AND CURRENT
      EKE=0.
      VDY=0.
      VDZ=0.
      VDF=0.
      DO 23 I=1,NO
      EKE=EKE+PSA(I)
      VX=PX(I)/PSA(I)
      VY=PY(I)/PSA(I)
      VZ=PZ(I)/PSA(I)
      VDF=VDF+VX
      VDY=VDY+VY
      VDZ=VDZ+VZ
 23   CONTINUE
      EKE=2.*(EKE-NO)*RTE
      IF(.NOT.MOVION) GO TO 706
      EKI=0.
      DO 707 J=NO1,NO2
      EKI=EKI+PSA(J)
  707 CONTINUE
      EKI=2.*(EKI-NO)*RTE
  706 CONTINUE
      VDF=VDF/NO
      VDY=VDY/NO
      VDZ=VDZ/NO
      EKI=EKI/RMASS
      EKT=EKE+EKI
      TE=EL+ET+EB+EKT
      CALL TEMP(TPEX,TPEY,TPEZ,TPIX,TPIY,TPIZ)
      CALL COUR(CEX,CEY,CEZ,CIX,CIY,CIZ)
C++
C++ CALCULATE NEW PERP. AND PARA. TEMPERATURE
C++
      IF (RESIST) THEN
      CALL TEMPR(TPE1,TPE2,TPAE,TPI1,TPI2,TPAI,
     &  CELEL,CELIO,NT,LPLX)
      ELSE
      CALL TEMPNR(TPE1,TPE2,TPAE,TPI1,TPI2,TPAI,
     &  CELEL,CELIO,NT,LPLX)
      ENDIF
C  STORE
      NC=NC+1
      TPLCE1(NC)=TPE1
      TPLCE2(NC)=TPE2
      TPALCE(NC)=TPAE
      TPLCI1(NC)=TPI1
      TPLCI2(NC)=TPI2
      TPALCI(NC)=TPAI
      EFLX(NC)=EXS
      EFLY(NC)=EYS
      EFTX(NC)=WXS
      EFTY(NC)=WYS
      EFTZ(NC)=WZS
      BFTX(NC)=BXS
      BFTY(NC)=BYS
      BFTZ(NC)=BZS
      EFL(NC)=EL
      EFT(NC)=ET
      BFT(NC)=EB
      EFIEL(NC)=EL+ET+EB
      EKINE(NC)=EKE
      EKINI(NC)=EKI
      EKIN(NC)=EKT
      ETOT(NC)=TE
      TEMPEX(NC)=TPEX
      TEMPEY(NC)=TPEY
      TEMPEZ(NC)=TPEZ
      TEMPIX(NC)=TPIX
      TEMPIY(NC)=TPIY
      TEMPIZ(NC)=TPIZ
      CUEX(NC)=CEX
      CUEY(NC)=CEY
      CUEZ(NC)=CEZ
      CUIX(NC)=CIX
      CUIY(NC)=CIY
      CUIZ(NC)=CIZ
      IF(CIX.LE.1.0E-5) GO TO 32
      RATIO(NC)=CEX/CIX
*
 32   CONTINUE

*VOCL LOOP,VI(I)
      DO 962 J=1,NCY
      DO 962 I=1,NCX
      ERX(I,J)=ETX(I,J)+ELX(I,J)
      ERY(I,J)=ETY(I,J)+ELY(I,J)
      ERZ(I,J)=ETZ(I,J)
      ERX1(I,J)=ETX(I,J)
      ERY1(I,J)=ETY(I,J)
  962 CONTINUE
C++
C++  TRANSFORM POTENTIEL FROM K-SPACE TO REAL SPACE
C++
      CALL NVCFT2(POT,N2,NCX,NCY,1,TABLE,WORK)
C++
C++  CALCULATE FREE SPACE CORRECTION FOR LONGITUDINAL
C++  FIELDS :THIS CORRECTION ADDED TO REAL SPACE
C++
      AELX=0.D0
      DO 233 N=1,NCY
      AELX=AELX+(POT(LXLPOT,N)-POT(LXRPOT,N))
  233 CONTINUE
      AELX=AELX/NCY
      AELX=AELX/(LXRPOT-LXLPOT)
C
C  TRANSFORM FIELDS FROM K-SPACE TO REAL SPACE
C
      CALL NVCFT2(ELX,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ELY,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(BRX,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(BRZ,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(BRY,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ERX,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ERY,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ERZ,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ERX1,N2,NCX,NCY,1,TABLE,WORK)
      CALL NVCFT2(ERY1,N2,NCX,NCY,1,TABLE,WORK)
C++
C++  CORRECT ELX FOR FREE SPACE
C++
      DO 150 J=1,NCY
      DO 150 I=1,NCX
      ELX(I,J)=ELX(I,J)+AELX
      ERX(I,J)=ERX(I,J)+AELX
  150 CONTINUE
C++
C++  STORE FIELDS COMPONENTS VERSUS X AND Y ON DISK (NO AVERAGING)
C++
      IF (NT.NE.NL) GO TO 111
      WRITE(30) ((ELX(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(31) ((ELY(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(32) ((ERX1(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(33) ((ERY1(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(34) ((ERZ(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(35) ((BRX(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(36) ((BRY(I,J),J=1,NCY),I=LPLX,NCX)
      WRITE(37) ((BRZ(I,J),J=1,NCY),I=LPLX,NCX)
      AT=FLOAT(NT)
      WRITE(39) AT

  111 CONTINUE
C++
C++ AVERAGE FIELDS VERSUS Y PLOTTED VERSUS X
C++  ELX,ELY,ETX,ETY,ETZ,BTX,BTY,BTZ
C++
      IF (NT.NE.NL) GOTO 412
      DO 515 I=1,NCX
      ELXYX(I)=0.
      ELYYX(I)=0.
      ETXYX(I)=0.
      ETYYX(I)=0.
      ETZYX(I)=0.
      BTXYX(I)=0.
      BTYYX(I)=0.
      BTZYX(I)=0.
      POTYX(I)=0.
  515 CONTINUE
*VOCL LOOP,VI(I)
      DO 516 J=1,NCY
      DO 516 I=1,NCX
      ELXYX(I)=ELXYX(I)+ELX(I,J)
      ELYYX(I)=ELYYX(I)+ELY(I,J)
      ETXYX(I)=ETXYX(I)+ERX1(I,J)
      ETYYX(I)=ETYYX(I)+ERY1(I,J)
      ETZYX(I)=ETZYX(I)+ERZ(I,J)
      BTXYX(I)=BTXYX(I)+BRX(I,J)
      BTYYX(I)=BTYYX(I)+BRY(I,J)
      BTZYX(I)=BTZYX(I)+BRZ(I,J)
      POTYX(I)=POTYX(I)+POT(I,J)
  516 CONTINUE
      DO 517 I=1,NCX
      ELXYX(I)=ELXYX(I)/NCY
      ELYYX(I)=ELYYX(I)/NCY
      ETXYX(I)=ETXYX(I)/NCY
      ETYYX(I)=ETYYX(I)/NCY
      ETZYX(I)=ETZYX(I)/NCY
      BTXYX(I)=BTXYX(I)/NCY
      BTYYX(I)=BTYYX(I)/NCY
      BTZYX(I)=BTZYX(I)/NCY
      POTYX(I)=POTYX(I)/NCY
  517 CONTINUE
      CALL GRAF1(ELXYX,LPLX,NCX,'ELX Y-X',0)
      CALL GRAF1(ELYYX,LPLX,NCX,'ELY Y-X',0)
      CALL GRAF1(ETXYX,LPLX,NCX,'ETX Y-X',0)
      CALL GRAF1(ETYYX,LPLX,NCX,'ETY Y-X',0)
      CALL GRAF1(ETZYX,LPLX,NCX,'ETZ Y-X',0)
      CALL GRAF1(BTXYX,LPLX,NCX,'BTX Y-X',0)
      CALL GRAF1(BTYYX,LPLX,NCX,'BTY Y-X',0)
      IF (.NOT.RESIST) CALL GRAF1(BTYYX,1,NCX,'BTY Y-X',0)
      CALL GRAF1(BTZYX,LPLX,NCX,'BTZ Y-X',0)
      IF (RESIST) CALL GRAF1(BTZYX,1,NCX,'BTZ Y-X',0)
      CALL GRAF1(POTYX,LPLX,NCX,'POT Y-X',0)

  412 CONTINUE

C DEFINE TOTAL B FIELD
C
      DO 961 J=1,NCY
      DO 961 I=1,NCX
      BRX(I,J)=BRX(I,J)+BXO
      BRY(I,J)=BRY(I,J)+BYO
      BRZ(I,J)=BRZ(I,J)+BZO
  961 CONTINUE

C
C  ERXAV, ERYAV INTEGRATED ALONG X AND PLOTTED VERSUS Y

C
C  STORE TEST PARTICLE ORBITS ON DISK
C
      NCT=NCT+1
      IF(NT.NE.NL) GO TO 7000

C++
C++  STORE PHASE SPACE INFORMATION ON DISK AT T=NL
C++
C++  TO PLOT ONLY FEW  POINTS
C++
      J=0
      DO 7 I=1,N4O
      J=J+MUS
      XG(I)=X(J)
      YG(I)=Y(J)
      PXG(I)=PX(J)
      PYG(I)=PY(J)
      PZG(I)=PZ(J)
  7   CONTINUE
      WRITE(25) (X(I),I=1,PART)
      WRITE(26) (Y(I),I=1,PART)
      WRITE(22) (PX(I),I=1,PART)
      WRITE(23) (PY(I),I=1,PART)
      WRITE(24) (PZ(I),I=1,PART)
*
      IF (MEMOP) THEN
      WRITE(50) XELSP,YELSP,XIOSP,YIOSP,NTSP
      ENDIF
*
C
C  PLOT DISTRIBUTIONS
C  PLOT B FIELD LINES
C
C  PARTICLES
C++
C++ TO PLOT IN OBLIQUE PROPAGATION
C++
      CALL OBLIC(COST,SINT)
      CALL RMINMX(PER1,1,N8O,1,ER1MIE,ER1MAE)
      CALL RMINMX(PER2,1,N8O,1,ER2MIE,ER2MAE)
      CALL RMINMX(PAR,1,N8O,1,ARMINE,ARMAXE)
      CALL RMINMX(PER1,N8O1,N4O,1,ER1MII,ER1MAI)
      CALL RMINMX(PER2,N8O1,N4O,1,ER2MII,ER2MAI)
      CALL RMINMX(PAR,N8O1,N4O,1,ARMINI,ARMAXI)
      CALL DIST(PER1,N8O,ER1MAE,ER1MIE,'P1 EL',LABL2)
      CALL DIST(PER2,N8O,ER2MAE,ER2MIE,'P2 EL',LABL2)
      CALL DIST(PAR,N8O,ARMAXE,ARMINE,'PAR EL',LABL2)
      CALL DIST(PER1(N8O1),N8O,ER1MAI,ER1MII,'P1 IO',LABL2)
      CALL DIST(PER2(N8O1),N8O,ER2MAI,ER2MII,'P2 IO',LABL2)
      CALL DIST(PAR(N8O1),N8O,ARMAXI,ARMINI,'PAR IO',LABL2)

      CALL PHASE(ER1MIE,ER1MII,ER1MAE,ER1MAI,ER2MIE
     1         ,ER2MII,ER2MAE,ER2MAI,NO,LABL2)
      CALL PHASEX(XMIN,XMAX,ER1MIE,ER1MII,ER1MAE,ER1MAI,ER2MIE
     1         ,ER2MII,ER2MAE,ER2MAI
     2         ,ARMINE,ARMINI,ARMAXE,ARMAXI,NO,LABL2)
      CALL PHASEY(YMIN,YMAX,ER1MIE,ER1MII,ER1MAE,ER1MAI,ER2MIE
     1         ,ER2MII,ER2MAE,ER2MAI
     2         ,ARMINE,ARMINI,ARMAXE,ARMAXI,NO,LABL2)

 7000 CONTINUE
C
C  PUSH ELECTRONS AND IONS, IF MOVION=TRUE
C
      SWITCH=.FALSE.
      BETA=BETE
      ALFA=ALFE

C++ I3,I4,I5 AND I4 ARE PARAMETER FOR THE VECTORISATION ON VP200
C++
      I8=I3
      I9=I4
   15 CONTINUE
      DO 40 M=I8,I9
      DO 40 J=1,2048
      I=J+(M-1)*2048
      LX=X(I)+1.5
      DX=(X(I)-LX+1.)*0.5
      LXL=IL(LX)
      LXR=IR(LX)
      LY=Y(I)+1.5
      DY=(Y(I)-LY+1)*0.5
      LYL=JL(LY)
      LYR=JR(LY)
      EXT=ERX(LX,LY)+DX*(ERX(LXR,LY)-ERX(LXL,LY))
     1+DY*(ERX(LX,LYR)-ERX(LX,LYL))
      EYT=ERY(LX,LY)+DX*(ERY(LXR,LY)-ERY(LXL,LY))
     1+DY*(ERY(LX,LYR)-ERY(LX,LYL))
      EZT=ERZ(LX,LY)+DX*(ERZ(LXR,LY)-ERZ(LXL,LY))
     1+DY*(ERZ(LX,LYR)-ERZ(LX,LYL))
      BXT=BRX(LX,LY)+DX*(BRX(LXR,LY)-BRX(LXL,LY))
     1+DY*(BRX(LX,LYR)-BRX(LX,LYL))
      BYT=BRY(LX,LY)+DX*(BRY(LXR,LY)-BRY(LXL,LY))
     1+DY*(BRY(LX,LYR)-BRY(LX,LYL))
      BZT=BRZ(LX,LY)+DX*(BRZ(LXR,LY)-BRZ(LXL,LY))
     1+DY*(BRZ(LX,LYR)-BRZ(LX,LYL))
      EXT=BETA*EXT
      EYT=BETA*EYT
      EZT=BETA*EZT
      PX(I)=PX(I)+EXT
      PY(I)=PY(I)+EYT
      PZ(I)=PZ(I)+EZT
      PT=PX(I)*PX(I)+PY(I)*PY(I)+PZ(I)*PZ(I)
      PRT=SQRT(1.+PT*CIN)
      PRTI=1./PRT
      BXT=BXT*PRTI
      BYT=BYT*PRTI
      BZT=BZT*PRTI
      B1=ALFA*BXT
      B2=B1*B1
      B3= ALFA*BYT
      B4=B3*B3
      B5= ALFA*BZT
      B6=B5*B5
      B7=B1*B3
      B8=B1*B5
      B9=B3*B5
      C1=1.+B2
      C2=B5+B7
      C3=B8-B3
      C4=1.+B2-B4-B6
      C5=2.*C2
      C6=2.*C3
      C7=B7-B5
      C8=1.+B4
      C9=B9+B1
      C10=2.*C7
      C11=1.+B4-B2-B6
      C12=2.*C9
      C13=B3+B8
      C14=B9-B1
      C15=1.+B6
      C16=2.*C13
      C17=2.*C14
      C18=1.+B6-B2-B4
      SUM=1./(1.+B2+B4+B6)
      TVX=SUM*(C4*PX(I)+C5*PY(I)+C6*PZ(I))
      TVY=SUM*(C10*PX(I)+C11*PY(I)+C12*PZ(I))
      TVZ=SUM*(C16*PX(I)+C17*PY(I)+C18*PZ(I))
      PX(I)=EXT+TVX
      PY(I)=EYT+TVY
      PZ(I)=EZT+TVZ
      PT=PX(I)*PX(I)+PY(I)*PY(I)+PZ(I)*PZ(I)
      PT=SQRT(1.+PT*CIN)
      PSA(I)=PT
   40 CONTINUE
      IF(MOVION) THEN
      IF(SWITCH) GO TO 14
      ALFA=ALFI
      BETA=BETI
C--   I2=NO2
C--   I1=NO1
C++
C++  PARAMETER FOR THE VECTORISATION ON VP200
      I8=I5
      I9=I6
      SWITCH=.TRUE.
      GO TO 15
   14 CONTINUE
      ENDIF
C++
C+++    ICI ON MEMORISE LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART
        IF (NT .EQ. NTPART) THEN
        NKPART=NKPART+1
      
         IF (NUMPAE .NE. 1) CALL TIMPA(NOMPAE,NUMPAE,NEMDI,NKPART,ELX,
     &ELY,ERX1,ERY1,ERZ,BRX,BRY,BRZ,TRAJEL)
        IF (NUMPAI .NE. 1) CALL TIMPA(NOMPAI,NUMPAI,NEMDI,NKPART, ELX,
     &ELY,ERX1,ERY1,ERZ,BRX,BRY,BRZ,TRAJIO)
        ENDIF
C++
C
C  TEST FOR END OF LOOP
C
      IF(NT.LT.NEND) GO TO 100
C++
C++
C+++    ICI ON MEMORISE LES DONNEES CONCERNANT
C+++    CERTAINES PARTICULES
C+++     CHANGING PART

      IF (NUMPAE .NE. 1) WRITE(51) NUMPAE,NEMDI,NUMFIC,NENDO,NEND,DT,
     &TRAJEL
      IF (NUMPAI .NE. 1) WRITE(61) NUMPAE,NEMDI,NUMFIC,NENDO,NEND,DT,
     &TRAJIO 
C++     
      WRITE(15) NT,X,Y
      WRITE(15) PX,PY,PZ
      WRITE(15) PSA
      WRITE(16) NT,ETX,ETY,ETZ
      WRITE(16) BTX,BTY,BTZ
      WRITE(17) NC,EFL,EFT,BFT,EFIEL,EKINE
     1  ,EKINI,EKIN,ETOT
     2  ,EFLX,EFLXEN,EFLY,EFTX,EFTY,EFTZ
     3  ,BFTX,BFTY,BFTZ
     4  ,TEMPEX,TEMPEY,TEMPEZ,TEMPIX
     5  ,TEMPIY,TEMPIZ
     6  ,CUEX,CUEY,CUEZ,CUIX,CUIY,CUIZ
     7  ,TPLCE1,TPLCE2,TPALCE,TPLCI1,TPLCI2,TPALCI
      CLOSE(15)
      CLOSE(16)
      CLOSE(17)
C++
      PRINT*,'LE RUN A ETE SAUVEGARDE JUSQU''AU TEMPS T=',T
C++
C
C PLOT TIME HISTORIES
C
C++
C++  IF HISTO=.FALSE. THEN TIME HISTORIES ARE NOT PLOTTED
C++
      IF (HISTO) THEN
      CALL GRAF1(EFL,1,NC,'EFL',0)
      CALL GRAF1(EFT,1,NC,'EFT',0)
      CALL GRAF1(BFT,1,NC,'BFT',0)
C++
C++
      CALL GRAF1(EFLX,1,NC,'WELX',0)
C--   CALL GRAF1(EFLXEN,1,NC,'WELX-PIS',0)
      CALL GRAF1(EFLY,1,NC,'WELY',0)
      CALL GRAF1(EFTX,1,NC,'WETX',0)
      CALL GRAF1(EFTY,1,NC,'WETY',0)
      CALL GRAF1(EFTZ,1,NC,'WETZ',0)
      CALL GRAF1(BFTX,1,NC,'WBTX',0)
      CALL GRAF1(BFTY,1,NC,'WBTY',0)
      CALL GRAF1(BFTZ,1,NC,'WBTZ',0)
C++
C++
      CALL GRAF1(EFIEL,1,NC,'EFIEL',0)
      CALL GRAF1(EKINE,1,NC,'EKINE',0)
      CALL GRAF1(EKINI,1,NC,'EKINI',0)
      CALL GRAF1(EKIN,1,NC,'EKIN',0)
      CALL GRAF1(ETOT,1,NC,'ETOT',0)
      CALL GRAF1(TEMPEX,1,NC,'TE X',0)
      CALL GRAF1(TEMPEY,1,NC,'TE Y',0)
      CALL GRAF1(TEMPEZ,1,NC,'TE Z',0)
      CALL GRAF1(TEMPIX,1,NC,'TI X',0)
      CALL GRAF1(TEMPIY,1,NC,'TI Y',0)
      CALL GRAF1(TEMPIZ,1,NC,'TI Z',0)
      CALL GRAF1(TPLCE1,1,NC,'TPERPE O',0)
      CALL GRAF1(TPLCE2,1,NC,'TPERPE I',0)
      CALL GRAF1(TPALCE,1,NC,'TPARA E',0)
      CALL GRAF1(TPLCI1,1,NC,'TPERPI O',0)
      CALL GRAF1(TPLCI2,1,NC,'TPERPI I',0)
      CALL GRAF1(TPALCI,1,NC,'TPARA IO',0)
      CALL GRAF1(CUEX,1,NC,'CUE X',0)
      CALL GRAF1(CUEY,1,NC,'CUE Y',0)
      CALL GRAF1(CUEZ,1,NC,'CUE Z',0)
      CALL GRAF1(CUIX,1,NC,'CUI X',0)
      CALL GRAF1(CUIY,1,NC,'CUI Y',0)
      CALL GRAF1(CUIZ,1,NC,'CUI Z',0)
      ENDIF
C++
      CELEL=CELEL-1
      CELIO=CELIO-1
      PRINT*,'TEMPS EGAL A ',T
      CALL GRDRES(CELEL,CELIO)
C+
C+    FERMETURE NAG GRAPHIQUE
C+
      CALL CLSGKS

C  CALCULATE THE COMPUTING TIME REQUIRED.
C
      TIME2 = SECOND()
      TIME2=TIME2-TIME1
      WRITE(6,1002) TIME2
 1002 FORMAT(' TOTAL CPU TIME : ',E12.3,1X,'SECONDS.')
      PRINT *,'ON A STOCK{ ',NUMPAG,' PAGES DANS LE M{TAFICHIER'
      STOP
      END
