V3FIT
getd3dcur.f
1 c********************************************************************
2 C taken from /c/efit/processor/summard/osf/summard.f 6/11/2002
3  subroutine getd3dcur(filenm,jj,ireflux,ierr)
4 c********************************************************************
5  use sparmd
6  use mcomd1
7  real*8 :: tavem(ntime),chipre(ntime),sibdrr(ntime)
8  character*50 histor,myname
9  character*50 filenm
10  character uday*10,qmflag*3,vernum*8
11  character(len=3) :: colnam(18), coilname(18)
12  integer ,dimension(18) :: mapcoil, indx
13 ! 1 f1a 2 f1b 3 f2a 4 f2b 5 f3a 6 f3b 7 f4a 8 f4b 9 f5a 10 f5b
14 ! 11 f6a 12 f6b 13 f7a 14 f7b 15 f8a 16 f8b 17 f9a 18 f9b 19 1oR 20 eca 21 ecb
15  data mapcoil/1,3,5,7,9,11,13,15,17,2,4,6,8,10,12,14,16,18/
16  data indx/1,10,2,11,3,12,4,13,5,14,6,15,7,16,8,17,9,18/
17  data colnam /'f1a','f1b','f2a','f2b','f3a','f3b','f4a','f4b',
18  . 'f5a','f5b','f6a','f6b','f7a','f7b','f8a','f8b','f9a','f9b'/
19  data coilname /
20  . 'f1a','f2a','f3a','f4a','f5a','f6a','f7a','f8a','f9a',
21  . 'f1b','f2b','f3b','f4b','f5b','f6b','f7b','f8b','f9b'/
22 
23  read(filenm,fmt='(x,i6)')mshot
24  open(unit=nhist,status='old',
25  . file=filenm,err=401)
26  myname=filenm
27  98 format (32hphys_data:[d3phys.diiid.efitd65],a13)
28  100 format (25hphys_data:[eqdsk_d3.work],a13)
29  102 format (26hphys_data:[eqdsk_d3.workm],a13)
30 c-----------------------------------------------------------------
31 c-- ascii format --
32 c-----------------------------------------------------------------
33  read (nhist,1055,err=4000,end=4000) uday,(mfvers(j),j=1,2)
34  ltime = 0
35  if (mshot.le.99999) then
36  read (nhist,1050,err=4000,end=4000) ishot,ltime
37  else
38  read (nhist,1053,err=4000,end=4000) ishot,ltime
39  endif
40  vernum=mfvers(2)(2:5)//mfvers(1)(1:2)//mfvers(1)(4:5)
41  read(vernum,fmt=6000)nvernum
42  6000 format (i8)
43  read (nhist,1040) (xtime,i=1,ltime)
44  1060 format (x,f7.2,10x,i5,11x,i5,1x,a3,1x,i3,1x,i3,1x,a3)
45  read (nhist,1060) xtime(jj),jflag,lflag,
46  . limloc(jj),mco2v,mco2r,qmflag
47  read (nhist,1040) tsaisq(jj),rcentr,bcentr(jj),pasmat(jj)
48  read (nhist,1040) cpasma(jj),rout(jj),zout(jj),aout(jj)
49  read (nhist,1040) eout(jj),doutu(jj),doutl(jj),vout(jj)
50  read (nhist,1040) rcurrt(jj),zcurrt(jj),qsta(jj),betat(jj)
51  read (nhist,1040) betap(jj),ali(jj),oleft(jj),oright(jj)
52  read (nhist,1040) otop(jj),obott(jj),qpsib(jj),vertn(jj)
53  read (nhist,1040) (rco2v(k,jj),k=1,mco2v)
54  read (nhist,1040) (dco2v(jj,k),k=1,mco2v)
55  read (nhist,1040) (rco2r(k,jj),k=1,mco2r)
56  read (nhist,1040) (dco2r(jj,k),k=1,mco2r)
57  read (nhist,1040) shearb(jj),bpolav(jj),s1(jj),s2(jj)
58  read (nhist,1040) s3(jj),qout(jj),olefs(jj),orighs(jj)
59  read (nhist,1040) otops(jj),sibdry(jj),areao(jj),wplasm(jj)
60  read (nhist,1040) terror(jj),elongm(jj),qqmagx(jj),cdflux(jj)
61  read (nhist,1040) alpha(jj),rttt(jj),psiref(jj),xndnt(jj)
62  read (nhist,1040) rseps(1,jj),zseps(1,jj),rseps(2,jj)
63  . ,zseps(2,jj)
64  read (nhist,1040) sepexp(jj),obots(jj),btaxp(jj),btaxv(jj)
65  read (nhist,1040) aaq1(jj),aaq2(jj),aaq3(jj),seplim(jj)
66  read (nhist,1040) rmagx(jj),zmagx(jj),simagx(jj),taumhd(jj)
67  read (nhist,1040,err=380) betapd(jj),betatd(jj),
68  . wplasmd(jj),diamag(jj)
69  read (nhist,1040,err=380) vloopt(jj),taudia(jj),qmerci(jj),
70  . tavem(jj)
71  if (nvernum.ge.970524) then
72  read (nhist,1041,err=380) nsilop0,magpri0,nfcoil0,nesum0
73  else
74  nsilop0 = nsilop
75  nfcoil0 = nfcoil
76  nesum0 = nesum
77  if (ishot.lt.91000) then
78  magpri0 = magpri67+magpri322
79  else
80  magpri0 = magpri
81  endif
82  endif
83  read (nhist,1040,err=380) (csilop(k,jj),k=1,nsilop0),
84  . (cmpr2(k,jj),k=1,magpri0)
85  read (nhist,1040,err=380) (ccbrsp(k,jj),k=1,nfcoil0)
86  read (nhist,1040,err=380) (eccurt(jj,k),k=1,nesum0)
87 c
88  500 continue
89  go to 380
90 c-----------------------------------------------------------------
91 c-- binary format --
92 c-----------------------------------------------------------------
93  4000 continue
94  close(unit=nhist)
95  open(unit=nhist,status='old',
96  . file=myname,err=401,form='unformatted')
97  read (nhist,err=5000) uday,(mfvers(j),j=1,2)
98  vernum=mfvers(2)(2:5)//mfvers(1)(1:2)//mfvers(1)(4:5)
99  read(vernum,fmt=6000)nvernum
100  read (nhist,err=5000) ishot,ltime
101  read (nhist,err=5000) (xtime,i=1,ltime)
102  read (nhist,err=5000) xtime(jj),jflag,lflag,limloc(jj),
103  . mco2v,mco2r,qmflag
104  read (nhist,err=5000) tsaisq(jj),rcentr,bcentr(jj),pasmat(jj)
105  read (nhist,err=5000) cpasma(jj),rout(jj),zout(jj),aout(jj)
106  read (nhist,err=5000) eout(jj),doutu(jj),doutl(jj),vout(jj)
107  read (nhist,err=5000) rcurrt(jj),zcurrt(jj),qsta(jj),betat(jj)
108  read (nhist,err=5000) betap(jj),ali(jj),oleft(jj),oright(jj)
109  read (nhist,err=5000) otop(jj),obott(jj),qpsib(jj),vertn(jj)
110  read (nhist,err=5000) (rco2v(k,jj),k=1,mco2v)
111  read (nhist,err=5000) (dco2v(jj,k),k=1,mco2v)
112  read (nhist,err=5000) (rco2r(k,jj),k=1,mco2r)
113  read (nhist,err=5000) (dco2r(jj,k),k=1,mco2r)
114  read (nhist,err=5000) shearb(jj),bpolav(jj),s1(jj),s2(jj)
115  read (nhist,err=5000) s3(jj),qout(jj),olefs(jj),orighs(jj)
116  read (nhist,err=5000) otops(jj),sibdry(jj),areao(jj),wplasm(jj)
117  read (nhist,err=5000) terror(jj),elongm(jj),qqmagx(jj),
118  . cdflux(jj)
119  read (nhist,err=5000) alpha(jj),rttt(jj),psiref(jj),xndnt(jj)
120  read (nhist,err=5000) rseps(1,jj),zseps(1,jj),rseps(2,jj)
121  . ,zseps(2,jj)
122  read (nhist,err=5000) sepexp(jj),obots(jj),btaxp(jj),btaxv(jj)
123  read (nhist,err=5000) aaq1(jj),aaq2(jj),aaq3(jj),seplim(jj)
124  read (nhist,err=5000) rmagx(jj),zmagx(jj),simagx(jj),taumhd(jj)
125  read (nhist,err=5000) betapd(jj),betatd(jj),
126  . wplasmd(jj),diamag(jj)
127  read (nhist,err=5001) vloopt(jj),taudia(jj),qmerci(jj),
128  . tavem(jj)
129  if (nvernum.ge.970524) then
130  read (nhist,err=5001) nsilop0,magpri0,nfcoil0,nesum0
131  else
132  nsilop0 = nsilop
133  nfcoil0 = nfcoil
134  nesum0 = nesum
135  if (ishot.lt.91000) then
136  magpri0 = magpri67+magpri322
137  else
138  magpri0 = magpri
139  endif
140  endif
141  read (nhist,err=5001) (csilop(k,jj),k=1,nsilop0),
142  . (cmpr2(k,jj),k=1,magpri0)
143  read (nhist,err=5001) (ccbrsp(k,jj),k=1,nfcoil0)
144  read (nhist,err=5001) (eccurt(jj,k),k=1,nesum0)
145 c
146  5001 read (nhist,err=5000,end=380) header
147  5000 continue
148 c
149  1020 format (1hm,i5,1h.,i3)
150  1040 format (1x,4e16.9)
151  1041 format (1x,4i5)
152  1050 format (1x,i5,11x,i5)
153  1053 format (1x,i6,11x,i5)
154  1055 format (1x,a10,2a5)
155  380 continue
156  sibdrr(jj)=sibdry(jj)-csilop(nslref,jj)+psiref(jj)
157  if (ireflux.eq.1) then
158  sibdry(jj)=sibdry(jj)-csilop(nslref,jj)+psiref(jj)
159  simagx(jj)=simagx(jj)-csilop(nslref,jj)+psiref(jj)
160  psiref(jj)=csilop(nslref,jj)
161  endif
162  iok=1; ierr=0
163  390 close(unit=nhist)
164  400 continue
165 !!!!!!!!!!!!!!!!!!!!!!!!!!! VMEC order !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166 ! 1 f1a 2 f1b 3 f2a 4 f2b 5 f3a 6 f3b 7 f4a 8 f4b 9 f5a 10 f5b
167 !11 f6a 12 f6b 13 f7a 14 f7b 15 f8a 16 f8b 17 f9a 18 f9b
168 ! 19 1oR 20 eca 21 ecb
169 !!!!!!!!!!!!!!!!!!!!!!!!!!! VMEC order !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
170  k=1
171  kk=mapcoil(k)
172  kk9=mapcoil(k+9)
173  m9=indx(k+9)
174  m=indx(k)
175  write(67, fmt='('' CURTOR = '',1pe14.6)')cpasma(jj)
176  bcoil=bcentr(jj)*(rcentr/100)/2.e-7
177  terminal='(2('' EXTCUR('',i2.2,'') = '',1pe14.6))'
178  do k=1,nfcoil0/2
179  k9=k+9
180  kk=mapcoil(k)
181  kk9=mapcoil(k+9)
182  m9=indx(k+9)
183  m=indx(k)
184  write(67,fmt=terminal)
185  . kk,ccbrsp(k,jj)/turnfc(k),
186  . kk9,ccbrsp(k9,jj)/turnfc(k9)
187  enddo
188  terminal='('' EXTCUR('',i2.2,'') = '',1pe14.6)'
189  k=19;write(67,fmt=terminal)k,bcoil
190  terminal='(2('' EXTCUR('',i2.2,'')='',1pe13.6))'
191  k=20; write(67,fmt=terminal)
192  . k,eccurt(jj,k-19),k+1,eccurt(jj,k-19+3)
193 ! . k+2,eccurt(jj,k-19+2),
194 ! . k=20,20+nesum0/3)
195  return
196 ! debugging below here
197  write(123, fmt='('' CURTOR = '',1pe14.6)')cpasma(jj)
198  bcoil=bcentr(jj)*(rcentr/100)/2.e-7
199  terminal='(2('' EXTCUR('',i2.2,'') = '',1pe14.6))'
200  do k=1,nfcoil0/2
201  k9=k+9
202  kk=mapcoil(k)
203  kk9=mapcoil(k+9)
204  m9=indx(k+1)
205  m=indx(k)
206  write(123,fmt=terminal)
207  . kk,ccbrsp(k,jj)/turnfc(k),
208  . kk9,ccbrsp(k9,jj)/turnfc(k9)
209  write(*,fmt='(a2,6i3.2,2(x,a3))')
210  . 'k=',k,k9,kk,kk9,m,m9,colnam(kk),colnam(kk9)
211  write(*,fmt=terminal)
212  . kk,ccbrsp(k,jj)/turnfc(k),
213  . kk9,ccbrsp(k9,jj)/turnfc(k9)
214 
215  enddo
216  terminal='('' EXTCUR('',i2.2,'') = '',1pe14.6)'
217  k=19;write(123,fmt=terminal)k,bcoil
218  terminal='(2('' EXTCUR('',i2.2,'')='',1pe13.6))'
219  k=20; write(123,fmt=terminal)
220  . k,eccurt(jj,k-19),k+1,eccurt(jj,k-19+3)
221 ! . k+2,eccurt(jj,k-19+2),
222 ! . k=20,20+nesum0/3)
223  write(123,*)"mapcoil",mapcoil
224  write(123,*)"turns",turnfc
225  do i=1,18
226  write(123,*)i,coilname(i),turnfc(i),ccbrsp(i,jj)/turnfc((i))
227  enddo
228  write(123,*)'i,mapcoil(i),colnam(i),coilname(mapcoil(i))
229  .turnfc(indx(i)), indx(i),(ccbrsp(indx(i),jj)/turnfc(indx(i)))'
230  do i=1,18
231  k=mapcoil(i)
232  m=indx(i)
233  write(123,*)i,k,colnam(i),' ',
234  .coilname(m) ,turnfc(m), m,
235  .real(ccbrsp(m,jj)/turnfc(m))
236  enddo
237 
238  write(123,*)"ccbrsp(1:18,jj)",ccbrsp(1:18,jj)
239  write(123,*)"eccurt",eccurt
240  write(123,*)
241  .'1 f1a 2 f1b 3 f2a 4 f2b 5 f3a 6 f3b 7 f4a 8 f4b 9 f5a 10 f5b'
242  .,'11 f6a 12 f6b 13 f7a 14 f7b 15 f8a 16 f8b 17 f9a 18 f9b'
243  .,' 19 1oR 20 eca 21 ecb'
244 
245  return
246 401 continue
247  print*,'>>>>>>>>>> ERROR opening a-file = ',trim(filenm)
248  end subroutine getd3dcur
249