V3FIT
gtovmi2.f
1 ! 2007 Ed Lazarus
2 ! Replaces older gtovmi based on MAPCODE
3 ! This mapping to equal arc cooordinates based on BALOO
4 ! Revisions:
5 ! December 2007 Correct the calculation of Ienclosed
6 ! Dec. 2007 Added flexibility for reducing number of
7 ! coefficioents in AM, AC, AI
8 ! " pmass_type='Akima_spline'"
9 ! : needs input arrays am_aux_f and am_aux_s
10 ! " pcurr_type='Akima_spline_Ip'"
11 ! : needs input arrays ac_aux_f and ac_aux_s
12 ! " piota_type='Akima_spline'"
13 ! : needs input arrays ai_aux_f and ai_aux_s
14  MODULE knots
15  USE precision
16  REAL(rprec), DIMENSION(:), ALLOCATABLE :: s2, f2
17  INTEGER :: n2
18  END MODULE knots
19  PROGRAM gtovmi2
20  USE precision
21  USE mapg_mod
22  USE gfile
23  USE mapout
24  USE polynomial
25  USE knots
26  IMPLICIT NONE
27  TYPE(GEQDSK) :: g1
28  TYPE(RSZS) :: g2
29  LOGICAL :: lexist=.false.
30  INTEGER :: numarg, iargc, id0, id1, id, np, pgopen, ier=0, m
31  INTEGER :: nlo, nfit, l1, l2, i, ix, ierr=0, ipos, mp=11
32  INTEGER :: kkac=99, kkai=99, kkam=99, count
33 ! LOGICAL :: lexist = .false.
34 ! INTEGER :: numarg, iargc, id0, id1, id, np, pgopen, ier, m
35 ! INTEGER :: nlo, nfit, l1, l2, i, ierr, ipos, mp = 11
36 ! INTEGER :: kkac = 99, kkai = 99, kkam = 99
37  LOGICAL :: exist, afile=.true.
38  CHARACTER*12 USErnm,pname*31,ident*172,formats*99
39  CHARACTER*12 idents1*72,idents2*72,formats1*80,formats2*80
40  CHARACTER(len=22) :: arg2
41  CHARACTER*24 datebuf, fdate, onechr*2, porder*2, cim*6
42  CHARACTER*90 geqdskfile, aeqdskfile
43  real*4, DIMENSION(:),allocatable :: xpl, xpl2, ypl, ypl2
44  REAL(rprec) :: cutoff = 3e-8, cutin, addv, phiedg, rcentr = 1.6955
45  REAL(rprec) :: tol = 1e-3, flux_fraction
46  REAL(rprec), DIMENSION(:),allocatable :: s , a0
47  REAL(rprec), POINTER :: yfit(:),yfit2(:)
48  REAL(rprec), DIMENSION(:), ALLOCATABLE :: xin,yin,b1,b2,b3,b4
49  INTEGER, DIMENSION(:),allocatable :: thinner
50  CHARACTER*90 dul1, formatstring*100
51  CHARACTER*63 vmecinput
52  CHARACTER, DIMENSION (100) :: line*60
53  CHARACTER(len=5) :: device='/cps'
54  CHARACTER*20 xlabel, ylabel, slabel*60, mlabel*70
55 c-----------------------------------------------------------------------
56 c mapperb etc. FROM ballooning code for tokamak equilibrium
57 c miller & linliu 7/96
58 c-----------------------------------------------------------------------
59  percenflux = 0.9990 !psiv(npsi)=psiaxis+(psilim-psiaxis)*percenflux
60  readeqdsk = .true.
61  alpsi = -0
62  alpsi = -1
63  filename = ""
64  numarg = iargc()
65  WRITE(6,'("CALL for g-file")')
66  IF (numarg > 0) THEN
67  CALL getarg(1,filename)
68  INQUIRE(file = trim(filename),exist = lexist)
69  ENDIF
70  IF(numarg == 0 .or. trim(filename).eq.'-h') then
71  WRITE(*,*) "(1) Arguments are g-file name ",
72  . "[, flux_fraction [,graphics-device],[klm]]"
73  WRITE(*,*)"(2) ",
74  . "A zero for flux_fraction will select the default."
75  WRITE(*,*)"(3) Devices uaualy available are: ",
76  . " '/xs','/xw/,'/ps','/cps'"
77  WRITE(*,*)"(4) Polynomial orders in format 3I2.2 ",
78  . "Alphebetical ordering: AC, AI, AM."
79  WRITE(*,*) "NEW:",
80  . " Use any of AC, AI, AM = 00 to invoke Akima splines!"
81  IF (trim(filename).eq.'-h') stop 'help'
82  ENDIF
83  IF(.NOT. lexist) THEN
84  WRITE(*,*)"Enter file name: "
85  READ(*,fmt='(a)')filename
86  WRITE(*,*)"Enter flux fraction: "
87  READ(*,fmt='(f6.14)')flux_fraction
88  IF((flux_fraction.ge.0.9).AND.(flux_fraction.lt.1.0))then
89  percenflux = flux_fraction
90  ELSE
91  print*,'NOT IN RANGE: flux_fraction=',percenflux
92  ENDIF
93  WRITE(*,*)"plot device: "
94  READ(*,fmt='(a)')device
95  INQUIRE(file = trim(filename),exist = lexist)
96  ENDIF
97  INQUIRE(file = trim(filename),exist = lexist)
98  IF(.NOT.lexist)print*,lexist, trim(filename)
99  IF(.NOT.lexist) stop "STOP on nofile"
100  IF (numarg > 1) THEN
101  CALL getarg(2,arg2)
102  READ(arg2,fmt='(f6.14)')flux_fraction
103  IF((flux_fraction.ge.0.9).AND.(flux_fraction.lt.1.0))then
104  percenflux = flux_fraction
105  ELSE
106  print*,'NOT IN RANGE: flux_fraction=',percenflux
107  ENDIF
108  ENDIF
109  if (numarg >= 3) CALL getarg (3, device)
110  if (numarg == 4) then
111  call getarg (4, cim)
112  read(cim,fmt='(3i2.2)')kkac, kkai, kkam
113  write(*,fmt='("kkac, kkai, kkam:",3(x,i2.2))')kkac, kkai, kkam
114  endif
115  geqdskfile = trim(filename)
116  slabel = trim(filename)
117  ipos = index(geqdskfile,'g',.true.)
118  aeqdskfile =
119  & geqdskfile(1:ipos-1)//'a'//geqdskfile(ipos+1:len(geqdskfile))
120  vmecinput = 'input.vmec.'//geqdskfile(ipos+1:)
121  if(kkac == 0 .OR. kkam == 0 .OR. kkai == 0)
122  & vmecinput = 'input.akima.'//geqdskfile(ipos+1:)
123  INQUIRE(file = trim(geqdskfile),exist = exist)
124  IF(.not.exist)stop 'no geqdskfile'
125  INQUIRE(file = trim(aeqdskfile),exist = exist)
126 ! IF(.not.exist)stop 'no aeqdskfile'
127  IF(.not.exist) afile=.false.
128  id0 = pgopen(trim(device))
129  if (id0.le.0)
130  . WRITE(*,*)'Fail to Open Graphics Device',
131  . trim(device),' with code ',id0
132  if (id0.le.0) stop 'PGPLOT device unavailable'
133  CALL pgpap(7.9,1)
134  CALL pgslw(2)
135  CALL pgscf(2)
136  CALL rdeqdsk(g1) !efit equilibrium READ eqdsk
137  WRITE(6,'("begin mapping")')
138  CALL mapperb
139  WRITE(6,'("finished mapping")')
140  CALL dvdpsi
141  CALL mapout_nc(g2)
142  WRITE(6,'("finished writing")')
143 
144 
145 c****************************************************************
146 c ****** initialize vmec input file *******
147  line(01) =
148  .' MGRID_FILE = ''/home/lazarus/mgrids/mgrid_d3d_ef.nc'''
149  line(02) =
150  . ' LFREEB = F LASYM = T'
151  line(03) =
152  . ' DELT = 1.00E-00'
153  line(04) =
154  . ' TCON0 = 2.00E+00'
155  line(05) =
156  . ' NFP = 1'
157  line(06) =
158  . ' NCURR = 01'
159  line(07) =
160  . ' NTOR = 0 NTHETA = 00 NZETA = 1'
161  line(12) =
162  . ' NS_ARRAY = 13 25 49 97'
163  line(11) =
164  . ' NITER_ARRAY = 1500 2500 3500 5000'
165  line(09) =
166  . ' NSTEP = 300'
167  line(10) =
168  . ' NVACSKIP = 3'
169  line(08) =
170  . ' GAMMA = 0.000000E+00'
171  line(13) =
172  . ' FTOL_ARRAY = 1.e-6 1.e-8 1.e-11 1.e-12'
173  WRITE(xlabel,fmt='(es14.8)')g2%fraction_bndry
174  line(14) =
175  . '! mapcode boundary fraction is '//trim(xlabel)
176  OPEN (unit = 67, file = vmecinput , status = 'UNKNOWN')
177  CALL getarg(0,ident)
178  l1 = 0
179  l2 = 0
180  DO i = len(ident),1,-1
181  IF(l2.eq.0.and.ident(i:i).ne.' ')l2 = i
182  IF(l1.eq.0.and.ident(i:i).eq.'/')l1 = i
183  END DO
184  l2 = min(l2,len(pname))
185  pname = ident(l1+1:l2)
186  CALL getlog(usernm)
187  datebuf = fdate()
188  formats = '(6h! file,x,a'
189  WRITE(onechr(1:2),fmt = '(i2.2)')len_trim(vmecinput)
190  formats = formats(1:len_trim(formats))//onechr(1:2)
191  formats1 = formats
192  formats = formats(1:len_trim(formats))//',14h generated by ,'
193  formats1 = formats1(1:len_trim(formats))//',14h generated by )'
194  formats = formats(1:len_trim(formats))//'/,2h! ,'
195  WRITE(onechr(1:2),fmt = '(i2.2)')len_trim(pname)
196  formats2 = '(2h! ,'
197  formats = formats(1:len_trim(formats))//'a'//onechr(1:2)
198  formats2 = formats2(1:len_trim(formats2))//'a'//onechr(1:2)
199  formats = formats(1:len_trim(formats))//',5h for ,'
200  formats2 = formats2(1:len_trim(formats2))//',5h for ,'
201  WRITE(onechr(1:2),fmt = '(i2.2)')len_trim(usernm)
202  formats = formats(1:len_trim(formats))//'a'//onechr(1:2)
203  formats2 = formats2(1:len_trim(formats2))//'a'//onechr(1:2)
204  formats = formats(1:len_trim(formats))//',4h on ,a24,1h!)'
205  formats2 = formats2(1:len_trim(formats2))//',4h on ,a24)'
206  WRITE(idents1,fmt = formats1)vmecinput
207  WRITE(idents2,fmt = formats2)pname,usernm,datebuf
208  line(15) = trim(idents1)
209  line(16) = trim(idents2)
210  WRITE(67,*)'&INDATA'
211  DO i = 1,13
212  WRITE(67, fmt='(a)')trim(line(i))
213  ENDDO
214 c****************************************************************
215  phiedg = g2%phi(g2%npsi)
216  IF (g1%bcentr < 0. ) phiedg = -phiedg
217  WRITE(67, fmt='('' PHIEDGE = '',es14.7)')phiedg
218  IF(afile) THEN
219  ierr = 0
220  CALL getd3dcur(aeqdskfile ,1,1,ierr)
221  IF(ierr .ne. 0) stop 'getd3dcur'
222  ELSE
223  write(67, fmt='('' CURTOR = '',1pe14.6)')g1%cpasma
224  ENDIF
225  CALL pgask(.true.)
226 CDESCUR j=1,itht1 points g2%rs(j,1),float(ij),g2%zs(j,1)
227  m=g2%npsi
228  ALLOCATE(s(m),yfit(m),yfit2(m),xin(m),yin(m),b4(mp))
229  ALLOCATE(thinner(m)); thinner=0
230  ALLOCATE(xpl(m), xpl2(m), ypl(m), ypl2(m))
231  s = g2%phi/g2%phi(m)
232  CALL pgpap(7.9,1.)
233  if(index(trim(device),'p')/=0) CALL pgpap(8.,1.)
234  CALL pgsubp(1,1)
235 ! if(INDEX(TRIM(device),'p')/=0) CALL pgsubp(2,2)
236  xlabel = 'rho'
237  xlabel = "\gr \(2240) \(2255)"// "s"
238  nfit = 11;cutin = cutoff*0.03_dbl;nlo = 0;b4 = 0
239  if(kkai > 0 .AND. kkai < nfit) nfit = kkai
240  yin = 1./g2%qsi
241  WRITE(porder,fmt='(i2.2)')nfit
242  mlabel = 'mapping & svd['//porder//'] of iota(s) vs sqrt(s)'
243  IF (kkai .EQ. 0) mlabel = trim(mlabel)//" - NOT USED"
244  ylabel = '1/qsi'
245  IF(.not.ALLOCATED(a0)) ALLOCATE(a0(nlo+1));a0 = 0.
246  CALL svdfit(nfit,nlo,cutin,a0,s,yin,m,b4)
247  yfit => polyval(b4,SIZE(b4),s,SIZE(s))
248  xpl = sqrt(s); ypl = yin;ypl2 = yfit
249  WRITE(slabel,fmt='("eqdsk = ",a," ; <% error>=",es9.2)')
250  . trim(geqdskfile), 200*sum(abs(yin-yfit))/sum(abs(yin+yfit))/nfit
251  CALL graf2pt
252  & (xpl,xpl,ypl,ypl2,m,xlabel,ylabel,' ',mlabel,slabel)
253  formatstring='('' AI = '',1pe14.6,2(/,2x,5(x,1pe14.6)))'
254  WRITE(67,fmt=formatstring)b4
255  formats=
256  .'(" AI= $",/," ["'
257  formats=trim(formats)//',5(1pe11.4,",")," $",/," "'
258  formats=trim(formats)//'5(1pe11.4,",")'
259  formats=trim(formats)//',"$",/," ",1pe11.4,"]")'
260  open(138,file="fort."//geqdskfile(ipos+1:),form="formatted")
261  write(138,*)"; ",trim(slabel)
262  write(138,fmt=trim(formats))(b4(i),i=1,11)
263  write(*,fmt=trim(formats))(b4(i),i=1,11)
264  ylabel='d(curint)/ds'
265  m=g2%npsi;xin=g2%phi;yin=g2%curintp
266  nfit=mp;cutin=cutoff*0.03_dbl;nlo=0;b4=0
267  if(kkac > 0 .AND. kkac < nfit) nfit=kkac
268  IF(.not.ALLOCATED(a0)) ALLOCATE(a0(nlo+1));a0=0.
269  CALL svdfit(nfit,nlo,cutin,a0,s,yin,m,b4)
270  yfit => polyval(b4,SIZE(b4),s,SIZE(s))
271  xpl = sqrt(s); ypl = yin; ypl2 = yfit
272  WRITE(porder,fmt='(i2.2)')nfit
273  xlabel = "\gr \(2240) \(2255)"// "s"
274  mlabel =
275  & 'mapping & svd['//porder//'] of dI\dtor\u(s)/ds vs sqrt(s)'
276  IF ( kkac .EQ. 0) mlabel = trim(mlabel)//" - NOT USED"
277  WRITE(slabel,fmt='("eqdsk = ",a," ; <% error>=",es9.2)')
278  . trim(geqdskfile), 200*sum(abs(yin-yfit))/sum(abs(yin+yfit))/nfit
279  CALL graf2pt
280  & (xpl,xpl,ypl,ypl2,m,xlabel,ylabel,' ',mlabel,slabel)
281  ylabel = 'd(curint)/ds'
282  CALL pgsci(4)
283  CALL pgpt1(0.1,real(maxval(yin))/5,5)
284  CALL pgtext (0.18, real(maxval(yin))/5, "mapped values")
285  CALL pgsci(2)
286  CALL pgpt1(0.1,real(maxval(yin))/4,4)
287  CALL pgtext (0.18, real(maxval(yin))/4, "polynomial fit")
288  formatstring = '('' AC = '',es14.7,2(/,2x,5(x,es14.7)))'
289  WRITE(67,fmt=formatstring)b4
290  formats =
291  .'(" AC= $",/," ["'
292  formats=trim(formats)//',5(1pe11.4,",")," $",/," "'
293  formats=trim(formats)//'5(1pe11.4,",")'
294  formats=trim(formats)//',"$",/," ",1pe11.4,"]")'
295  write(138,*)"; ",trim(slabel)
296  write(138,fmt=trim(formats))(b4(i),i=1,11)
297  write(*,fmt=trim(formats))(b4(i)/b4(1),i=1,11)
298  write(*,*)"Ienclosed=",real(itor(size(itor)))
299  yfit2=>polyint(b4,SIZE(b4),s,SIZE(s)) ! save integral
300  ylabel='curint'
301  m=g2%npsi;xin=g2%phi;yin=g2%curint
302  nfit=mp;cutin=cutoff*0.03_dbl;nlo=1;b4=0
303  IF(.not.ALLOCATED(a0)) ALLOCATE(a0(nlo+1));a0=0.
304  CALL svdfit(nfit,nlo,cutin,a0,s,yin,m,b4)
305  yfit => polyval(b4,SIZE(b4),s,SIZE(s))
306  xpl = sqrt(s); ypl = yin; ypl2 = yfit
307  WRITE(porder,fmt='(i2.2)')nfit
308  mlabel = 'mapping & svd['//porder//'] of Itor(s) vs sqrt(s)'
309  slabel = trim(filename)
310  WRITE(slabel,fmt='(a," DFF=",1pe12.6)')trim(geqdskfile),
311  . 2*sum(abs(yfit2-yfit))/sum(abs(yfit2+yfit))
312  CALL graf2pt
313  & (xpl,xpl,ypl,ypl2,m,xlabel,ylabel,' ',mlabel,slabel)
314  ypl=yfit2
315  call pgsci(8)
316  call pgslw(9)
317  if(index(trim(device),'p')/=0)call pgslw(2)
318  call pgsls(1)
319  call pgline(m,xpl,ypl)
320  mlabel="Solid Line is analytic integral of AC array"
321  call pgslw(1)
322  CALL pgmtxt('T',0.3724,0.5,0.5,trim(mlabel))
323  call pgslw(2)
324  ylabel='pres'
325  nfit=mp;cutin=cutoff*0.03_dbl;nlo=0;b4=0
326  if(kkam > 0 .AND. kkam < nfit) nfit=kkam
327  yin=g2%pressure
328  IF(.not.ALLOCATED(a0)) ALLOCATE(a0(nlo+1));a0=0.
329  CALL svdfit(nfit,nlo,cutin,a0,s,yin,m,b4)
330  yfit => polyval(b4,SIZE(b4),s,SIZE(s))
331  addv = sum(b4)
332  IF(addv .lt. 0 .AND. nfit == 11) THEN
333  b4 = b4-addv
334  ENDIF
335  yfit => polyval(b4,SIZE(b4),s,SIZE(s))
336  xpl = sqrt(s); ypl = yin; ypl2 = yfit
337  WRITE(porder,fmt='(i2.2)')nfit
338  mlabel = 'mapping & svd['//porder//'] of pressure(s) vs sqrt(s)'
339  IF (kkam .EQ. 0) mlabel = trim(mlabel)//" - NOT USED"
340  WRITE(slabel,fmt='("eqdsk = ",a," ; <% error>=",es9.2)')
341  . trim(geqdskfile), 200*sum(abs(yin-yfit))/sum(abs(yin+yfit))/nfit
342  CALL graf2pt
343  & (xpl,xpl,ypl,ypl2,m,xlabel,ylabel,' ',mlabel,slabel)
344  CALL pgsci(4)
345  CALL pgpt1(0.1,real(maxval(yin))/5,5)
346  CALL pgtext (0.18, real(maxval(yin))/5, "mapped values")
347  CALL pgsci(2)
348  CALL pgpt1(0.1,real(maxval(yin))/4,4)
349  CALL pgtext (0.18, real(maxval(yin))/4, "polynomial fit")
350  formatstring = '('' AM = '',es14.7,2(/,2x,5(x,es14.7)))'
351  WRITE(67,fmt=formatstring)b4
352 
353  formats =
354  .'(" AM= $",/," ["'
355  formats = trim(formats)//',5(1pe11.4,",")," $",/," "'
356  formats = trim(formats)//'5(1pe11.4,",")'
357  formats = trim(formats)//',"$",/," ",1pe11.4,"]")'
358  WRITE(138,*)"; ",trim(slabel)
359  WRITE(138,fmt=trim(formats))(b4(i),i=1,11)
360  WRITE(138,*)"p[0]=",b4(1)," & p[1]=",sum(b4)
361  WRITE(*,fmt=trim(formats))(b4(i)/b4(1),i=1,11)
362  WRITE(*,*)"p[0]=",b4(1)," & p[1]=",sum(b4)
363  IF ( kkam .EQ. 0 ) THEN
364  yin=g2%pressure
365  ier = 0
366  ylabel = 'pressure'
367  CALL doakima(s, yin, m, ier, trim(ylabel))
368  if(ier < 0) print*, '# knots :',-ier
369  formatstring = '(" pmass_type=''Akima_spline''")'
370  WRITE(67,fmt=trim(formatstring))
371  formatstring = '('' am_aux_s = '',es14.7,/,(2x,4es15.7)))'
372  WRITE(67,fmt=trim(formatstring))s2(1:n2)
373  formatstring = '('' am_aux_f = '',es14.7,/,(2x,4es15.7)))'
374  WRITE(67,fmt=trim(formatstring))f2(1:n2)
375  DEALLOCATE(s2,f2)
376  n2 = 0
377  ENDIF
378  IF ( kkac .EQ. 0 ) THEN
379  yin=g2%curint
380  ier = 0
381  ylabel = 'Itor'
382  CALL doakima(s, yin, m, ier, trim(ylabel))
383  if(ier < 0) print*, '# knots :', -ier
384  formatstring = '(" pcurr_type=''Akima_spline_Ip''")'
385  WRITE(67,fmt=trim(formatstring))
386  formatstring = '('' ac_aux_s = '',es14.7,/,(2x,4es15.7)))'
387  WRITE(67,fmt=trim(formatstring))s2(1:n2)
388  formatstring = '('' ac_aux_f = '',es14.7,/,(2x,4es15.7)))'
389  WRITE(67,fmt=trim(formatstring))f2(1:n2)
390  DEALLOCATE(s2,f2)
391  n2 = 0
392  ENDIF
393  IF ( kkai .EQ. 0 ) THEN
394  yin = 1./g2%qsi
395  ier = 0
396  ylabel = 'iota'
397  CALL doakima(s, yin, m, ier, trim(ylabel))
398  if(ier < 0) print*, '# knots :', -ier
399  formatstring = '(" piota_type=''Akima_spline''")'
400  WRITE(67,fmt=trim(formatstring))
401  formatstring = '('' ai_aux_s = '',es14.7,/,(2x,4es15.7)))'
402  WRITE(67,fmt=trim(formatstring))s2(1:n2)
403  formatstring = '('' ai_aux_f = '',es14.7,/,(2x,4es15.7)))'
404  WRITE(67,fmt=trim(formatstring))f2(1:n2)
405  DEALLOCATE(s2,f2)
406  n2 = 0
407  ENDIF
408  CALL pgsch(2.)
409  WRITE (67,fmt='(a,f8.4)') ' RAXIS = ',g1%rmaxis
410  WRITE (67,fmt='(a,f8.4)') ' ZAXIS = ',g1%zmaxis
411  CALL descur_sub(g1,g2)
412  CALL surfaces_plot(g1,g2,device)
413  CALL pgend
414  WRITE(67,fmt='(1h/)')
415  WRITE(67,fmt='(1h/)')
416  DO i=14,16
417  WRITE(67, fmt='(a)')trim(line(i))
418  ENDDO
419  WRITE(67,*)'!Coil ordering for VMEC differs from that for EFIT'
420  WRITE(67,*)'! 1: f1a 2: f1b'
421  WRITE(67,*)'! 3: f2a 4: f2b '
422  WRITE(67,*)'! 5: f3a 6: f3b '
423  WRITE(67,*)'! 7: f4a 8: f4b '
424  WRITE(67,*)'! 9: f5a 10: f5b '
425  WRITE(67,*)'! 11: f6a 12: f6b '
426  WRITE(67,*)'! 13: f7a 14: f7b '
427  WRITE(67,*)'! 15: f8a 16: f8b '
428  WRITE(67,*)'! 17: f9a 18: f9b '
429  WRITE(67,*)'! 19: 1oR '
430  WRITE(67,*)'! 20: eca 21: ecb'
431  WRITE(67,*)'! 22: iu330 23: iu270'
432  WRITE(67,*)'! 24: iu210 25: iu150'
433  WRITE(67,*)'! 26: iu90 27: iu30'
434  WRITE(67,*)'! 28: il330 29: il270'
435  WRITE(67,*)'! 30: il210 31: il150'
436  WRITE(67,*)'! 32: il90 33: il30'
437  WRITE(67,*)'! 34: ccoil79'
438  WRITE(67,*)'! 35: ccoil139'
439  WRITE(67,*)'! 36: ccoil199'
440 
441  formats = '(" end")'
442  WRITE(138,fmt=trim(formats))
443 46 stop 'normal termination'
444  END PROGRAM gtovmi2
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17