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

Last change on this file since 99 was 80, checked in by lemeur, 12 years ago

prise en compte d'un path pour les fichiers

File size: 45.5 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-----------------------------------------------------------------------
235
236c modif F. Touze, G. Le Meur 15/11/2012
237c---->program parmela
238      subroutine parmela(work_path, lonnom)
239c-----------------------------------------------------------------------
240c label n --> n*100
241c go to  7000 --> 7008
242c do     8000 --> 8017
243c format 9000 --> 9025
244c-----------------------------------------------------------------------
245c
246      include 'param_sz.h'
247      include 'var_char.h'
248      include 'coordcom.h'
249      include 'constcom.h'
250      include 'debug.h'
251      include 'bfieldcom.h'
252      include 'flagcom.h'
253      include 'misccom.h'
254      include 'ncordscom.h'
255      include 'outcom.h'
256      include 'phizoutcom.h'
257      include 'pipecom.h'
258      include 'sccom.h'
259      include 'syscom.h'
260      include 'tstepcom.h'
261      include 'ucom.h'
262      include 'upawcom.h'
263
264      character*132 work_path
265      character*256 nom_input
266      character*256 nomfic
267      character*256 fileelec
268      character*256 ficpoiss
269      integer lonnom
270
271      character*25 kvers
272      common/aswap/acor(7,imaa)
273      common/back/bmax,byz0,byzs,pzowmin,pzowmax,iback
274      common/cathcur/rcur
275      common/chars/title
276      common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7
277      common/ct/var,rp,dzmin,dzmax,wrm,dwbm,xrmsbg,yrmsbg,
278     1xsmax,xpsmax,ysmax,ypsmax,istat,nebm
279      common/ddat/ctime,ddat
280      common/errors/eraln(5,0:lmx)
281      common/image/jj,ij
282      common/intype/intype ! mcd
283      common/jones/ajones,zjones,zcath(imaa) ! mcd
284      common/out6/nbphase,nbremesh
285      common/outbuf/outcor(8,imaa*imb)
286      common/outpar/noutp,outpar(0:lmx,10),moutp(10)
287      common/part/nparticle(imaa)
288      common/pnum/ipnum(imaa)
289      common/save/zrsav
290      common/sbload1/tsbparm(5,50)
291      common/tapes/ifl(2)
292      common/version/kvers
293      common /wtstep/wtstep,nstep ! mcd
294c--------------------------------------------------------------------------
295
296      character Istart_date*9,Iend_date*9,Istart_time*8,Iend_time*8
297      real*4    time_begin
298
299c_with cernlib
300c      call timest(1.0e9)
301c      call timed(t0)
302c
303
304      write(*,*) ' Go Parmela'
305
306
307c_for ibm
308c      time_begin = secnds(0.0)
309
310      call datime(Istart_date,Istart_time)
311c*
312c---set constants
313c
314       kvers='V503 on LINUX'
315c       kvers='V503 on OSF'
316c
317      pi=4.*atan(1.)
318      twopi=2.*pi
319      radian=pi/180.
320      clight=29979.29   ! 10.e6 cm/s
321      erest=.511        ! MeV
322      brhof=1704.51     ! gauss.cm moc/e
323c     brhof = (9.109389E-31*299792458/1.602177E-19)*(100)*(10000) !1704.5
324      zlimit=1000000.   ! cm
325c---set flags
326      idebug=0          ! for debuging in subroutine ouvrir
327      ifld=0            ! flag for external electric field
328      poiflag=.false.   ! flag for external magnetic field
329      nel=0             ! number of element
330      ip=0              ! flag of output type
331      isc=0             ! flag for space charge calculation
332      istat=0           ! flag for statistics
333      ifl(1)=0
334      ifl(2)=0
335      iback=0
336      ltype=1           ! by default disk and washer cells linac
337      ifoclal=0         ! bm bzlal
338      npipe=0           ! mcd
339      nstep=0           ! mcd
340      wtstep=0.         ! mcd
341      ajones=0.         ! mcd
342      nout=0            ! mcd
343      do 8001 n=1,10    ! mcd
344      moutp(n)=0        ! mcd
345      do 8001 i=1,lmx   ! mcd
346      outpar(i,n)=0.    ! mcd
347 8001 continue          ! mcd
348      icfield = 0 ! cbm
349      ictrw= 0    ! cbm
350      nbphase=0   ! cbm nombre de phase en sortie
351      nbremesh=0  ! cbm nombre de remesh
352c----set polarization to zero
353      do 7999 n=1,imaa
354         ksi1(n)=0.
355         ksi2(n)=0.
356         ksi3(n)=0.
357 7999 continue
358c---
359c---
360      nom_input = nomfic(work_path(1:lonnom),'parmin')
361      print *, ' parmela : fichier input : ', nom_input
362      open(unit=nin,file=nom_input,status='old',
363     *     access='sequential',form='formatted')
364
365
366c$$$      open(unit=nin,file='parmin',status='old',
367c$$$     *     access='sequential',form='formatted')
368      open(unit=nnout,file='parmdebug',status='unknown',
369     *        access='sequential',form='formatted')
370      call ouvrir(lonnom, work_path)
371c---
372      call setwall
373c---set misalignments to zero.
374      do 8002 n=1,100
375      do 8002 i=1,5
376      eraln(i,n)=0.
377 8002 continue
378c---read input data
379c---the read(jj)  reads one card with a label and all
380c---subsequent cards without labels.  the stop subroutine compares the label
381c---with a list of words and assigns a value to jj corresponding to the
382c---position of the word in the dictionary.  the word list and corresponding
383c---values of jj are:
384c---(drift     ,1), (solenoid  ,2), (quad      ,3), (bend      ,4),
385c---(buncher   ,5), (chopper   ,6), (cell      ,7), (tank      ,8),
386c---(trwave    ,9), (coil      ,10),(run       ,11),(input     ,12),
387c---(output    ,13),(title     ,14),(scheff    ,15),(zout      ,16),
388c---(adjust    ,17),(start     ,18),(restart   ,19),(continue  ,20),
389c---(save      ,21),(end       ,22),(zlimit    ,23),(errors    ,24),
390c---(change    ,25),(rotate    ,26),(sbload    ,27),(cfield    ,28),
391c---(dpout     ,29),(cathode   ,30),(design    ,31),(pipe      ,32),
392c---(foclal    ,33),(backb     ,34),(wiggler   ,35),(alpham    ,36),
393c---(stat      ,37),(poisson   ,38),(sextupole ,39)
394c---the subroutine reads numbers from the card in a free format.
395c---the values are stored in vv, and there are nn of them.
396 99   continue
397cbm      call cread(jj)
398      call cread
399      go to (100,100,100,100,100,100,100,100,100,1000,
400     1       1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,
401     2       2100,2200,2300,2400,2500,100,100,2800,100,100,
402     3       3100,3200,3300,3400,100,100,3700,3800,100),jj
403c---after transport element or linac cell
404c---after drift:     vv(n)=l,aper,iout
405c---after solenoid:  vv(n)=l,aper,iout,h
406c---after quad:      vv(n)=l,aper,iout,hp
407c---after chopper:   vv(n)=0,aper,iout,freq,phi,dphi
408c---after cathode:   vv(n)=0,aper,iout,radius of curvate of cathode
409c---after wiggler:   vv(n)=l,aper,iout,periods,steps,b0,ky/kw
410c---after sextupole: vv(n)=l,aper,iout,Bpole/ape**2
411 100  continue
412      if(nn.gt.11.and.jj.ne.7.and.jj.ne.9.and.jj.ne.27)go to 7000
413      if(nel.ge.lmx)go to 7001
414      if(nn.lt.2)go to 7002
415      nel=nel+1
416      ntype(nel)=jj
417      do 8003 n=1,nn
418      el(n,nel)=0.      ! mcd
419      if(n.gt.nn)go to 8003
420      el(n,nel)=vv(n)
421 8003 continue
422      zloc(nel)=el(1,nel)
423      if(nel.gt.1)zloc(nel)=zloc(nel)+zloc(nel-1)
424      if (ip.eq.1 .or. ip.eq.2) then
425        iout = vv(3)
426        write(nnout,9001) nel,jj,vv(1),zloc(nel),vv(2),iout
427 9001   format (
428     . ' element no. [ne] =',t40,i8/
429     . ' element type [ntype(ne)] =',t40,i8/
430     . ' element length [el(1,ne)] =',t40,1pg12.4,' cm'/
431     . ' z coord at downstream end [zloc(ne)] =',t40,g12.4,' cm'/
432     . ' radial aperture at exit [el(2,ne)] =',t40,g12.4,' cm'/
433     . ' output at end of element [el(3,ne)] =',t40,i8)
434        if((jj.eq.3) .or. (jj.eq.39)) write(nnout,9002) vv(4)
435 9002   format(
436     . ' magnetic field gradient = ',t40,g12.3,' gauss/cm')
437      endif
438      if(jj.eq.2) ifld=1
439      if(jj.eq.26) go to 2600
440      if(jj.eq.29) go to 99
441      if(jj.eq.27) go to 2700
442      if(jj.eq.30) rcur=vv(4)
443      if(jj.eq.35) go to 3500
444      if(jj.eq.36) go to 3600
445      if(jj.eq.39) go to 99
446      go to (99,99,99,400,500,99,700,800,900), jj
447c---after bend:  vv(n)=l,aper,iout,wr,alpha,beta1,beta2,psi1,psi2,r1,r2
448  400 continue
449      do 8004 n=5,9
450      el(n,nel)=el(n,nel)*radian
451 8004 continue
452      go to 99
453c---after buncher:  vv(n)=0.,aper,iout,wmax,freq,phi,w
454  500 continue
455      w=w0
456      if(nn.ge.7)w=el(7,nel)
457      g=w/erest
458      bg=sqrt(g*(2.+g))
459      cay=twopi*el(5,nel)/(bg*clight)
460      el(8,nel)=cay
461      go to 99
462c---after trwave;
463  900 continue
464      ictrw = ictrw+1
465      call trwfield
466      go to 99
467c---after cell:  vv(n)=cl,aper,iout,phi,e0,nc,(dwtmax)
468  700 continue
469      el(5,nel)=.01*el(5,nel)
470      if (nn.lt.8) el(8,nel)=0.
471      phideg(nel)=el(4,nel)
472c      phi=el(4,nel)*radian this is now defined in cellfld
473c      el(9,nel)=sin(phi)    and this
474c      el(10,nel)=cos(phi)  and this also
475c---if unspecified, set maximum step size in cell to 10 degrees.
476      if(nn.lt.7)el(7,nel)=10.
477c---get field information for this cell
478      cl=el(1,nel)
479      if (el(8,nel).ne.0.) cl=2.*cl
480      nc=el(6,nel)
481c      iel(6,nel)=nc
482      if (ip.eq.1 .or. ip.eq.2) then
483        write(nnout,9003)vv(4),vv(5),nc,el(7,nel),el(11,nel),el(8,nel)
484 9003   format (
485     . ' phase of cell [phi0] =',t40,g12.3,' degrees'/
486     . ' average axial electric field [el(5,ne)] =',t40,g12.3,' MV/m'/
487     . ' cell no. [nc] =',t40,i8/
488     . ' max step size [dwtmax] =',t40,g12.3,' degrees'/
489     . ' axial static magnetic field =',t40,g12.3,' gauss'/
490     . ' cell symmetry [el(8,ne)] =',t40,g12.3)
491        if (el(8,nel).eq.0.) write(nnout,9004)
492 9004   format (' cell is symmetric about midplane in z')
493        if (el(8,nel).gt.0.) write(nnout,9005)
494 9005   format (' half cell, with zero field at entrance')
495        if (el(8,nel).lt.0.) write(nnout,9006)
496 9006   format (' half cell, with zero field at exit')
497      endif
498      call cellfld(cl,nc)
499      go to 99
500c---after tank:  tl,aper,iout,e0t,nc,p,phi
501  800 continue
502      el(8,nel)=el(1,nel)/el(5,nel)
503c---convert e0t from mv/m to mv/cm
504      el(4,nel)=.01*el(4,nel)
505c---calculate beta and k**2 for the tank
506      beta=2.*el(8,nel)/(el(6,nel)*wavel)
507      cl=el(8,nel)
508      caysq=((wavel/(2.*cl))**2 - 1.)*(twopi/wavel)**2
509      el(9,nel)=beta
510      el(10,nel)=caysq
511      go to 99
512c---after field:  vv parameters depend on field subroutine
513 1000 continue
514      icfield=icfield+1
515      call field
516c---define field structure after coil
517c--- vv(n)=z,r,ni,(zbeg,zend,nptcbu,order)
518c---in ( ) required on first coil card ignored on rest
519      go to 99
520c---after errors.  generate errors in geners
521 2400 continue
522      call geners
523      go to 99
524c---after run:  vv(n)=irun,ip,freq,z0,w0,ltype
525 1100 continue
526c     process run card         ! mcd
527      if (nn.ge.1) irun=vv(1)  ! mcd
528c--debug
529c-- Si ip negatif sur carte run on a la version de debug
530      if (nn.ge.2) then
531         if(vv(2).lt.0.) then
532            idebug=1
533            write(*,*) '-----------------------------------'
534            write(*,*) 'You have ask the version with debug'
535            write(*,*) '-----------------------------------'
536            write(nnout,*) '-----------------------------------'
537            write(nnout,*) 'You have ask the version with debug'
538            write(nnout,*) '-----------------------------------'
539         else
540         endif
541         ip=abs(vv(2))
542      endif
543c--debug
544      if (nn.ge.3) freq=vv(3)  ! mcd
545      if (nn.lt.4) go to 99 
546      if (nn.ge.4) z0=vv(4)    ! mcd
547      if (nn.lt.5) go to 99
548      if (nn.ge.5) w0=vv(5)    ! mcd
549      if (nn.ge.6) ltype=vv(6) ! mcd
550      if(freq.gt.0.) wavel=clight/freq
551      write(nnout,9007)irun
552 9007 format (/' '/
553     . ' run no. =',t40,i8/)
554      if((ip.gt.4).and.(ip.ne.999)) then
555      ip=1
556      write(nnout,9008)ip
557 9008 format(
558     .' print option modified =',t40,i8)
559      else
560      write(nnout,9009)ip
561 9009 format(
562     . ' print option =',t40,i8)
563      endif
564c---
565      if (ip.eq.3) then
566      open(unit=nuch,file='parmtrch',access='sequential',
567     *     status='unknown')
568      write(nuch,*)'  nc  sym  r  z  ez  er  bt  '
569      endif
570      if (ip.eq.999) open(unit=nimp,file='parmimp',status='unknown')
571c---
572      write(nnout,9010)freq,wavel,z0,w0,ltype
573 9010 format(
574     . ' rf frequency =',t40,f9.2,'    MHz'/
575     . ' wavelength =',t40,1pg12.3,' cm'/
576     , ' z0 =',t40,g12.3,' cm'/
577     , ' reference particle energy =',t40,g12.3,' MeV'/
578     . ' linac type =',t40,i8)
579      do 8005 i=1,4
580      corout7(i,1)=0.
581      cor(i,1)=0.
582 8005 continue
583      cor(5,1)=z0
584      corout7(5,1)=z0
585      corout7(6,1)=w0
586      g=w0/erest
587      gam(1)=g+1
588      cor(6,1)=sqrt(g*(2.+g))
589      cord(6,1)=cor(6,1)
590      if (npoints.lt.1)npoints=1
591      ngood=npoints
592      nparticle(1)=1
593      phizero(1)=0.
594      if(nn.lt.6)go to 99
595      go to 99
596c---after input:  vv parameters depend on input subroutine
597 1200 continue
598      call input
599      ngood=npoints
600      do 8006 i=1,ngood
601      nparticle(i)=i
602      ipnum(i)=i
603 8006 continue
604      go to 99
605c---after output:  vv parameters stored in optcon array
606c---after output:  vv parameters stored in outpar array
607c
608 1300 continue
609      noutp=vv(1)+1
610      if(noutp.gt.10 .or. noutp.lt.1) then
611      write(ndiag,9011)
612 9011 format(' invalid OUTPUT type')
613      go to 99
614      endif
615      moutp(noutp)=nn
616      do 8007 n=1,nn
617      optcon(n)=vv(n)
618      outpar(n,noutp)=vv(n)
619 8007  continue
620      if(vv(1).gt.1..and.vv(1).lt.5.)then
621      imbf=imb-1
622      else
623      imbf=imb
624      endif
625      go to 99
626c---after title:  read title information from next card
627 1400 continue
628      read(nin,9012)title
629 9012 format(a80)
630      write(nnout,*)title
631      go to 99
632c---after scheff:  vv(1)=beam current; vv(2)= program number
633c--- rest depend on space-charge subroutines
634 1500 continue
635      do 8008 n=1,12
636      sce(n)=0.
637      if (n.le.nn) sce(n)=vv(n)
638 8008 continue
639      beami=vv(1)
640      iprog=vv(2)
641      if(isc.eq.0) then
642        if(iprog.eq.1) call scheff1(0.,-99,0.)
643        if(iprog.eq.2) call scheff2(0.,-99,0.)
644        if(iprog.eq.3) call scheff3(0.,-99,0.)
645         isc=1
646      else          ! isc
647      endif         ! isc
648      go to 99
649c---after zout.  print locations of elements
650 1600 continue
651      if(ifoclal.eq.1) then
652        call zoutlal
653      else
654        call zout
655      endif
656      go to 99
657c---after adjust:  call arbitrary subroutine named adjust
658 1700 continue
659      call adjust
660      go to 99
661c---after change.  call change subroutine
662 2500 continue
663      call change
664      go to 99
665c---after rotate
666 2600 continue
667      el(4,nel)=el(4,nel)*radian
668      go to 99
669c---  after sbload paramaters are l,aperature,iout,zlen,ncells,a1,a2,a3
670c---  ,a4,a5,ta1,ta2,ta3,ta4,ta5. if ncells is negative then the paramaters
671c---  ta1 to ta4 are required. ta5 is optional. if ncells is negative the
672c---  transvers wake will be applied also.
673 2700 continue
674      if(nn.gt.11)then
675      numsb=numsb+1
676      tsbparm(1,numsb)=vv(11)
677      el(11,nel)=numsb
678      tsbparm(2,numsb)=vv(12)
679      tsbparm(3,numsb)=vv(13)
680      tsbparm(4,numsb)=vv(14)
681      if(el(5,nel).lt.0..and.nn.lt.14)go to 7002
682      if(nn.gt.14)tsbparm(5,numsb)=vv(15)
683      endif
684      go to 99
685c---after start: start dynamics calculations
686 1800 continue
687      if(nn.lt.5)go to 7002
688      wt0=vv(1)
689c---  For input type 10 and 11, WTSTEP is the half width of the laser pulse.
690c     The clock must be set back by this amount to generate the time
691c     distribution of photoelectrons.
692      wt0=wt0-wtstep
693      dwt=vv(2)
694      nsteps=vv(3)
695      nsc=vv(4)
696      nout=vv(5)
697      go to 7003
698c---after restart.  set the starting time (phase), initialize cord
699c---from cor, and find largest element number (len) that any particle is in.
700 1900 if(nn.lt.4)go to 7002
701      dwt=vv(1)
702      nsteps=vv(2)
703      nsc=vv(3)
704      nout=vv(4)
705c       -----------------
706      if(jj.eq.19) then                         
707                npoints=-vv(6)
708      endif
709c       -------------------------
710      if(nn.gt.5.and.(vv(6).gt.0..and.vv(6).le.26))then
711        i=vv(6)-1
712      call orest1(i)
713      else
714      call orest2
715      endif
716      read(nsav) cor,wt0,ngs
717      close(nsav)
718      ngood=ngs
719C       ----------------------------
720      IF(NGS.NE.NPOINTS) then
721        call appendparm
722        print *,'ngs=',ngs,'npoints=',npoints
723        stop ' Abnormal stop ngs.ne.npoints main '
724      endif
725C       -------------------------------
726 7003 continue
727      if (ip.eq.1 .or. ip.eq.2)
728     .   write (nnout,9013) wt0,dwt,nsteps,vv(5),nsc,nout
729 9013 format (
730     . ' starting phase =',t40,1pg12.3,' degrees'/
731     . ' integration step size =',t40,g12.3,' degrees'/
732     . ' no. of steps in integration =',t40,i8/
733     . ' shift of reference particle =',t40,g12.3,' cm'/
734     . ' call SCHEFF every',t40,i8,'     step(s)'/
735     . ' call OUTPUT every',t40,i8,'     step(s)'/)
736      wt=wt0
737      do 8009 i=1,maxbuf
738      neb(i)=0
739      nbuf(i)=0
740 8009 continue
741c---if nn.gt.4, the 5th parameter is dist to shift ref particle
742c   only for restart. not for start
743      if (jj.eq.19.and.nn.gt.4) then
744        cor(5,1)=cor(5,1)+vv(5)
745      endif
746c       --------------------------------------------------
747c--------- vv(5)= dz
748      if(jj.eq.19.and.vv(5).ne.0.) then
749            cor(1,1)=0.  ! for energy change of ref part
750            cor(2,1)=0.
751            cor(3,1)=0.
752            cor(4,1)=0.
753            print *,'vv(7) =',vv(7)
754            cor(6,1)=sqrt((1.+vv(7)/0.511)**2-1.)
755      endif
756c       ---------------------------------------------------
757      do 8010 n=1,ngood
758      do 8011 i=1,6
759      cord(i,n)=cor(i,n)
760 8011  continue
761      gam(n)=sqrt(1.+cord(2,n)**2+cord(4,n)**2+cord(6,n)**2)
762c---find which element each particle is in
763      ne=0
764      z=cord(5,n)
765      if(z.lt.0.)go to 7004
766      do 8012 i=1,nel
767      ne=i
768      if(z.lt.zloc(i))go to 7004
769 8012 continue
770      ne=nel+1
771 7004 cord(7,n)=ne
772c     if generating from a photocathode, put ne=0
773      if(jj.eq.19) go to 8010           
774      if (intype.eq.10 .or. intype.eq.11 .or. intype.eq.15) 
775     .   cord(7,n)=0.
776 8010 continue
777c---if nn.gt.4, the 5th parameter is dist to shift ref particle
778c   only for restart. not for start
779c---start dynamics calculation
780      call pardyn
781      go to 99
782c---after continue:  used to change step size and continue computations
783 2000 continue
784      if(nn.lt.4)go to 7002
785      dwt=vv(1)
786      nsteps=vv(2)
787      nsc=vv(3)
788      nout=vv(4)
789c---if nn.gt.4, the 5th parameter is dist to shift ref particle
790      if (nn.gt.4)cord(5,1)=cord(5,1)+vv(5)
791      call pardyn
792      go to 99
793c---after save
794c---save cord in cor
795 2100 continue
796      ngs=ngood
797      do 8013 n=1,ngood
798      do 8013 i=1,6
799      cor(i,n)=cord(i,n)
800 8013 continue
801      write(nnout,*) 'zr of save ',zrsav
802      wt0=wt
803      if(nn.ge.1.and.(vv(1).gt.0..and.vv(1).le.26))then
804      i=vv(1)-1
805      call osave1(i)
806      if(ifl(2).ne.0) then
807      close(nsav)
808      ifl(2)=0
809      endif
810      if(idebug.eq.0) then
811        close(unit=nnout,status='delete')
812      else
813        close(unit=nnout,status='keep')
814      endif
815      close(nnout)
816      open(unit=nnout,file='parmout',access='sequential',
817     *     status='unknown',form='formatted')
818 7005 continue
819      read(nnout,'(10a8)',err=7006)tem
820      go to 7005
821 7006 continue
822      backspace (nnout)
823      else
824      open (unit=nsav,file='savecor',access='sequential',
825     *      form='unformatted',status='new')
826      endif
827      write(nsav) cor,wt0,ngs,bgs
828      close(nsav)
829      go to 99
830c---after zlimit:  lower limit on lagging particles
831 2300 continue
832      zlimit=vv(1)
833      go to 99
834c---after cfield:  read fields from sfo7 as written on tape7
835c--- first parameter on cfield is starting cell number
836c--- next parameter is ending cell number and third parameter is
837c--- increment of cell number as in a do loop
838c---line after cfield is filename of sfo7 data
839 2800 continue
840C modif glm 29/10/2012      read(nin,'(a12)')filesuper
841      read(nin,*) filesuper
842      print *, ' j''ai lu le fichier champ : ', filesuper
843      write(nnout,'(1x,a12)')filesuper
844      if(nn.lt.1)then
845      nc=1
846      ne=1
847      ns=1
848      elseif(nn.eq.1)then
849      nc=vv(1)
850      ne=nc
851      ns=1
852      elseif(nn.eq.2)then
853      nc=vv(1)
854      ne=vv(2)
855      ns=1
856      elseif(nn.ge.3)then
857      nc=vv(1)
858      ne=vv(2)
859      ns=vv(3)
860      endif
861
862c$$$      call cellfeld(filesuper,nc,ne,ns)
863      fileelec = nomfic(work_path(1:lonnom),filesuper)
864      print *, ' j''envoie lu le fichier champ : ', fileelec
865      call cellfeld(fileelec,nc,ne,ns)   ! modif glm 16/11/2012
866      go to 99
867c---after design
868 3100 continue
869      call design
870      go to 99
871c---after pipe, vv(n) = radius,zlow,zhigh
872 3200 continue
873      if (npipe.eq.100) then
874        write(nnout,9014)
875 9014   format (' Too many PIPE cards')
876        go to 99
877      endif
878      if (vv(1).le.0.) then
879        write(nnout,9015)
880 9015   format (' Bad radius on PIPE card')
881        go to 99
882      endif
883      npipe = npipe + 1
884      rpipe(npipe) = vv(1)
885      zlpipe(npipe) = vv(2)
886      zhpipe(npipe) = vv(3)
887      if (ip.eq.1 .or. ip.eq.2)
888     .write(nnout,9016) npipe,vv(1),vv(2),vv(3)
889 9016 format (
890     . ' Beam pipe for point-by-point image calculation'/
891     . ' pipe no. = ',t40,i8/
892     . ' pipe radius = ',t40,1pg12.3,' cm'/
893     . ' zlow = ',t40,g12.3,' cm'/
894     . ' zhigh = ',t40,g12.3,' cm')
895      go to 99
896c---after foclal vv(n)=zmin,zmax,dzz,npchb,facbz,opt
897 3300 continue
898      iopt=vv(6)
899      if(iopt.eq.0) then
900         ifld=0
901         ifoclal=0
902         read(nin,'(a12)') filebz
903         write(nnout,*) 'foclal not used'
904      else
905         ifld=1
906         ifoclal=1
907         read(nin,'(a12)') filebz
908         call fieldlal(filebz)
909      endif
910      go to 99
911c---after backb back bombardment
912 3400 continue
913      iback=1
914      bmax=vv(1)
915      byz0=vv(2)
916      byzs=vv(3)
917      iope=0
918      if(nn.gt.3) then
919           pzowmin=vv(4)
920           pzowmax=vv(5)
921      else
922           pzowmin=0.
923           pzowmax=0.
924      endif
925      if(pzowmax.ne.0.) iope=1
926      call opeback(iope)
927      go to 99
928c---after wiggler l,aper,iout,periods,steps,b0,kx,ky,kw,iprint
929 3500 continue
930      if(nn.lt.7) go to 7002
931      xlam=el(1,nel)/el(4,nel)
932      cayw=twopi/xlam
933      cayy=el(7,nel)*cayw
934      if(el(7,nel).le.0. .or. el(7,nel).gt.1.)then
935         write(nnout,9018)
936 9018    format(' error: ky/kw on wiggler must be between 0 and 1')
937         jj=100 ! for cread
938         go to 99
939      endif
940      cayx=0.
941      cayxsq=cayw**2-cayy**2
942      if(cayxsq.gt.0.) cayx=sqrt(cayxsq)
943      el(7,nel)=cayx
944      el(8,nel)=cayy
945      el(9,nel)=cayw
946      if(ip.eq.1 .or. ip.eq.2)
947     .write(nnout,9019)el(4,nel),el(5,nel),el(6,nel),el(7,nel)
948 9019 format(
949     .' number of wiggler periods =',t40,g12.3/
950     .' step of wiggler =',t40,g12.3/
951     .' wiggler field =',t40,g12.3,' gauss'/
952     .' ratio of Ky to Kw =',t40,g12.3)
953      go to 99
954c---after alpha magnet
955 3600 continue
956      if(ip.eq.1 .or. ip.eq.2)
957     .write(nnout,9020) el(1,nel),el(4,nel),el(5,nel),el(7,nel)
958 9020 format(
959     .' distance travelled by the r.p. =',t40,g12.3,' cm'/
960     .' reference particle energy =',t40,g12.3,' MeV'/
961     .' alpha magnet angle =',t40,g12.3,' Degrees'/
962     .' absolute field gradient =',t40,g12.3,' Gauss/cm')
963      go to 99
964c---after stat
965 3700 continue
966      var=vv(1)
967      istat=1
968      open(unit=nstat,file='parmstat',form='unformatted',
969     *     status='unknown')
970      write(nstat)ij
971      go to 99
972c---after Poisson
973 3800 continue
974C modif glm 29/10/2012      read (nin,'(a12)')filepois
975      read (nin,*)filepois
976      write(nnout,'(1x,a12)')filepois
977      if(nn.ge.1)then
978      zoffset=vv(1)
979      else
980      zoffset=0.
981      endif
982      if(nn.ge.2)then
983      rmult=vv(2)
984      else
985      rmult=1.
986      endif
987      ficpoiss = nomfic(work_path(1:lonnom),filepois)
988      call poisson(ficpoiss,zoffset,rmult)
989      go to 99
990c
991cbm retrait de nfieldb 04/02/2010 apres modif de coil
992c--- after nfieldb number of points used for the magnetic field
993c    < 10000
994c 4000 continue
995c      nptcbu=vv(1)
996c      write(nnout,*)
997c     *' Nombre de points maximal pour le champ magnetique ',nptcb
998c      write(nnout,*)
999c     *' Number of used values for magnetic field',nptcbu
1000c      if(nptcbu.gt.nptcb) then
1001c        write(nnout,*)
1002c     *  ' Number of points for magnetic field gt ',nptcb
1003c      endif
1004c      go to 99 
1005c---after-end output any used buffer
1006 2200 continue
1007c---check output buffers.
1008c---element neo?
1009      inb=imaa*imbf/npoints
1010      if(inb.gt.maxbuf)inb=maxbuf
1011      do 8014 i=1,inb
1012      if(nbuf(i).eq.0)go to 8014
1013      neo=neb(i)
1014c---yes.  empty buffer.
1015      call output(neo,nbuf(i),outcor(1,(i-1)*npoints+1),npoints)
1016 8014 continue
1017      do 8015 j=0,nel
1018      sum=0.
1019      sumi=0.
1020      do 8016 i=ngood+1,npoints
1021      if(cord(7,i).eq.j)then
1022         sum=sum+1
1023         sumi=sumi+(gam(i)-1.)*erest
1024      endif
1025 8016 continue
1026      if(sum.ne.0.)then
1027         sumi=sumi/sum
1028         write(nnout,9021) j,sum/npoints,sumi
1029 9021 format(' fraction lost in element',i3,' = ',f7.5,
1030     *' average energy =',f9.5,' MeV')
1031      endif
1032 8015 continue
1033      if(ngood.gt.2)then
1034         sumi=0.
1035         do 8017 i=1,ngood
1036         sumi=sumi+(gam(i)-1.)*erest
1037 8017    continue
1038         sumi=sumi/ngood
1039         write(nnout,9022) float(ngood)/npoints,sumi
1040 9022    format(' fraction good at end of calculation = ',f7.5,
1041     *   ' average energy = ',f9.5,' MeV ')
1042      endif
1043c-----
1044cbm 18/07/2k      if(lolun1.ne.0) call closehisto
1045c-----
1046      write(ndiag,*) '   '
1047      write(ndiag,*) '  Statistique sur ce run'
1048      write(ndiag,*) '   '
1049
1050      call datime(Iend_date,Iend_time)
1051     
1052c-with cernlib
1053c      call timed(t1)
1054c      CPU_time=t1
1055c_vms/unix
1056c_vms
1057c      write(ndiag,698) 
1058c     * Istart_date(1:4),Istart_date(5:6),Istart_date(7:8),Istart_time,
1059c     * Iend_date(1:4),Iend_date(5:6),Iend_date(7:8),Iend_time,
1060c     * secnds(time_begin),CPU_time
1061c 698  format(//,' --  on ALPHA VMS -- LAL --',28('-')/
1062c     * '   Start Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
1063c     * '   End   Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
1064c     * '   Time used for calculation = ',f10.2,' [sec]',/
1065c     * '   CPU Time                  = ',f10.2,' [sec]',/
1066c     * ' -------------------------------',20('-'))
1067c_unix
1068      if (idebug.eq.0) then
1069      else
1070        write(ndiag,*) '          ---- Version de DEBUG ----- '
1071      endif
1072      write(ndiag,698) kvers,
1073     *   Istart_date,Istart_time,
1074     *   Iend_date, Iend_time
1075c     *   ,secnds(time_begin),CPU_time
1076 698  format(//,' -- LAL -- ',a,16('-')/
1077     *          '   Start Execution : ',a9,' ',a8,/
1078     *          '   End   Execution : ',a9,' ',a8,/
1079     *          ' -------------------------------',20('-'))
1080c les 2 lignes supprimees apres suppression de la cernlib
1081c retrait des appel a la fonction secnds et a la subroutine timed
1082c     *          '   Time used for calculation = ',f10.2,' [sec]',/
1083c     *          '   CPU Time                  = ',f10.2,' [sec]',/
1084
1085c_vms/unix
1086c------------------------------------------------------------------------
1087      call appendparm
1088      if(idebug.eq.0) then
1089        close(unit=nnout,status='delete')
1090      else
1091        close(unit=nnout,status='keep')
1092        write(*,*) '  '
1093        write(*,*) ' CLOSE statement keep file parmdebug '
1094        write(*,*) '  '
1095      endif
1096      print *, ' parmela normal end'
1097      stop ' normal end'
1098c---error messages
1099 7000 continue
1100      write(nnout,9023)
1101 9023 format(' too many parameters on element card')
1102      jj=100
1103      go to 99
1104 7001 continue
1105      write(nnout,9024) lmx
1106 9024 format(' too many elements max is ',i4)
1107      go to 2200
1108 7002 write(nnout,9025)
1109 9025 format(' not enough parameters')
1110      jj=100
1111      go to 99
1112      end
1113
1114c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
1115      subroutine cellfld(cl,nc)
1116c     cell initialization
1117c----------------------------------------------------------------------
1118      save
1119c
1120      include 'param_sz.h'
1121      include 'cfldscom.h'
1122      include 'constcom.h'
1123      include 'misccom.h'
1124      include 'syscom.h'
1125      include 'ucom.h'
1126c
1127      common/jones/ajones,zjones,zcath(imaa)
1128      common/fcoef/ncoef(100),noff(100),ncer(100)
1129c--------------------------------------------------------------------------
1130c*
1131      hcll(nc)=.5*cl
1132      beta=2.*cl/wavel
1133      cellfreq(nc)=freq
1134      lctype=ltype
1135      if(nn.ge.9)then
1136      if(vv(9).ne.0.)cellfreq(nc)=vv(9)
1137      beta=beta*cellfreq(nc)/freq
1138      endif
1139      if(nn.ge.10.and.vv(10).ne.0.)lctype=vv(10)
1140      phi=el(4,nel)*radian
1141      el(9,nel)=sin(phi)
1142      el(10,nel)=cos(phi)
1143      el(4,nel)=cellfreq(nc)
1144      do 30 i = 1,numcf
1145      fc(i,nc) = 0.
1146   30 continue
1147      if (ip.eq.1 .or. ip.eq.2) then
1148        write(nnout,20) lctype,cellfreq(nc),hcll(nc)
1149   20   format (
1150     . ' cell type =',t40,i8/
1151     . ' cell frequency =',t40,1pg12.3,' MHz'/
1152     . ' "half-cell" length =',t40,g12.3,' cm')
1153        if (nel.eq.1.and.ajones.gt.0.) then
1154          write(nnout,21)
1155   21     format (' rf fields from analytic form of M. Jones')
1156          return
1157        endif
1158      endif
1159      if(nn.gt.11) then
1160      np=nn
1161      ncoef(nel)=vv(12)
1162      noff(nel)=vv(13)
1163      ncer(nel)=vv(14)
1164      nmax=14+ncoef(nel)
1165      if(np.gt.nmax)np=nmax
1166      do 10 i=15,nmax
1167      fc(i-14,nc)=vv(i)
1168   10 continue
1169        if (ip.eq.1 .or. ip.eq.2) then
1170          write(nnout,39) ncoef(nel),noff(nel)
1171   39     format (' No. of Fourier coefs =',t40,i8/
1172     .    ' odd/even cosines =',t40,i8)
1173          write(nnout,40) nmax
1174   40     format (' Using vv(15 to ',i2,') as Fourier coefficents of'
1175     .    ,' this cell')
1176          write(nnout,41) ncer(nel)
1177   41     format (' No. of Fourier coefs used in er calculation =',
1178     .    t50,i8)
1179        endif
1180        go to 100
1181      endif
1182c        get fourier coefficients for special cells from tables
1183      ncoef(nel) = 14
1184      noff(nel)  =  1
1185      ncer(nel)  = 14
1186      iflag(nc)=1
1187      if(iflag(nc).ne.0)return
1188      if(lctype.eq.8.and.abs(hcll(nc)-2.21).lt..01)beta=1.02
1189      if(lctype.eq.6)call rtmfc(beta,fc(1,nc))
1190      if(lctype.eq.2)call sccfc(beta,fc(1,nc))
1191      if(lctype.ne.1.and.lctype.ne.2)call genfc(beta,fc(1,nc))
1192      return
1193c         make lookup table for rf field
1194  100 continue
1195      call genfld(nc)
1196      return
1197      end
1198c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
Note: See TracBrowser for help on using the repository browser.