source: PSPA/parmelaPSPA/trunk/parmelamainv503.f @ 73

Last change on this file since 73 was 73, checked in by touze, 12 years ago

ajout pprincipal
File size: 44.7 KB
Line 
1c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
2c WARNING                                                              c
3c Warning for change of System and Operating system we have to do      c
4c one  change of the variable kvers in main program ligne 282          c
5c utilised in ddate                                                    c
6c----------------------------------------------------------------------c
7c                                                                      c
8c                            VERSION V5.03                             c
9c                                                                      c
10c     October 90  Orsay    V4.00                                       c
11c     Release Mai 92       V4.10                                       c
12c     Release June 92      V4.20                                       c
13c     Release October 92   V4.21                                       c
14c     Release December 92  V4.22                                       c
15c     Release July 93      V4.30                                       c
16c     Release March 94     V4.31                                       c
17c     Release April 94     V4.32                                       c
18c     Release June  95     V5.00                                       c
19c     Release June  96     V5.01                                       c
20c     Release July  96     V5.02    VMS / UNIX                         c
21c     Release January 97   V5.03    VMS / UNIX                         c
22c     Release September 2000 V5.03  Multi plateforme sans cernlib      c
23c     Release April 2010   V5.03    devient version officielle LAL     c
24c                                                                      c
25c                                                                      c
26c     B. Mouton  LAL Bat 200 91405 Orsay France                        c
27c                                                                      c
28c     MOUTON@LAL.IN2P3.FR                                              c
29c                                                                      c
30c     Version ALPHA                                                    c
31c     This version of PARMELA is a pachwork of 3 versions              c
32c     one of from Lloyd Young (version of July 88) of Los Alamos,      c
33c     another of KT Mc Donald of Brookhaven and one of SLAC with       c
34c     some local modifications, with modifications and news            c
35c     routines from H Liu of IHEP Beijing,with news routines from      c
36c     CEBAF.                                                           c
37c                                                                      c
38c---PhAse and Radial Motion in Electron Linear Accelerators            c
39c---(and associated transport systems)                                 c
40c                                                                      c
41c   on VAX VMS you have to increase your Paging file quota to 20000    c
42c                                                                      c
43c   $for/noop parmelav503                                              c
44c   $link/exe=parmelav5 parmvelav503,'lib$'                            c
45c                                                                      c
46c   # Makefile for PARMELA cas un seul fichier                         c
47c   # on Alpha  OSF1                                                   c
48c   #                                                                  c
49c   parmela : parmelav5.f                                              c
50c           f77 -g -O0 parmelav5.f -o parmelav5  ${CERNLIBS}           c
51c                                                                      c
52c   # on HP                                                            c
53c   parmela : parmelav5.f                                              c
54c     f77 +ppu -K -g -w +E1 -O parmelav503.f -o parmelav5hp ${CERNLIBS}c
55c                                                                      c
56c----------------------------------------------------------------------c
57c                                                                      c
58c   Pour le cas de la version eclatee 07/2000                          c
59c   #!//bin/ksh                                                        c
60c   # gmake sans cernlib                                               c
61c   programe="parmela"                                                 c
62c                                                                      c
63c   for files in *.f                                                   c
64c   do                                                                 c
65c   objname="${files:%.f}"                                             c
66c   echo $objname                                                      c
67c   objname=${objname}.o                                               c
68c   liste=${liste}" ${files}"                                          c
69c   #echo $liste                                                       c
70c   echo ${objname}:${files} >>temp1                                   c
71c   echo "\t \$(FF) -c \$(FFLAGS) "${files} >>temp1                    c
72c   echo " ">>temp1                                                    c
73c   done                                                               c
74c                                                                      c
75c   echo "${programe}: \$(OBJ)" >>temp1                                c
76c   #echo "\t \$(FF)  \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)"      c
77c          \${CERNLIBS} >>temp1                                        c
78c   echo "\t \$(FF)  \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)"       c
79c          >>temp1                                                     c
80c   echo " " >>temp1                                                   c
81c   echo "clean:" >>temp1                                              c
82c   echo "\t rm *.o">>temp1                                            c
83c                                                                      c
84c   # ecriture dans le fichier de makefile definitif                   c
85c                                                                      c
86c   ## Pour HP                                                         c
87c   echo FF=f77 >>makefile                                             c
88c   #echo FFLAGS=+ppu +U77 -K -g -w +E1 -O >>makefile                  c
89c   echo FFLAGS= -K -w +E1 +Oall >>makefile                            c
90c   echo " " >>makefile                                                c
91c   echo "VER=alonehpopt" >>makefile                                   c
92c                                                                      c
93c   ## Pour OSF                                                        c
94c   #echo FF=f77 >>makefile                                            c
95c   #echo FFLAGS= -g -O2 >>makefile                                    c
96c                                                                      c
97c   echo " " >>makefile                                                c
98c   echo SRC=${liste} >>makefile                                       c
99c   echo " " >>makefile                                                c
100c   echo "OBJ=\$(SRC:.f=.o)" >>makefile                                c
101c   echo " " >>makefile                                                c
102c   cat temp1 >> makefile                                              c
103c   rm -f temp1                                                        c
104c   #c'est fini                                                        c
105c----------------------------------------------------------------------c
106c                                                                      c
107c          modification des call scheff pour input0 02/92              c
108c          modification dans emit90 pour kk=0 05/92 --> version 4.10   c
109c          modification dans elips pour negative square root 05/92     c
110c          outlal, emit90 and elips are in double precision   05/92    c
111c          out6 is in double precision 06/92                           c
112c          version 4.20  02/06/92                                      c
113c          modification dans out6dp boucle 9996 remplace cor par cord  c
114c          version 4.21                                                c
115c          version 4.22 ip.eq.999 for imp.dat cell and trwave          c
116c          modification dans cellfeld rcir en rcur  04/93              c
117c          version 4.30 suppression of calls hbook                     c
118c          version 4.31 24/03/94                                       c
119c          pvol is in double precision pvoldp 03/94                    c
120c          pvoldp supressed icall 03/94                                c
121c          intro de cord(6,1)=cor(6,1) in main after 125 continue 03/94c
122c          scheff2 loop 20 i=2,nr ----> i=2,nr1 03/94                  c
123c          in out6dp change in format 220 zmax = f8.4 in               c
124c          zmax = f8.2 03/94                                           c
125c          version 4.32 18/04/94                                       c
126c          introduction scaling factor on foclal card 04/94            c
127c          modif /pcord/ intro nupar 05/94                             c
128c          nspl=10000  ----> 10001          11/94                      c
129c          nsplc=1497 ----> 30005=10001*3+2 11/94                      c
130c          in zoutlal replaced in data bcd : design by trwave 11/94    c
131c          renumber subroutine input 27/03/95                          c
132c          input 11 27/03/95                                           c
133c          histogram1 and histogram2 28/03/95                          c
134c          renumber main program 29/03/95                              c
135c          subroutine poisson 30/03/95                                 c
136c          input 12 31/03/95                                           c
137c          subroutine sextupole 11/04/95                               c
138c          modification of space charge card 11/04/95                  c
139c               new parameter iprog ;                                  c
140c               scheff beami iprog point and subroutine dependent      c
141c          modification extension infld in sffld 01/06/95              c
142c          zout introduction de angle pour signe arctangente 04/09/95  c
143c          modification input11 write ngood et tmax=2.5 sigma 12/09/95 c
144c          modif dans input11 du calcul de degmax                      c
145c                         remplace tmax par tfac   28/09/95            c
146c          intro dans parmela routine gpindp plus dans cernlib 28/09/95c
147c          deplacement de l'open de parmtrch cas ip=3 19/02/96         c
148c          modification input 4 26/06/96                               c
149c                               - 8 parameters one first card input 4  c
150c          modification input11 05/07/96 two news parameters angle of  c
151c            injection of laser light and distane focal point cathode  c
152c          modification input10 23/07/96 two news parameters angle of  c
153c            injection of laser light and distane focal point cathode  c
154c          13/12/96 modif zout pour trwave (phase)                     c
155c          13/12/96 modif fieldlal pour nchp.lt.nspl                   c
156c          13/01/97 modif fieldlal pour free format                    c
157c          13/01/97 modif common/image/                                c
158c          25/02/97 intro dans subroutine input pour input15           c
159c            dppi,dpepsout double precision et modif de l'appel gpindp c
160c          03/09/97 mise sous Unix des statistiques elapsed cpu time   c
161c          04/09/97 version pour HP (c_hp)                             c
162c          14/11/97 nettoyage des variables non utilisees              c
163c                   intro du call date_and_time sous VMS pour an 2000  c
164c                   tous les open avec status='scratch' sans filename  c
165c          03/02/98 modifation de zout et zoutlal pour cas solenoid    c
166c          18/09/98 introduction de iadp a la place ngood dans out6dp  c
167c                   write(ndes2), write(ndes1)                         c
168c          29/03/00 reprise version eclate a partir du source non      c
169c                   non eclate du 31/05/99  pour correction d'un bug   c
170c                   dans une version ephemere dite v503b               c
171c          18/07/00 abandon de la cernlib                              c
172c                   suppression de call histo dans main et input       c
173c          20/07/00 Modification de ddate.f pour essayer une methode   c
174c                   universelle de date, modification du main          c
175c                   introduction de kvers, introduction de la          c
176c                   subroutine datime.f machine dependante             c
177c                   suppression de l'utilisation de la cernlib         c
178c          13/09/00 Introduction des include 'filename.h'              c
179c                   Les fichiers .h  contiennent les parameters et les c
180c                   common.                                            c
181c                   Modification des noms de variable dans scwall      c
182c                   Modification npoinz en npoints dans scheff2        c
183c                   Declaration real zloc, gam, beami, brhof, radian , c
184c                   twopi, z0 dans out6dp                              c
185c                   Dans les appels a scheff ajouter 2 variables mt    c
186c                        et wtsauv                                     c
187c                   Modification de swap pour /pcord/                  c
188c                   Suppression des instructions                       c
189c                      dimension ac(14) cc(9) bc(1+(1+lmx)*13          c
190c                      equivalence (ac(1),pi) (bc(1),nel) (cc(1),wt0)  c
191c                   Ajouter mt aux arguments de la subroutine output   c
192c                   Debug version modif main pardyn ouvrir trwave      c
193c                   new common debug , new possibility for ip          c
194c                   negative value means keep file parmdebug           c
195c                   common/debug/idebug et tous les idebug             c
196c    octobre 2002   augmentation du nombre de cartes cell numc 20-->60 c
197c                                                                      c
198c    29/10/2002     nom de l'exe sous OSF et HP-UX : parmela           c
199c                                                                      c
200c            nom de l'exe sous OSF et HP-UX : parmela                  c
201c                                                                      c
202c    19/04/2007    introduction de la polarisation ksi1,ksi2,ksi3      c
203c                  modification lecture input4                         c
204c                  modification out6dp  et sortie                      c
205c                  18 variables dans ndes1 et ndes2 dans out6          c
206c                                                                      c
207c Last revision 14/11/2007 mise sous Linux                             c
208c             utilisation du debugger                                  c
209c             gdb /exp/sera/mouton/Parmela/V503/Linux/parmelalinux     c
210c             pour cela faire cp toto.in parmin avant de lancer gdb    c
211c                                                                      c
212c    14/11/2007     modification des open de file.f, pour utilisation  c
213c                   du disque local de la machine /tmp                 c
214c                   dans le parmela.sh il faut :                       c
215c                                      mv /tmp/parmdesz toto.desz      c
216c Last revision 18/06/2008                                             c
217c             un seul source de parmela v503   compatible OSF LINUX    c
218c             un seul nom d'exe parmela                                c
219c                                                                      c
220c     Release octobre 2009  V5.03                                      c
221c         modification du calcul de dzz dans fieldlal                  c
222c              dzz=z(2)-z(1) ----> dzz=(zmax-zmin)/(nchp-1)            c
223c         augmentation du nombre de points pour le  champ magnetique   c
224c         nspl 1001 ---> 10001                                         c
225c         nouveau parametre nptcb = 10000 dans param_sz.h              c
226c         creation d'une nouvelle carte nfieldb avec comme parametre   c
227c          nptcbu nombre de points que l'on veut utiliser pour le      c
228c          champ magnetique < ou = 10000                               c
229c    04/02/2010 retrait de la carte nfieldb                            c
230c               modification de la 1er carte coil                      c
231c               des tests dans field, solnoid si nb pt champ b > nptcb c
232c                                                                      c
233c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
234c-----------------------------------------------------------------------
235c---->program parmela
236      subroutine parmela()
237c-----------------------------------------------------------------------
238c label n --> n*100
239c go to  7000 --> 7008
240c do     8000 --> 8017
241c format 9000 --> 9025
242c-----------------------------------------------------------------------
243c
244      include 'param_sz.h'
245      include 'var_char.h'
246      include 'coordcom.h'
247      include 'constcom.h'
248      include 'debug.h'
249      include 'bfieldcom.h'
250      include 'flagcom.h'
251      include 'misccom.h'
252      include 'ncordscom.h'
253      include 'outcom.h'
254      include 'phizoutcom.h'
255      include 'pipecom.h'
256      include 'sccom.h'
257      include 'syscom.h'
258      include 'tstepcom.h'
259      include 'ucom.h'
260      include 'upawcom.h'
261      character*25 kvers
262      common/aswap/acor(7,imaa)
263      common/back/bmax,byz0,byzs,pzowmin,pzowmax,iback
264      common/cathcur/rcur
265      common/chars/title
266      common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7
267      common/ct/var,rp,dzmin,dzmax,wrm,dwbm,xrmsbg,yrmsbg,
268     1xsmax,xpsmax,ysmax,ypsmax,istat,nebm
269      common/ddat/ctime,ddat
270      common/errors/eraln(5,0:lmx)
271      common/image/jj,ij
272      common/intype/intype ! mcd
273      common/jones/ajones,zjones,zcath(imaa) ! mcd
274      common/out6/nbphase,nbremesh
275      common/outbuf/outcor(8,imaa*imb)
276      common/outpar/noutp,outpar(0:lmx,10),moutp(10)
277      common/part/nparticle(imaa)
278      common/pnum/ipnum(imaa)
279      common/save/zrsav
280      common/sbload1/tsbparm(5,50)
281      common/tapes/ifl(2)
282      common/version/kvers
283      common /wtstep/wtstep,nstep ! mcd
284c--------------------------------------------------------------------------
285      character Istart_date*9,Iend_date*9,Istart_time*8,Iend_time*8
286      real*4    time_begin
287
288c_with cernlib
289c      call timest(1.0e9)
290c      call timed(t0)
291c
292
293      write(*,*) ' Go Parmela'
294c_for ibm
295c      time_begin = secnds(0.0)
296
297      call datime(Istart_date,Istart_time)
298c*
299c---set constants
300c
301       kvers='V503 on LINUX'
302c       kvers='V503 on OSF'
303c
304      pi=4.*atan(1.)
305      twopi=2.*pi
306      radian=pi/180.
307      clight=29979.29   ! 10.e6 cm/s
308      erest=.511        ! MeV
309      brhof=1704.51     ! gauss.cm moc/e
310c     brhof = (9.109389E-31*299792458/1.602177E-19)*(100)*(10000) !1704.5
311      zlimit=1000000.   ! cm
312c---set flags
313      idebug=0          ! for debuging in subroutine ouvrir
314      ifld=0            ! flag for external electric field
315      poiflag=.false.   ! flag for external magnetic field
316      nel=0             ! number of element
317      ip=0              ! flag of output type
318      isc=0             ! flag for space charge calculation
319      istat=0           ! flag for statistics
320      ifl(1)=0
321      ifl(2)=0
322      iback=0
323      ltype=1           ! by default disk and washer cells linac
324      ifoclal=0         ! bm bzlal
325      npipe=0           ! mcd
326      nstep=0           ! mcd
327      wtstep=0.         ! mcd
328      ajones=0.         ! mcd
329      nout=0            ! mcd
330      do 8001 n=1,10    ! mcd
331      moutp(n)=0        ! mcd
332      do 8001 i=1,lmx   ! mcd
333      outpar(i,n)=0.    ! mcd
334 8001 continue          ! mcd
335      icfield = 0 ! cbm
336      ictrw= 0    ! cbm
337      nbphase=0   ! cbm nombre de phase en sortie
338      nbremesh=0  ! cbm nombre de remesh
339c----set polarization to zero
340      do 7999 n=1,imaa
341         ksi1(n)=0.
342         ksi2(n)=0.
343         ksi3(n)=0.
344 7999 continue
345c---
346c---
347      open(unit=nin,file='parmin',status='old',
348     *     access='sequential',form='formatted')
349      open(unit=nnout,file='parmdebug',status='unknown',
350     *        access='sequential',form='formatted')
351      call ouvrir
352c---
353      call setwall
354c---set misalignments to zero.
355      do 8002 n=1,100
356      do 8002 i=1,5
357      eraln(i,n)=0.
358 8002 continue
359c---read input data
360c---the read(jj)  reads one card with a label and all
361c---subsequent cards without labels.  the stop subroutine compares the label
362c---with a list of words and assigns a value to jj corresponding to the
363c---position of the word in the dictionary.  the word list and corresponding
364c---values of jj are:
365c---(drift     ,1), (solenoid  ,2), (quad      ,3), (bend      ,4),
366c---(buncher   ,5), (chopper   ,6), (cell      ,7), (tank      ,8),
367c---(trwave    ,9), (coil      ,10),(run       ,11),(input     ,12),
368c---(output    ,13),(title     ,14),(scheff    ,15),(zout      ,16),
369c---(adjust    ,17),(start     ,18),(restart   ,19),(continue  ,20),
370c---(save      ,21),(end       ,22),(zlimit    ,23),(errors    ,24),
371c---(change    ,25),(rotate    ,26),(sbload    ,27),(cfield    ,28),
372c---(dpout     ,29),(cathode   ,30),(design    ,31),(pipe      ,32),
373c---(foclal    ,33),(backb     ,34),(wiggler   ,35),(alpham    ,36),
374c---(stat      ,37),(poisson   ,38),(sextupole ,39)
375c---the subroutine reads numbers from the card in a free format.
376c---the values are stored in vv, and there are nn of them.
377 99   continue
378cbm      call cread(jj)
379      call cread
380      go to (100,100,100,100,100,100,100,100,100,1000,
381     1       1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,
382     2       2100,2200,2300,2400,2500,100,100,2800,100,100,
383     3       3100,3200,3300,3400,100,100,3700,3800,100),jj
384c---after transport element or linac cell
385c---after drift:     vv(n)=l,aper,iout
386c---after solenoid:  vv(n)=l,aper,iout,h
387c---after quad:      vv(n)=l,aper,iout,hp
388c---after chopper:   vv(n)=0,aper,iout,freq,phi,dphi
389c---after cathode:   vv(n)=0,aper,iout,radius of curvate of cathode
390c---after wiggler:   vv(n)=l,aper,iout,periods,steps,b0,ky/kw
391c---after sextupole: vv(n)=l,aper,iout,Bpole/ape**2
392 100  continue
393      if(nn.gt.11.and.jj.ne.7.and.jj.ne.9.and.jj.ne.27)go to 7000
394      if(nel.ge.lmx)go to 7001
395      if(nn.lt.2)go to 7002
396      nel=nel+1
397      ntype(nel)=jj
398      do 8003 n=1,nn
399      el(n,nel)=0.      ! mcd
400      if(n.gt.nn)go to 8003
401      el(n,nel)=vv(n)
402 8003 continue
403      zloc(nel)=el(1,nel)
404      if(nel.gt.1)zloc(nel)=zloc(nel)+zloc(nel-1)
405      if (ip.eq.1 .or. ip.eq.2) then
406        iout = vv(3)
407        write(nnout,9001) nel,jj,vv(1),zloc(nel),vv(2),iout
408 9001   format (
409     . ' element no. [ne] =',t40,i8/
410     . ' element type [ntype(ne)] =',t40,i8/
411     . ' element length [el(1,ne)] =',t40,1pg12.4,' cm'/
412     . ' z coord at downstream end [zloc(ne)] =',t40,g12.4,' cm'/
413     . ' radial aperture at exit [el(2,ne)] =',t40,g12.4,' cm'/
414     . ' output at end of element [el(3,ne)] =',t40,i8)
415        if((jj.eq.3) .or. (jj.eq.39)) write(nnout,9002) vv(4)
416 9002   format(
417     . ' magnetic field gradient = ',t40,g12.3,' gauss/cm')
418      endif
419      if(jj.eq.2) ifld=1
420      if(jj.eq.26) go to 2600
421      if(jj.eq.29) go to 99
422      if(jj.eq.27) go to 2700
423      if(jj.eq.30) rcur=vv(4)
424      if(jj.eq.35) go to 3500
425      if(jj.eq.36) go to 3600
426      if(jj.eq.39) go to 99
427      go to (99,99,99,400,500,99,700,800,900), jj
428c---after bend:  vv(n)=l,aper,iout,wr,alpha,beta1,beta2,psi1,psi2,r1,r2
429  400 continue
430      do 8004 n=5,9
431      el(n,nel)=el(n,nel)*radian
432 8004 continue
433      go to 99
434c---after buncher:  vv(n)=0.,aper,iout,wmax,freq,phi,w
435  500 continue
436      w=w0
437      if(nn.ge.7)w=el(7,nel)
438      g=w/erest
439      bg=sqrt(g*(2.+g))
440      cay=twopi*el(5,nel)/(bg*clight)
441      el(8,nel)=cay
442      go to 99
443c---after trwave;
444  900 continue
445      ictrw = ictrw+1
446      call trwfield
447      go to 99
448c---after cell:  vv(n)=cl,aper,iout,phi,e0,nc,(dwtmax)
449  700 continue
450      el(5,nel)=.01*el(5,nel)
451      if (nn.lt.8) el(8,nel)=0.
452      phideg(nel)=el(4,nel)
453c      phi=el(4,nel)*radian this is now defined in cellfld
454c      el(9,nel)=sin(phi)    and this
455c      el(10,nel)=cos(phi)  and this also
456c---if unspecified, set maximum step size in cell to 10 degrees.
457      if(nn.lt.7)el(7,nel)=10.
458c---get field information for this cell
459      cl=el(1,nel)
460      if (el(8,nel).ne.0.) cl=2.*cl
461      nc=el(6,nel)
462c      iel(6,nel)=nc
463      if (ip.eq.1 .or. ip.eq.2) then
464        write(nnout,9003)vv(4),vv(5),nc,el(7,nel),el(11,nel),el(8,nel)
465 9003   format (
466     . ' phase of cell [phi0] =',t40,g12.3,' degrees'/
467     . ' average axial electric field [el(5,ne)] =',t40,g12.3,' MV/m'/
468     . ' cell no. [nc] =',t40,i8/
469     . ' max step size [dwtmax] =',t40,g12.3,' degrees'/
470     . ' axial static magnetic field =',t40,g12.3,' gauss'/
471     . ' cell symmetry [el(8,ne)] =',t40,g12.3)
472        if (el(8,nel).eq.0.) write(nnout,9004)
473 9004   format (' cell is symmetric about midplane in z')
474        if (el(8,nel).gt.0.) write(nnout,9005)
475 9005   format (' half cell, with zero field at entrance')
476        if (el(8,nel).lt.0.) write(nnout,9006)
477 9006   format (' half cell, with zero field at exit')
478      endif
479      call cellfld(cl,nc)
480      go to 99
481c---after tank:  tl,aper,iout,e0t,nc,p,phi
482  800 continue
483      el(8,nel)=el(1,nel)/el(5,nel)
484c---convert e0t from mv/m to mv/cm
485      el(4,nel)=.01*el(4,nel)
486c---calculate beta and k**2 for the tank
487      beta=2.*el(8,nel)/(el(6,nel)*wavel)
488      cl=el(8,nel)
489      caysq=((wavel/(2.*cl))**2 - 1.)*(twopi/wavel)**2
490      el(9,nel)=beta
491      el(10,nel)=caysq
492      go to 99
493c---after field:  vv parameters depend on field subroutine
494 1000 continue
495      icfield=icfield+1
496      call field
497c---define field structure after coil
498c--- vv(n)=z,r,ni,(zbeg,zend,nptcbu,order)
499c---in ( ) required on first coil card ignored on rest
500      go to 99
501c---after errors.  generate errors in geners
502 2400 continue
503      call geners
504      go to 99
505c---after run:  vv(n)=irun,ip,freq,z0,w0,ltype
506 1100 continue
507c     process run card         ! mcd
508      if (nn.ge.1) irun=vv(1)  ! mcd
509c--debug
510c-- Si ip negatif sur carte run on a la version de debug
511      if (nn.ge.2) then
512         if(vv(2).lt.0.) then
513            idebug=1
514            write(*,*) '-----------------------------------'
515            write(*,*) 'You have ask the version with debug'
516            write(*,*) '-----------------------------------'
517            write(nnout,*) '-----------------------------------'
518            write(nnout,*) 'You have ask the version with debug'
519            write(nnout,*) '-----------------------------------'
520         else
521         endif
522         ip=abs(vv(2))
523      endif
524c--debug
525      if (nn.ge.3) freq=vv(3)  ! mcd
526      if (nn.lt.4) go to 99 
527      if (nn.ge.4) z0=vv(4)    ! mcd
528      if (nn.lt.5) go to 99
529      if (nn.ge.5) w0=vv(5)    ! mcd
530      if (nn.ge.6) ltype=vv(6) ! mcd
531      if(freq.gt.0.) wavel=clight/freq
532      write(nnout,9007)irun
533 9007 format (/' '/
534     . ' run no. =',t40,i8/)
535      if((ip.gt.4).and.(ip.ne.999)) then
536      ip=1
537      write(nnout,9008)ip
538 9008 format(
539     .' print option modified =',t40,i8)
540      else
541      write(nnout,9009)ip
542 9009 format(
543     . ' print option =',t40,i8)
544      endif
545c---
546      if (ip.eq.3) then
547      open(unit=nuch,file='parmtrch',access='sequential',
548     *     status='unknown')
549      write(nuch,*)'  nc  sym  r  z  ez  er  bt  '
550      endif
551      if (ip.eq.999) open(unit=nimp,file='parmimp',status='unknown')
552c---
553      write(nnout,9010)freq,wavel,z0,w0,ltype
554 9010 format(
555     . ' rf frequency =',t40,f9.2,'    MHz'/
556     . ' wavelength =',t40,1pg12.3,' cm'/
557     , ' z0 =',t40,g12.3,' cm'/
558     , ' reference particle energy =',t40,g12.3,' MeV'/
559     . ' linac type =',t40,i8)
560      do 8005 i=1,4
561      corout7(i,1)=0.
562      cor(i,1)=0.
563 8005 continue
564      cor(5,1)=z0
565      corout7(5,1)=z0
566      corout7(6,1)=w0
567      g=w0/erest
568      gam(1)=g+1
569      cor(6,1)=sqrt(g*(2.+g))
570      cord(6,1)=cor(6,1)
571      if (npoints.lt.1)npoints=1
572      ngood=npoints
573      nparticle(1)=1
574      phizero(1)=0.
575      if(nn.lt.6)go to 99
576      go to 99
577c---after input:  vv parameters depend on input subroutine
578 1200 continue
579      call input
580      ngood=npoints
581      do 8006 i=1,ngood
582      nparticle(i)=i
583      ipnum(i)=i
584 8006 continue
585      go to 99
586c---after output:  vv parameters stored in optcon array
587c---after output:  vv parameters stored in outpar array
588c
589 1300 continue
590      noutp=vv(1)+1
591      if(noutp.gt.10 .or. noutp.lt.1) then
592      write(ndiag,9011)
593 9011 format(' invalid OUTPUT type')
594      go to 99
595      endif
596      moutp(noutp)=nn
597      do 8007 n=1,nn
598      optcon(n)=vv(n)
599      outpar(n,noutp)=vv(n)
600 8007  continue
601      if(vv(1).gt.1..and.vv(1).lt.5.)then
602      imbf=imb-1
603      else
604      imbf=imb
605      endif
606      go to 99
607c---after title:  read title information from next card
608 1400 continue
609      read(nin,9012)title
610 9012 format(a80)
611      write(nnout,*)title
612      go to 99
613c---after scheff:  vv(1)=beam current; vv(2)= program number
614c--- rest depend on space-charge subroutines
615 1500 continue
616      do 8008 n=1,12
617      sce(n)=0.
618      if (n.le.nn) sce(n)=vv(n)
619 8008 continue
620      beami=vv(1)
621      iprog=vv(2)
622      if(isc.eq.0) then
623        if(iprog.eq.1) call scheff1(0.,-99,0.)
624        if(iprog.eq.2) call scheff2(0.,-99,0.)
625        if(iprog.eq.3) call scheff3(0.,-99,0.)
626         isc=1
627      else          ! isc
628      endif         ! isc
629      go to 99
630c---after zout.  print locations of elements
631 1600 continue
632      if(ifoclal.eq.1) then
633        call zoutlal
634      else
635        call zout
636      endif
637      go to 99
638c---after adjust:  call arbitrary subroutine named adjust
639 1700 continue
640      call adjust
641      go to 99
642c---after change.  call change subroutine
643 2500 continue
644      call change
645      go to 99
646c---after rotate
647 2600 continue
648      el(4,nel)=el(4,nel)*radian
649      go to 99
650c---  after sbload paramaters are l,aperature,iout,zlen,ncells,a1,a2,a3
651c---  ,a4,a5,ta1,ta2,ta3,ta4,ta5. if ncells is negative then the paramaters
652c---  ta1 to ta4 are required. ta5 is optional. if ncells is negative the
653c---  transvers wake will be applied also.
654 2700 continue
655      if(nn.gt.11)then
656      numsb=numsb+1
657      tsbparm(1,numsb)=vv(11)
658      el(11,nel)=numsb
659      tsbparm(2,numsb)=vv(12)
660      tsbparm(3,numsb)=vv(13)
661      tsbparm(4,numsb)=vv(14)
662      if(el(5,nel).lt.0..and.nn.lt.14)go to 7002
663      if(nn.gt.14)tsbparm(5,numsb)=vv(15)
664      endif
665      go to 99
666c---after start: start dynamics calculations
667 1800 continue
668      if(nn.lt.5)go to 7002
669      wt0=vv(1)
670c---  For input type 10 and 11, WTSTEP is the half width of the laser pulse.
671c     The clock must be set back by this amount to generate the time
672c     distribution of photoelectrons.
673      wt0=wt0-wtstep
674      dwt=vv(2)
675      nsteps=vv(3)
676      nsc=vv(4)
677      nout=vv(5)
678      go to 7003
679c---after restart.  set the starting time (phase), initialize cord
680c---from cor, and find largest element number (len) that any particle is in.
681 1900 if(nn.lt.4)go to 7002
682      dwt=vv(1)
683      nsteps=vv(2)
684      nsc=vv(3)
685      nout=vv(4)
686c       -----------------
687      if(jj.eq.19) then                         
688                npoints=-vv(6)
689      endif
690c       -------------------------
691      if(nn.gt.5.and.(vv(6).gt.0..and.vv(6).le.26))then
692        i=vv(6)-1
693      call orest1(i)
694      else
695      call orest2
696      endif
697      read(nsav) cor,wt0,ngs
698      close(nsav)
699      ngood=ngs
700C       ----------------------------
701      IF(NGS.NE.NPOINTS) then
702        call appendparm
703        print *,'ngs=',ngs,'npoints=',npoints
704        stop ' Abnormal stop ngs.ne.npoints main '
705      endif
706C       -------------------------------
707 7003 continue
708      if (ip.eq.1 .or. ip.eq.2)
709     .   write (nnout,9013) wt0,dwt,nsteps,vv(5),nsc,nout
710 9013 format (
711     . ' starting phase =',t40,1pg12.3,' degrees'/
712     . ' integration step size =',t40,g12.3,' degrees'/
713     . ' no. of steps in integration =',t40,i8/
714     . ' shift of reference particle =',t40,g12.3,' cm'/
715     . ' call SCHEFF every',t40,i8,'     step(s)'/
716     . ' call OUTPUT every',t40,i8,'     step(s)'/)
717      wt=wt0
718      do 8009 i=1,maxbuf
719      neb(i)=0
720      nbuf(i)=0
721 8009 continue
722c---if nn.gt.4, the 5th parameter is dist to shift ref particle
723c   only for restart. not for start
724      if (jj.eq.19.and.nn.gt.4) then
725        cor(5,1)=cor(5,1)+vv(5)
726      endif
727c       --------------------------------------------------
728c--------- vv(5)= dz
729      if(jj.eq.19.and.vv(5).ne.0.) then
730            cor(1,1)=0.  ! for energy change of ref part
731            cor(2,1)=0.
732            cor(3,1)=0.
733            cor(4,1)=0.
734            print *,'vv(7) =',vv(7)
735            cor(6,1)=sqrt((1.+vv(7)/0.511)**2-1.)
736      endif
737c       ---------------------------------------------------
738      do 8010 n=1,ngood
739      do 8011 i=1,6
740      cord(i,n)=cor(i,n)
741 8011  continue
742      gam(n)=sqrt(1.+cord(2,n)**2+cord(4,n)**2+cord(6,n)**2)
743c---find which element each particle is in
744      ne=0
745      z=cord(5,n)
746      if(z.lt.0.)go to 7004
747      do 8012 i=1,nel
748      ne=i
749      if(z.lt.zloc(i))go to 7004
750 8012 continue
751      ne=nel+1
752 7004 cord(7,n)=ne
753c     if generating from a photocathode, put ne=0
754      if(jj.eq.19) go to 8010           
755      if (intype.eq.10 .or. intype.eq.11 .or. intype.eq.15) 
756     .   cord(7,n)=0.
757 8010 continue
758c---if nn.gt.4, the 5th parameter is dist to shift ref particle
759c   only for restart. not for start
760c---start dynamics calculation
761      call pardyn
762      go to 99
763c---after continue:  used to change step size and continue computations
764 2000 continue
765      if(nn.lt.4)go to 7002
766      dwt=vv(1)
767      nsteps=vv(2)
768      nsc=vv(3)
769      nout=vv(4)
770c---if nn.gt.4, the 5th parameter is dist to shift ref particle
771      if (nn.gt.4)cord(5,1)=cord(5,1)+vv(5)
772      call pardyn
773      go to 99
774c---after save
775c---save cord in cor
776 2100 continue
777      ngs=ngood
778      do 8013 n=1,ngood
779      do 8013 i=1,6
780      cor(i,n)=cord(i,n)
781 8013 continue
782      write(nnout,*) 'zr of save ',zrsav
783      wt0=wt
784      if(nn.ge.1.and.(vv(1).gt.0..and.vv(1).le.26))then
785      i=vv(1)-1
786      call osave1(i)
787      if(ifl(2).ne.0) then
788      close(nsav)
789      ifl(2)=0
790      endif
791      if(idebug.eq.0) then
792        close(unit=nnout,status='delete')
793      else
794        close(unit=nnout,status='keep')
795      endif
796      close(nnout)
797      open(unit=nnout,file='parmout',access='sequential',
798     *     status='unknown',form='formatted')
799 7005 continue
800      read(nnout,'(10a8)',err=7006)tem
801      go to 7005
802 7006 continue
803      backspace (nnout)
804      else
805      open (unit=nsav,file='savecor',access='sequential',
806     *      form='unformatted',status='new')
807      endif
808      write(nsav) cor,wt0,ngs,bgs
809      close(nsav)
810      go to 99
811c---after zlimit:  lower limit on lagging particles
812 2300 continue
813      zlimit=vv(1)
814      go to 99
815c---after cfield:  read fields from sfo7 as written on tape7
816c--- first parameter on cfield is starting cell number
817c--- next parameter is ending cell number and third parameter is
818c--- increment of cell number as in a do loop
819c---line after cfield is filename of sfo7 data
820 2800 continue
821C modif glm 29/10/2012      read(nin,'(a12)')filesuper
822      read(nin,*) filesuper
823      write(nnout,'(1x,a12)')filesuper
824      if(nn.lt.1)then
825      nc=1
826      ne=1
827      ns=1
828      elseif(nn.eq.1)then
829      nc=vv(1)
830      ne=nc
831      ns=1
832      elseif(nn.eq.2)then
833      nc=vv(1)
834      ne=vv(2)
835      ns=1
836      elseif(nn.ge.3)then
837      nc=vv(1)
838      ne=vv(2)
839      ns=vv(3)
840      endif
841      call cellfeld(filesuper,nc,ne,ns)
842      go to 99
843c---after design
844 3100 continue
845      call design
846      go to 99
847c---after pipe, vv(n) = radius,zlow,zhigh
848 3200 continue
849      if (npipe.eq.100) then
850        write(nnout,9014)
851 9014   format (' Too many PIPE cards')
852        go to 99
853      endif
854      if (vv(1).le.0.) then
855        write(nnout,9015)
856 9015   format (' Bad radius on PIPE card')
857        go to 99
858      endif
859      npipe = npipe + 1
860      rpipe(npipe) = vv(1)
861      zlpipe(npipe) = vv(2)
862      zhpipe(npipe) = vv(3)
863      if (ip.eq.1 .or. ip.eq.2)
864     .write(nnout,9016) npipe,vv(1),vv(2),vv(3)
865 9016 format (
866     . ' Beam pipe for point-by-point image calculation'/
867     . ' pipe no. = ',t40,i8/
868     . ' pipe radius = ',t40,1pg12.3,' cm'/
869     . ' zlow = ',t40,g12.3,' cm'/
870     . ' zhigh = ',t40,g12.3,' cm')
871      go to 99
872c---after foclal vv(n)=zmin,zmax,dzz,npchb,facbz,opt
873 3300 continue
874      iopt=vv(6)
875      if(iopt.eq.0) then
876         ifld=0
877         ifoclal=0
878         read(nin,'(a12)') filebz
879         write(nnout,*) 'foclal not used'
880      else
881         ifld=1
882         ifoclal=1
883         read(nin,'(a12)') filebz
884         call fieldlal(filebz)
885      endif
886      go to 99
887c---after backb back bombardment
888 3400 continue
889      iback=1
890      bmax=vv(1)
891      byz0=vv(2)
892      byzs=vv(3)
893      iope=0
894      if(nn.gt.3) then
895           pzowmin=vv(4)
896           pzowmax=vv(5)
897      else
898           pzowmin=0.
899           pzowmax=0.
900      endif
901      if(pzowmax.ne.0.) iope=1
902      call opeback(iope)
903      go to 99
904c---after wiggler l,aper,iout,periods,steps,b0,kx,ky,kw,iprint
905 3500 continue
906      if(nn.lt.7) go to 7002
907      xlam=el(1,nel)/el(4,nel)
908      cayw=twopi/xlam
909      cayy=el(7,nel)*cayw
910      if(el(7,nel).le.0. .or. el(7,nel).gt.1.)then
911         write(nnout,9018)
912 9018    format(' error: ky/kw on wiggler must be between 0 and 1')
913         jj=100 ! for cread
914         go to 99
915      endif
916      cayx=0.
917      cayxsq=cayw**2-cayy**2
918      if(cayxsq.gt.0.) cayx=sqrt(cayxsq)
919      el(7,nel)=cayx
920      el(8,nel)=cayy
921      el(9,nel)=cayw
922      if(ip.eq.1 .or. ip.eq.2)
923     .write(nnout,9019)el(4,nel),el(5,nel),el(6,nel),el(7,nel)
924 9019 format(
925     .' number of wiggler periods =',t40,g12.3/
926     .' step of wiggler =',t40,g12.3/
927     .' wiggler field =',t40,g12.3,' gauss'/
928     .' ratio of Ky to Kw =',t40,g12.3)
929      go to 99
930c---after alpha magnet
931 3600 continue
932      if(ip.eq.1 .or. ip.eq.2)
933     .write(nnout,9020) el(1,nel),el(4,nel),el(5,nel),el(7,nel)
934 9020 format(
935     .' distance travelled by the r.p. =',t40,g12.3,' cm'/
936     .' reference particle energy =',t40,g12.3,' MeV'/
937     .' alpha magnet angle =',t40,g12.3,' Degrees'/
938     .' absolute field gradient =',t40,g12.3,' Gauss/cm')
939      go to 99
940c---after stat
941 3700 continue
942      var=vv(1)
943      istat=1
944      open(unit=nstat,file='parmstat',form='unformatted',
945     *     status='unknown')
946      write(nstat)ij
947      go to 99
948c---after Poisson
949 3800 continue
950C modif glm 29/10/2012      read (nin,'(a12)')filepois
951      read (nin,*)filepois
952      write(nnout,'(1x,a12)')filepois
953      if(nn.ge.1)then
954      zoffset=vv(1)
955      else
956      zoffset=0.
957      endif
958      if(nn.ge.2)then
959      rmult=vv(2)
960      else
961      rmult=1.
962      endif
963      call poisson(filepois,zoffset,rmult)
964      go to 99
965c
966cbm retrait de nfieldb 04/02/2010 apres modif de coil
967c--- after nfieldb number of points used for the magnetic field
968c    < 10000
969c 4000 continue
970c      nptcbu=vv(1)
971c      write(nnout,*)
972c     *' Nombre de points maximal pour le champ magnetique ',nptcb
973c      write(nnout,*)
974c     *' Number of used values for magnetic field',nptcbu
975c      if(nptcbu.gt.nptcb) then
976c        write(nnout,*)
977c     *  ' Number of points for magnetic field gt ',nptcb
978c      endif
979c      go to 99 
980c---after-end output any used buffer
981 2200 continue
982c---check output buffers.
983c---element neo?
984      inb=imaa*imbf/npoints
985      if(inb.gt.maxbuf)inb=maxbuf
986      do 8014 i=1,inb
987      if(nbuf(i).eq.0)go to 8014
988      neo=neb(i)
989c---yes.  empty buffer.
990      call output(neo,nbuf(i),outcor(1,(i-1)*npoints+1),npoints)
991 8014 continue
992      do 8015 j=0,nel
993      sum=0.
994      sumi=0.
995      do 8016 i=ngood+1,npoints
996      if(cord(7,i).eq.j)then
997         sum=sum+1
998         sumi=sumi+(gam(i)-1.)*erest
999      endif
1000 8016 continue
1001      if(sum.ne.0.)then
1002         sumi=sumi/sum
1003         write(nnout,9021) j,sum/npoints,sumi
1004 9021 format(' fraction lost in element',i3,' = ',f7.5,
1005     *' average energy =',f9.5,' MeV')
1006      endif
1007 8015 continue
1008      if(ngood.gt.2)then
1009         sumi=0.
1010         do 8017 i=1,ngood
1011         sumi=sumi+(gam(i)-1.)*erest
1012 8017    continue
1013         sumi=sumi/ngood
1014         write(nnout,9022) float(ngood)/npoints,sumi
1015 9022    format(' fraction good at end of calculation = ',f7.5,
1016     *   ' average energy = ',f9.5,' MeV ')
1017      endif
1018c-----
1019cbm 18/07/2k      if(lolun1.ne.0) call closehisto
1020c-----
1021      write(ndiag,*) '   '
1022      write(ndiag,*) '  Statistique sur ce run'
1023      write(ndiag,*) '   '
1024
1025      call datime(Iend_date,Iend_time)
1026     
1027c-with cernlib
1028c      call timed(t1)
1029c      CPU_time=t1
1030c_vms/unix
1031c_vms
1032c      write(ndiag,698) 
1033c     * Istart_date(1:4),Istart_date(5:6),Istart_date(7:8),Istart_time,
1034c     * Iend_date(1:4),Iend_date(5:6),Iend_date(7:8),Iend_time,
1035c     * secnds(time_begin),CPU_time
1036c 698  format(//,' --  on ALPHA VMS -- LAL --',28('-')/
1037c     * '   Start Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
1038c     * '   End   Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
1039c     * '   Time used for calculation = ',f10.2,' [sec]',/
1040c     * '   CPU Time                  = ',f10.2,' [sec]',/
1041c     * ' -------------------------------',20('-'))
1042c_unix
1043      if (idebug.eq.0) then
1044      else
1045        write(ndiag,*) '          ---- Version de DEBUG ----- '
1046      endif
1047      write(ndiag,698) kvers,
1048     *   Istart_date,Istart_time,
1049     *   Iend_date, Iend_time
1050c     *   ,secnds(time_begin),CPU_time
1051 698  format(//,' -- LAL -- ',a,16('-')/
1052     *          '   Start Execution : ',a9,' ',a8,/
1053     *          '   End   Execution : ',a9,' ',a8,/
1054     *          ' -------------------------------',20('-'))
1055c les 2 lignes supprimees apres suppression de la cernlib
1056c retrait des appel a la fonction secnds et a la subroutine timed
1057c     *          '   Time used for calculation = ',f10.2,' [sec]',/
1058c     *          '   CPU Time                  = ',f10.2,' [sec]',/
1059
1060c_vms/unix
1061c------------------------------------------------------------------------
1062      call appendparm
1063      if(idebug.eq.0) then
1064        close(unit=nnout,status='delete')
1065      else
1066        close(unit=nnout,status='keep')
1067        write(*,*) '  '
1068        write(*,*) ' CLOSE statement keep file parmdebug '
1069        write(*,*) '  '
1070      endif
1071      print *, ' parmela normal end'
1072      stop ' normal end'
1073c---error messages
1074 7000 continue
1075      write(nnout,9023)
1076 9023 format(' too many parameters on element card')
1077      jj=100
1078      go to 99
1079 7001 continue
1080      write(nnout,9024) lmx
1081 9024 format(' too many elements max is ',i4)
1082      go to 2200
1083 7002 write(nnout,9025)
1084 9025 format(' not enough parameters')
1085      jj=100
1086      go to 99
1087      end
1088
1089c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
1090      subroutine cellfld(cl,nc)
1091c     cell initialization
1092c----------------------------------------------------------------------
1093      save
1094c
1095      include 'param_sz.h'
1096      include 'cfldscom.h'
1097      include 'constcom.h'
1098      include 'misccom.h'
1099      include 'syscom.h'
1100      include 'ucom.h'
1101c
1102      common/jones/ajones,zjones,zcath(imaa)
1103      common/fcoef/ncoef(100),noff(100),ncer(100)
1104c--------------------------------------------------------------------------
1105c*
1106      hcll(nc)=.5*cl
1107      beta=2.*cl/wavel
1108      cellfreq(nc)=freq
1109      lctype=ltype
1110      if(nn.ge.9)then
1111      if(vv(9).ne.0.)cellfreq(nc)=vv(9)
1112      beta=beta*cellfreq(nc)/freq
1113      endif
1114      if(nn.ge.10.and.vv(10).ne.0.)lctype=vv(10)
1115      phi=el(4,nel)*radian
1116      el(9,nel)=sin(phi)
1117      el(10,nel)=cos(phi)
1118      el(4,nel)=cellfreq(nc)
1119      do 30 i = 1,numcf
1120      fc(i,nc) = 0.
1121   30 continue
1122      if (ip.eq.1 .or. ip.eq.2) then
1123        write(nnout,20) lctype,cellfreq(nc),hcll(nc)
1124   20   format (
1125     . ' cell type =',t40,i8/
1126     . ' cell frequency =',t40,1pg12.3,' MHz'/
1127     . ' "half-cell" length =',t40,g12.3,' cm')
1128        if (nel.eq.1.and.ajones.gt.0.) then
1129          write(nnout,21)
1130   21     format (' rf fields from analytic form of M. Jones')
1131          return
1132        endif
1133      endif
1134      if(nn.gt.11) then
1135      np=nn
1136      ncoef(nel)=vv(12)
1137      noff(nel)=vv(13)
1138      ncer(nel)=vv(14)
1139      nmax=14+ncoef(nel)
1140      if(np.gt.nmax)np=nmax
1141      do 10 i=15,nmax
1142      fc(i-14,nc)=vv(i)
1143   10 continue
1144        if (ip.eq.1 .or. ip.eq.2) then
1145          write(nnout,39) ncoef(nel),noff(nel)
1146   39     format (' No. of Fourier coefs =',t40,i8/
1147     .    ' odd/even cosines =',t40,i8)
1148          write(nnout,40) nmax
1149   40     format (' Using vv(15 to ',i2,') as Fourier coefficents of'
1150     .    ,' this cell')
1151          write(nnout,41) ncer(nel)
1152   41     format (' No. of Fourier coefs used in er calculation =',
1153     .    t50,i8)
1154        endif
1155        go to 100
1156      endif
1157c        get fourier coefficients for special cells from tables
1158      ncoef(nel) = 14
1159      noff(nel)  =  1
1160      ncer(nel)  = 14
1161      iflag(nc)=1
1162      if(iflag(nc).ne.0)return
1163      if(lctype.eq.8.and.abs(hcll(nc)-2.21).lt..01)beta=1.02
1164      if(lctype.eq.6)call rtmfc(beta,fc(1,nc))
1165      if(lctype.eq.2)call sccfc(beta,fc(1,nc))
1166      if(lctype.ne.1.and.lctype.ne.2)call genfc(beta,fc(1,nc))
1167      return
1168c         make lookup table for rf field
1169  100 continue
1170      call genfld(nc)
1171      return
1172      end
1173c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
Note: See TracBrowser for help on using the repository browser.