V3FIT
mapout_nc.f
1  SUBROUTINE mapout_nc(g2)
2  USE netcdf
3 
4  USE mapg_mod
5  USE mapout
6  IMPLICIT NONE
7  TYPE(rszs) :: g2
8  REAL(rprec),dimension(:), allocatable :: w1,zxsq,rxinv
9  REAL(rprec) rsurf(kthet), zsurf(kthet), bpsurf(kthet), dloB(npsi)
10  REAL(rprec) spare1(kthet), spare2(kthet), lssurf(kthet)
11  REAL(rprec) dv(npsi), dl(npsi), xcentr(npsi), zcentr(npsi)
12  REAL(rprec) arclength, two, half, twopi, one, mtwopi
13  REAL(rprec) rcinv(kthet), y2(kthet), stemp, thetav, towb, mu0
14  REAL(rprec) rmid(2*npsi-1), jmid(2*npsi-1) , bpmid(2*npsi-1)
15  REAL(rprec) rmin, rmax, zmin, zmax, ru, rd, r0, z0, elong, a,
16  . triang, indent, square, squarelast, err, errlast, cntrlngth,
17  . cntrlngth_r, cntrlngth_z, t, rnow, znow, dg
18  INTEGER :: i, j, ntwopi, ier, ncid, k, row, col, pgopen, id0
19  INTEGER, DIMENSION(1) :: lu, ld
20  INTEGER, DIMENSION(2) :: idxid
21  INTEGER :: iterg, nstep, status, twopsim1, location=0
22  INTEGER :: PsiDimID, id_npsi, id_nthet, id_dpsi, id_percenflux
23  INTEGER :: TwoPsiM1ID, id_psiaxis, id_psilim, id_psiv, id_qsfin
24  INTEGER :: ThtDimID, id_jovR , id_tflx, id_pprime, id_rmid
25  INTEGER :: PsiThtID, id_vprime, id_press, id_iprime, id_volume
26  INTEGER :: id_eqdsk, id_ffprime, id_fval, id_itor, id_jmid
27  INTEGER :: id_bpmid, id_xs, id_zs, id_arcsur, id_bps, id_jtor
28  INTEGER :: id
29  CHARACTER(len=99) :: label
30  CHARACTER*63 cdffile,formatstring
31 c-----------------------------------------------------------------------
32 c 0.0 prepare to WRITE output to netcdf
33 c-----------------------------------------------------------------------
34 ! cdfttttttt(ncid,varnam, val [,ier]) --
35 ! LENgth and TYPE are inferred from val
36 ! cdf_setatt(ncid,varnam, longnam [,units='ZZ'[,ier]]) --
37 ! set long name attributes/units
38 ! cdf_WRITE(ncid,varnam,val[,ier]) -- WRITE variable val
39 c-------------- NetCDF labels & attributes -----------------------------
40  CHARACTER(LEN=*), PARAMETER ::
41  1 vn_npsi = 'npsi' ,
42  2 vn_nthet = 'nthet' ,
43  3 vn_dpsi = 'dpsi' ,
44  3 vn_percenflux = 'percenflux' ,
45  4 vn_psiaxis = 'psiaxis' ,
46  5 vn_psilim = 'psilim' ,
47  6 vn_psiv = 'pflx' ,
48  7 vn_qsf = 'qsf' ,
49  7 vn_tflx = 'tflx' ,
50  9 vn_pprime = 'pprime' ,
51  9 vn_vprime = 'vprime' ,
52  9 vn_press = 'press' ,
53  9 vn_iprime = 'iprime' , ! was bsquarav
54  9 vn_itor = 'itor' , ! was kappanav
55  9 vn_jovr = 'jovR' , ! was bsqprimeav
56  9 vn_volume = 'volume' , ! removed jacob
57  f vn_f = 'fpol',
58  g vn_ffprime = 'ffprime',
59  h vn_bp = 'bp',
60  h vn_rmid = 'rmid',
61  h vn_jmid = 'jmid',
62  h vn_bpmid = 'bpmid'
63 
64  CHARACTER(LEN=*), PARAMETER ::
65  6 ln_psiv = 'polodial flux' ,
66  7 ln_qsf = 'safety factor' ,
67  7 ln_tflx = 'toroidal flux' ,
68  9 ln_pprime = 'dp/dpsi' ,
69  9 ln_vprime = 'dV/dpsi'
70  CHARACTER(LEN=*), PARAMETER ::
71  8 vn_jtor = 'jtor',
72  9 vn_xs = 'xs',
73  a vn_zs = 'zs',
74  b vn_arcsur ='arcsur',
75  k vn_eqdsk = 'eqdskname'
76 
77  CHARACTER(LEN=*), PARAMETER ::
78  1 ln_npsi = 'redial grid points',
79  2 ln_nthet = 'poloidal grid',
80  3 ln_dpsi = 'dpsi',
81  4 ln_psiaxis = 'poloidal at axis',
82  5 ln_psilim = 'poloidal at boundaryx ',
83  6 ln_kapn = 'normal curvature',
84  8 ln_bsqrd = 'B^2',
85  9 ln_xs = 'x on grid',
86  a ln_zs = 'y on grid',
87  b ln_arcsur ='arc LENth (theta)'
88 c-----------------------------------------------------------------------
89 ! REAL(rprec), EXTERNAL, shiftx, deriv
90 ! REAL(rprec), EXTERNAL :: deriv
91 c-----------------------------------------------------------------------
92 c magnetic well quantities calculated in wellc
93 c-----------------------------------------------------------------------
94  INTERFACE deriv
95  FUNCTION deriv(xx,yy,n)
96  IMPLICIT NONE
97  INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
98  INTEGER :: n,n2
99  REAL(rprec) :: xx(n),yy(n)
100  REAL(rprec),DIMENSION(n) :: x, y, dydx, x01, x02, x12
101  REAL(rprec) :: deriv(n)
102  INTERFACE shiftx
103  END INTERFACE shiftx
104  END FUNCTION deriv
105  END INTERFACE deriv
106  INTERFACE shiftx
107  FUNCTION shiftx(x,n,m)
108  IMPLICIT NONE
109  INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
110  INTEGER :: m, n
111  REAL(rprec) :: x(n)
112  REAL(rprec),DIMENSION(n) :: x1
113  REAL(rprec) :: shiftx(n)
114  END FUNCTION shiftx
115  END INTERFACE shiftx
116 c
117 ! get units!
118 ! xs,zs are in cm, |B| is Tesla, p' is Pa/Wb/rad
119  IF(len_trim(filename) .eq. 0) filename='eqdsk'
120  cdffile="map_"//
121  &filename(1+index(trim(filename),"/",.true.):len_trim(filename))
122  &//".nc"
123 
124 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126 ! cdf_define here - complete ALL define's before ANY cdf_write:
127  ncid = 0
128  dpsi=1./(npsi-1.)
129  twopsim1 = 2*npsi - 1
130  status = nf90_create
131  1 (path = trim(cdffile), cmode = nf90_clobber, ncid = ncid)
132  location=location+1
133  IF (status /= nf90_noerr) print*," AT:", location
134  IF (status /= nf90_noerr) CALL handle_err1(status)
135  status = nf90_def_dim(ncid, "Psi", npsi, psidimid)
136  location=location+1
137  IF (status /= nf90_noerr) print*," AT:", location
138  IF (status /= nf90_noerr) CALL handle_err1(status)
139  status = nf90_def_dim(ncid, "TwoPsiM1", twopsim1, twopsim1id)
140  location=location+1
141  IF (status /= nf90_noerr) print*," AT:", location
142  IF (status /= nf90_noerr) CALL handle_err1(status)
143  status = nf90_def_dim(ncid, "Theta", nthet, thtdimid)
144  location=location+1
145  IF (status /= nf90_noerr) print*," AT:", location
146  IF (status /= nf90_noerr) CALL handle_err1(status)
147  idxid = (/psidimid, thtdimid/)
148  status = nf90_def_var(ncid,vn_eqdsk,nf90_char, id_eqdsk)
149  location=location+1
150  IF (status /= nf90_noerr) print*," AT:", location
151  IF (status /= nf90_noerr) CALL handle_err1(status)
152  status = nf90_def_var(ncid,vn_npsi, nf90_int, id_npsi)
153  location=location+1
154  IF (status /= nf90_noerr) print*," AT:", location
155  IF (status /= nf90_noerr) CALL handle_err1(status)
156  status = nf90_def_var(ncid,vn_nthet, nf90_int, id_nthet)
157  location=location+1
158  IF (status /= nf90_noerr) print*," AT:", location
159  IF (status /= nf90_noerr) CALL handle_err1(status)
160  status = nf90_def_var(ncid,vn_dpsi, nf90_double, id_dpsi)
161  location=location+1
162  IF (status /= nf90_noerr) print*," AT:", location
163  IF (status /= nf90_noerr) CALL handle_err1(status)
164  status = nf90_def_var
165  5 (ncid,vn_percenflux, nf90_double, id_percenflux)
166  location=location+1
167  IF (status /= nf90_noerr) print*," AT:", location
168  IF (status /= nf90_noerr) CALL handle_err1(status)
169  status = nf90_def_var(ncid,vn_psiaxis, nf90_double, id_psiaxis)
170  location=location+1
171  IF (status /= nf90_noerr) print*," AT:", location
172  IF (status /= nf90_noerr) CALL handle_err1(status)
173  status = nf90_def_var(ncid,vn_psilim, nf90_double, id_psilim)
174  location=location+1
175  IF (status /= nf90_noerr) print*," AT:", location
176  IF (status /= nf90_noerr) CALL handle_err1(status)
177  status = nf90_def_var
178  5 (ncid,vn_psiv, nf90_double, psidimid, id_psiv)
179  location=location+1
180  IF (status /= nf90_noerr) print*," AT:", location
181  IF (status /= nf90_noerr) CALL handle_err1(status)
182  status = nf90_def_var
183  5 (ncid,vn_qsf, nf90_double, psidimid, id_qsfin)
184  location=location+1
185  IF (status /= nf90_noerr) print*," AT:", location
186  IF (status /= nf90_noerr) CALL handle_err1(status)
187  status = nf90_def_var
188  5 (ncid,vn_tflx, nf90_double, psidimid, id_tflx)
189  location=location+1
190  IF (status /= nf90_noerr) print*," AT:", location
191  IF (status /= nf90_noerr) CALL handle_err1(status)
192  status = nf90_def_var
193  5 (ncid,vn_pprime, nf90_double, psidimid, id_pprime)
194  location=location+1
195  IF (status /= nf90_noerr) print*," AT:", location
196  IF (status /= nf90_noerr) CALL handle_err1(status)
197  status = nf90_def_var
198  5 (ncid,vn_vprime, nf90_double, psidimid, id_vprime)
199  location=location+1
200  IF (status /= nf90_noerr) print*," AT:", location
201  IF (status /= nf90_noerr) CALL handle_err1(status)
202  status = nf90_def_var
203  5 (ncid,vn_press, nf90_double, psidimid, id_press)
204  location=location+1
205  IF (status /= nf90_noerr) print*," AT:", location
206  IF (status /= nf90_noerr) CALL handle_err1(status)
207  status = nf90_def_var
208  5 (ncid,vn_iprime, nf90_double, psidimid, id_iprime)
209  location=location+1
210  IF (status /= nf90_noerr) print*," AT:", location
211  IF (status /= nf90_noerr) CALL handle_err1(status)
212  status = nf90_def_var
213  5 (ncid,vn_jovr, nf90_double, psidimid, id_jovr)
214  location=location+1
215  IF (status /= nf90_noerr) print*," AT:", location
216  IF (status /= nf90_noerr) CALL handle_err1(status)
217  status = nf90_def_var
218  5 (ncid,vn_volume, nf90_double, psidimid, id_volume)
219  location=location+1
220  IF (status /= nf90_noerr) print*," AT:", location
221  IF (status /= nf90_noerr) CALL handle_err1(status)
222  status = nf90_def_var
223  5 (ncid,vn_ffprime, nf90_double, psidimid, id_ffprime)
224  location=location+1
225  IF (status /= nf90_noerr) print*," AT:", location
226  IF (status /= nf90_noerr) CALL handle_err1(status)
227  status = nf90_def_var
228  5 (ncid,vn_f, nf90_double, psidimid, id_fval)
229  location=location+1
230  IF (status /= nf90_noerr) print*," AT:", location
231  IF (status /= nf90_noerr) CALL handle_err1(status)
232  status = nf90_def_var
233  5 (ncid,vn_itor, nf90_double, psidimid, id_itor)
234  location=location+1
235  IF (status /= nf90_noerr) print*," AT:", location
236  IF (status /= nf90_noerr) CALL handle_err1(status)
237  status = nf90_def_var
238  5 (ncid,vn_rmid, nf90_double, twopsim1id, id_rmid)
239  location=location+1
240  IF (status /= nf90_noerr) print*," AT:", location
241  IF (status /= nf90_noerr) CALL handle_err1(status)
242  status = nf90_def_var
243  5 (ncid,vn_jmid, nf90_double, twopsim1id, id_jmid)
244  location=location+1
245  IF (status /= nf90_noerr) print*," AT:", location
246  IF (status /= nf90_noerr) CALL handle_err1(status)
247  status = nf90_def_var
248  5 (ncid,vn_bpmid, nf90_double, twopsim1id, id_bpmid)
249  location=location+1
250  IF (status /= nf90_noerr) print*," AT:", location
251  IF (status /= nf90_noerr) CALL handle_err1(status)
252  status = nf90_def_var
253  5 (ncid,vn_xs, nf90_double, idxid ,id_xs)
254  location=location+1
255  IF (status /= nf90_noerr) print*," AT:", location
256  IF (status /= nf90_noerr) CALL handle_err1(status)
257  status = nf90_def_var
258  5 (ncid,vn_zs, nf90_double, idxid ,id_zs)
259  location=location+1
260  IF (status /= nf90_noerr) print*," AT:", location
261  IF (status /= nf90_noerr) CALL handle_err1(status)
262  status = nf90_def_var
263  5 (ncid,vn_bp, nf90_double, idxid , id_bps)
264  location=location+1
265  IF (status /= nf90_noerr) print*," AT:", location
266  IF (status /= nf90_noerr) CALL handle_err1(status)
267  status = nf90_def_var(ncid,vn_arcsur, nf90_double,
268  5 idxid ,id_arcsur)
269  location=location+1
270  IF (status /= nf90_noerr) print*," AT:", location
271  IF (status /= nf90_noerr) CALL handle_err1(status)
272  status = nf90_def_var(ncid,vn_jtor, nf90_double,
273  5 idxid , id_jtor)
274  location=location+1
275  IF (status /= nf90_noerr) print*," AT:", location
276  IF (status /= nf90_noerr) CALL handle_err1(status)
277 ! END DEFINITIONS
278  status = nf90_enddef(ncid)
279  location=location+1
280  IF (status /= nf90_noerr) print*," AT:", location
281  IF (status /= nf90_noerr) CALL handle_err1(status)
282 
283 !!!!!!
284 ! cdf_WRITE here - convert to MKS - complete ALL and CLOSE:
285 ! psi=2 pi psiv
286 ! d/dpsi = 1/[2 pi (psilim-psiaxis)]
287 
288  pi=acos(-1.)
289  two = 2
290  one=1
291  twopi = two*pi
292  half=one/two
293  mu0=pi*4e-7
294  row=npsi
295  col=nthet
296  jtor=0
297  DO j=1,row
298  DO k=1,col
299  jtor(j,k)=xs(j,k)*pprime(j)+mu0*ffprime(j)/xs(j,k)
300  ENDDO
301  ENDDO
302  dv=(psiv-psiv(1))/(psiv(size(psiv))-psiv(1))
303  formatstring='(a,1pe10.3,x,1pe10.3)'
304 ! print trim(formatstring),'minval(real(pprime 2 pi))=',
305 ! . minval(real(pprime*twopi))
306 
307 ! print trim(formatstring),'minval(real(ffprime u0))=',
308 ! . minval(real(ffprime*twopi/mu0))
309  DO j=1,row
310  DO k=1,col
311  jtor(j,k)=xs(j,k)*pprime(j)*twopi+twopi/mu0*ffprime(j)/xs(j,k)
312  ENDDO
313  ENDDO
314  jtor=-jtor/twopi
315  rmid=(/xs(row:1:-1,col/2),xs(2:row,1)/)
316  jmid=(/jtor(row:1:-1,col/2),jtor(2:row,1)/)
317  goto 99997 ! semi-permanent bypass
318  call pgsubp(1,1)
319  call graf1pt(real(rmid),real(jmid),size(rmid),'R','J [A/m\u3\d]',
320  . 'Midplane Local Current Density',"eqdsk="//trim(filename))
321  CALL pgsci(4)
322  label="(Informatinal only.)"
323  CALL pgtext(1.2,0.,trim(label))
324 99997 continue ! keep for debugging
325  bpmid=(/bps(row:1:-1,col/2),bps(2:row,1)/)
326  jovr=0; dl=0; dlob=0
327  DO j=2,row
328  dl(j)=(arcsur(j,col/2)-arcsur(j,col/2-1))
329  ENDDO
330  DO j=2,row ! cint(dl/B) vs psi
331  DO k=1,col
332  dlob(j)=dlob(j)+dl(j)/bps(j,k)
333  ENDDO
334  ENDDO
335  DO j=2,row ! cint(dl J/BR) vs psi
336  DO k=1,col
337  jovr(j)=jovr(j)+dl(j)*jtor(j,k)/bps(j,k)/xs(j,k)
338  ENDDO
339  ENDDO
340  jovr(1)=jtor(1,1)/xs(1,1)
341  jovr=jovr/dlob
342  dv=0
343  k=SIZE(psiv)
344  dv=(vprime +shiftx(vprime,SIZE(vprime),1))/
345  . 2*(psiv-shiftx(psiv,SIZE(psiv),1))
346  dv(1)=0
347 ! PRINT*,"V=",SUM(dv)
348  xcentr=0
349  zcentr=0
350  DO j=2,row
351  DO k=1,col
352  xcentr(j)=xcentr(j)+xs(j,k)*dl(j)/arcsur(j,col)
353  zcentr(j)=xcentr(j)+zs(j,k)*dl(j)/arcsur(j,col)
354  ENDDO
355  ENDDO
356  xcentr(1)=xs(1,1)
357  zcentr(1)=zs(1,1)
358  itor=0
359  DO j=2,row
360  itor(j)=itor(j-1)+dv(j)*jovr(j)/twopi
361  ENDDO
362 ! print formatstring,"Ienclosed=",real(itor(size(itor)))
363  iprime=deriv(tflx/tflx(SIZE(tflx)),itor,SIZE(itor))! d/ds NOT d/dphi
364  status = nf90_put_var(ncid, id_eqdsk,filename)
365  location=location+1
366  IF (status /= nf90_noerr) print*," AT:", location
367  IF (status /= nf90_noerr) CALL handle_err1(status)
368  status = nf90_put_var(ncid, id_npsi,npsi)
369  location=location+1
370  IF (status /= nf90_noerr) print*," AT:", location
371  IF (status /= nf90_noerr) CALL handle_err1(status)
372  status = nf90_put_var(ncid, id_nthet,nthet)
373  location=location+1
374  IF (status /= nf90_noerr) print*," AT:", location
375  IF (status /= nf90_noerr) CALL handle_err1(status)
376  status = nf90_put_var(ncid, id_dpsi,dpsi)
377  location=location+1
378  IF (status /= nf90_noerr) print*," AT:", location
379  IF (status /= nf90_noerr) CALL handle_err1(status)
380  status = nf90_put_var(ncid, id_percenflux,percenflux)
381  location=location+1
382  IF (status /= nf90_noerr) print*," AT:", location
383  IF (status /= nf90_noerr) CALL handle_err1(status)
384  status = nf90_put_var(ncid, id_psiaxis,psiaxis*twopi) ! Wb
385  location=location+1
386  IF (status /= nf90_noerr) print*," AT:", location
387  IF (status /= nf90_noerr) CALL handle_err1(status)
388  status = nf90_put_var(ncid, id_psilim,psilim*twopi) ! Wb
389  location=location+1
390  IF (status /= nf90_noerr) print*," AT:", location
391  IF (status /= nf90_noerr) CALL handle_err1(status)
392  status = nf90_put_var(ncid, id_psiv,psiv(1:npsi)*twopi) ! Wb
393  location=location+1
394  IF (status /= nf90_noerr) print*," AT:", location
395  IF (status /= nf90_noerr) CALL handle_err1(status)
396  status = nf90_put_var(ncid, id_qsfin,qsfin(1:npsi))
397  location=location+1
398  IF (status /= nf90_noerr) print*," AT:", location
399  IF (status /= nf90_noerr) CALL handle_err1(status)
400  status = nf90_put_var(ncid, id_tflx,tflx(1:npsi)*twopi)
401  location=location+1
402  IF (status /= nf90_noerr) print*," AT:", location
403  IF (status /= nf90_noerr) CALL handle_err1(status)
404  status = nf90_put_var(ncid, id_pprime,pprime(1:npsi)*twopi)
405  location=location+1
406  IF (status /= nf90_noerr) print*," AT:", location
407  IF (status /= nf90_noerr) CALL handle_err1(status)
408  status = nf90_put_var(ncid, id_vprime,vprime(1:npsi)/twopi)
409  location=location+1
410  IF (status /= nf90_noerr) print*," AT:", location
411  IF (status /= nf90_noerr) CALL handle_err1(status)
412  status = nf90_put_var(ncid, id_press,press(1:npsi))
413  location=location+1
414  IF (status /= nf90_noerr) print*," AT:", location
415  IF (status /= nf90_noerr) CALL handle_err1(status)
416  status = nf90_put_var(ncid, id_iprime,iprime(1:npsi))
417  location=location+1
418  IF (status /= nf90_noerr) print*," AT:", location
419  IF (status /= nf90_noerr) CALL handle_err1(status)
420  status = nf90_put_var(ncid, id_itor,itor(1:npsi))
421  location=location+1
422  IF (status /= nf90_noerr) print*," AT:", location
423  IF (status /= nf90_noerr) CALL handle_err1(status)
424  status = nf90_put_var(ncid, id_jovr,jovr(1:npsi)/twopi)
425  location=location+1
426  IF (status /= nf90_noerr) print*," AT:", location
427  IF (status /= nf90_noerr) CALL handle_err1(status)
428  status = nf90_put_var(ncid, id_volume,volume(1:npsi))
429  location=location+1
430  IF (status /= nf90_noerr) print*," AT:", location
431  IF (status /= nf90_noerr) CALL handle_err1(status)
432  status = nf90_put_var(ncid, id_ffprime,
433  & (one/mu0)**2*ffprime(1:npsi)*twopi)
434  location=location+1
435  IF (status /= nf90_noerr) print*," AT:", location
436  IF (status /= nf90_noerr) CALL handle_err1(status)
437  status = nf90_put_var(ncid, id_fval,one/mu0*fval(1:npsi))
438  location=location+1
439  IF (status /= nf90_noerr) print*," AT:", location
440  IF (status /= nf90_noerr) CALL handle_err1(status)
441  status = nf90_put_var(ncid, id_rmid,rmid(1:2*npsi-1))
442  location=location+1
443  IF (status /= nf90_noerr) print*," AT:", location
444  IF (status /= nf90_noerr) CALL handle_err1(status)
445  status = nf90_put_var(ncid, id_jmid,jmid(1:2*npsi-1))
446  location=location+1
447  IF (status /= nf90_noerr) print*," AT:", location
448  IF (status /= nf90_noerr) CALL handle_err1(status)
449  status = nf90_put_var(ncid, id_bpmid,bpmid(1:2*npsi-1))
450  location=location+1
451  IF (status /= nf90_noerr) print*," AT:", location
452  IF (status /= nf90_noerr) CALL handle_err1(status)
453  status = nf90_put_var(ncid, id_xs,xs(1:npsi,1:nthet))
454  location=location+1
455  IF (status /= nf90_noerr) print*," AT:", location
456  IF (status /= nf90_noerr) CALL handle_err1(status)
457  status = nf90_put_var(ncid, id_zs,zs(1:npsi,1:nthet))
458  location=location+1
459  IF (status /= nf90_noerr) print*," AT:", location
460  IF (status /= nf90_noerr) CALL handle_err1(status)
461  status = nf90_put_var(ncid, id_bps,bps(1:npsi,1:nthet))
462  location=location+1
463  IF (status /= nf90_noerr) print*," AT:", location
464  IF (status /= nf90_noerr) CALL handle_err1(status)
465  status = nf90_put_var(ncid, id_arcsur,arcsur(1:npsi,1:nthet))
466  location=location+1
467  IF (status /= nf90_noerr) print*," AT:", location
468  IF (status /= nf90_noerr) CALL handle_err1(status)
469  status = nf90_put_var(ncid, id_jtor,jtor(1:npsi,1:nthet))
470  location=location+1
471  IF (status /= nf90_noerr) print*," AT:", location
472  IF (status /= nf90_noerr) CALL handle_err1(status)
473  status = nf90_close(ncid)
474  location=location+1
475  IF (status /= nf90_noerr) print*," AT:", location
476  IF (status /= nf90_noerr) CALL handle_err1(status)
477  CALL rszs_allocate(npsi,nthet,g2)
478  g2%fraction_bndry=percenflux
479  g2%npsi=npsi
480  g2%nthet=nthet
481  g2%raxis=xs(1,1)
482  g2%zaxis=zs(1,1)
483  g2%rcentr=xcentr
484  g2%zcentr=zcentr
485  g2%vplas=0
486  g2%vplas=volume
487  g2%area=g2%vplas/g2%rcentr
488  g2%rs=xs
489  g2%zs=zs
490  g2%arcsur=arcsur
491  DO j=1,npsi
492  g2%aminor(j)=(g2%rs(j,1)-g2%rs(j,nthet/2))/2
493  ENDDO
494  g2%psival=psiv*twopi
495  g2%phi=tflx*twopi
496  g2%qsi=qsfin
497  g2%vprime=vprime/twopi
498  g2%pprime=pprime/twopi
499  g2%ffprime=(one/mu0)**2*ffprime*twopi
500  g2%fpol=one/mu0*fval
501  g2%curavg=jovr/twopi
502  g2%curint=itor
503  g2%curintp=iprime
504  g2%pressure=press
505  DO j=1,npsi
506  zmax=maxval(zs(j,1:nthet))
507  zmin=minval(zs(j,1:nthet))
508  rmax=maxval(xs(j,1:nthet))
509  rmin=minval(xs(j,1:nthet))
510  g2%elong(j)=(zmax-zmin)/2/g2%aminor(j)
511  lu(1)=maxloc(zs(j,1:nthet),1)
512  ld(1)=minloc(zs(j,1:nthet),1)
513  ru=xs(j,lu(1))
514  rd=xs(j,ld(1))
515  g2%triang(j)=min((g2%rcentr(j)-(ru+rd)/2)/g2%aminor(j),1.0_dbl)
516  ru=minval(xs(j,nthet/4:nthet/2))
517  rd=minval(xs(j,nthet/2:3*nthet/4))
518  g2%tnedni(j)=
519  . ((xs(j,1)+xs(j,nthet/2))/2-(ru+rd)/2-g2%aminor(j))/g2%aminor(j)
520  ENDDO
521  allocate(w1(g2%npsi),rxinv(g2%npsi),zxsq(g2%npsi))
522  do j=1,npsi
523  dx=arcsur(j,2)-arcsur(j,1)
524  rxinv(j)=0
525  zxsq(j)=0
526  do i=1,nthet
527  rxinv(j)=rxinv(j)+dx/xs(j,i)
528  zxsq(j)=zxsq(j)+dx*zs(j,i)**2
529  ENDDO
530  rxinv(j)=rxinv(j)/arcsur(j,nthet)
531  zxsq(j)=zxsq(j)/arcsur(j,nthet)
532  ENDDO
533 
534 !fit square
535  DO j=1,npsi
536  r0=g2%rcentr(j)
537  z0=g2%zcentr(j)
538  elong=g2%elong(j)
539  triang=g2%triang(j)
540  indent=g2%tnedni(j)
541  a=g2%aminor(j)
542  square=-.14
543  squarelast=square
544  errlast=1e10
545  err=1e9
546  iterg=0
547  nstep=1290
548  DO 10 WHILE
549  &((err.le.errlast) .and. (square.lt.1.) .and. (iterg.lt.100))
550  errlast=err
551  iterg=iterg+1
552  squarelast=square
553  square=square+.015
554  cntrlngth=0
555  cntrlngth_r=0
556  cntrlngth_z=0
557  DO i=1,nstep
558  t=float(i-1)*2*pi/(nstep-1)
559  dg=2*pi/(nstep-1)*sqrt(
560  & (-3*a*cos(t)**2*cos(t + sin(t)*triang)*sin(t)*indent +
561  & a*sin(t + sin(t)*triang)*(-1 - cos(t)*triang)*
562  & (1 + cos(t)**3*indent))**2 +
563  & a**2*cos(t + sin(t)*square)**2*(1 + cos(t)*square)**2*elong**2)
564  rnow=r0+a*cos(t+sin(t)*triang)*(1+cos(t)**3*indent)
565  znow=z0+a*elong*sin(t+square*sin(t))
566  cntrlngth_z=cntrlngth_z+dg*znow**2
567  cntrlngth=cntrlngth+dg
568  ENDDO
569  err=100*abs(cntrlngth_z/cntrlngth-zxsq(j))/zxsq(j)
570  10 CONTINUE
571  g2%square(j)=squarelast
572  ENDDO !j loop for squareness
573 
574 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
575 ! Here g1 & g2 structures have been filled. Once the code works
576 ! there should be a big deallocation that covers everything ALLOCATED in
577 ! mapg_MOD and in rdeqdsk
578 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
579 
580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581  END SUBROUTINE mapout_nc
582  subroutine handle_err1(status)
583  USE netcdf
584  integer, intent ( in) :: status
585  if(status /= nf90_noerr)
586  5 print *, trim(nf90_strerror(status))
587  end subroutine handle_err1
588  SUBROUTINE handle_err12(status)
589  IMPLICIT none
590  INTEGER, INTENT(in) :: status
591  INTEGER :: nf_noerr=0
592  IF (status .ne. nf_noerr) THEN
593  WRITE(*,10) status
594  ENDIF
595 10 FORMAT('% ',i4,"--E-- A netCDF error has occurred in")
596  END SUBROUTINE handle_err12
597 
598 !!!** Structure <40337470>, 18 tags, LENgth=805032, data LENgth=805032, refs=1:
599 !!! EQDSKNAME BYTE Array[40]
600 !!! NPSI LONG 129
601 !!! NTHT LONG 129
602 !!! DPSI DOUBLE 0.0078125000
603 !!! PERCENTFLUX DOUBLE 0.9995
604 !!! PSIAXIS DOUBLE -1.8060079
605 !!! PSILIM DOUBLE 0.20937294
606 !!! PSIV DOUBLE Array[129]
607 !!! QSF DOUBLE Array[129]
608 !!! TFLX DOUBLE Array[129]
609 !!! PPRIME DOUBLE Array[129]
610 !!! VPRIME DOUBLE Array[129]
611 !!! PRESS DOUBLE Array[129]
612 !!! iprime DOUBLE Array[129]
613 !!! itor DOUBLE Array[129]
614 !!! jovR DOUBLE Array[129]
615 !!! VOLUME DOUBLE Array[129]
616 !!! FFPRIME DOUBLE Array[129]
617 !!! FPOL DOUBLE Array[129]
618 !!! jtor DOUBLE Array[129, 129]
619 !!! XS DOUBLE Array[129, 129]
620 !!! ZS DOUBLE Array[129, 129]
621 !!! BP DOUBLE Array[129, 129]
622 !!! ARCSUR DOUBLE Array[129, 129]
623 
deriv
Definition: mapout_nc.f:94
shiftx
Definition: mapout_nc.f:102
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17