1 SUBROUTINE mapout_nc(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
29 CHARACTER(len=99) :: label
30 CHARACTER*63 cdffile,formatstring
40 CHARACTER(LEN=*),
PARAMETER ::
42 2 vn_nthet =
'nthet' ,
44 3 vn_percenflux =
'percenflux' ,
45 4 vn_psiaxis =
'psiaxis' ,
46 5 vn_psilim =
'psilim' ,
50 9 vn_pprime =
'pprime' ,
51 9 vn_vprime =
'vprime' ,
52 9 vn_press =
'press' ,
53 9 vn_iprime =
'iprime' ,
56 9 vn_volume =
'volume' ,
58 g vn_ffprime =
'ffprime',
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 ::
74 b vn_arcsur =
'arcsur',
75 k vn_eqdsk =
'eqdskname'
77 CHARACTER(LEN=*),
PARAMETER ::
78 1 ln_npsi =
'redial grid points',
79 2 ln_nthet =
'poloidal grid',
81 4 ln_psiaxis =
'poloidal at axis',
82 5 ln_psilim =
'poloidal at boundaryx ',
83 6 ln_kapn =
'normal curvature',
85 9 ln_xs =
'x on grid',
86 a ln_zs =
'y on grid',
87 b ln_arcsur =
'arc LENth (theta)'
95 FUNCTION deriv(xx,yy,n)
97 INTEGER,
PARAMETER :: rprec = selected_real_kind(12,100)
99 REAL(rprec) :: xx(n),yy(n)
100 REAL(rprec),
DIMENSION(n) :: x, y, dydx, x01, x02, x12
101 REAL(rprec) :: deriv(n)
109 INTEGER,
PARAMETER :: rprec = selected_real_kind(12,100)
112 REAL(rprec),
DIMENSION(n) :: x1
113 REAL(rprec) :: shiftx(n)
119 IF(len_trim(filename) .eq. 0) filename=
'eqdsk'
121 &filename(1+index(trim(filename),
"/",.true.):len_trim(filename))
129 twopsim1 = 2*npsi - 1
131 1 (path = trim(cdffile), cmode = nf90_clobber, ncid = ncid)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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,
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,
275 IF (status /= nf90_noerr) print*,
" AT:", location
276 IF (status /= nf90_noerr)
CALL handle_err1(status)
278 status = nf90_enddef(ncid)
280 IF (status /= nf90_noerr) print*,
" AT:", location
281 IF (status /= nf90_noerr)
CALL handle_err1(status)
299 jtor(j,k)=xs(j,k)*pprime(j)+mu0*ffprime(j)/xs(j,k)
302 dv=(psiv-psiv(1))/(psiv(
size(psiv))-psiv(1))
303 formatstring=
'(a,1pe10.3,x,1pe10.3)'
311 jtor(j,k)=xs(j,k)*pprime(j)*twopi+twopi/mu0*ffprime(j)/xs(j,k)
315 rmid=(/xs(row:1:-1,col/2),xs(2:row,1)/)
316 jmid=(/jtor(row:1:-1,col/2),jtor(2:row,1)/)
319 call graf1pt(real(rmid),real(jmid),
size(rmid),
'R',
'J [A/m\u3\d]',
320 .
'Midplane Local Current Density',
"eqdsk="//trim(filename))
322 label=
"(Informatinal only.)"
323 CALL pgtext(1.2,0.,trim(label))
325 bpmid=(/bps(row:1:-1,col/2),bps(2:row,1)/)
328 dl(j)=(arcsur(j,col/2)-arcsur(j,col/2-1))
332 dlob(j)=dlob(j)+dl(j)/bps(j,k)
337 jovr(j)=jovr(j)+dl(j)*jtor(j,k)/bps(j,k)/xs(j,k)
340 jovr(1)=jtor(1,1)/xs(1,1)
344 dv=(vprime +
shiftx(vprime,
SIZE(vprime),1))/
345 . 2*(psiv-
shiftx(psiv,
SIZE(psiv),1))
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)
360 itor(j)=itor(j-1)+dv(j)*jovr(j)/twopi
363 iprime=
deriv(tflx/tflx(
SIZE(tflx)),itor,
SIZE(itor))
364 status = nf90_put_var(ncid, id_eqdsk,filename)
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)
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)
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)
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)
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)
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)
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)
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))
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)
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)
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)
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))
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))
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))
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)
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))
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)
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))
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))
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))
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))
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))
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))
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))
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))
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))
471 IF (status /= nf90_noerr) print*,
" AT:", location
472 IF (status /= nf90_noerr)
CALL handle_err1(status)
473 status = nf90_close(ncid)
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
487 g2%area=g2%vplas/g2%rcentr
492 g2%aminor(j)=(g2%rs(j,1)-g2%rs(j,nthet/2))/2
497 g2%vprime=vprime/twopi
498 g2%pprime=pprime/twopi
499 g2%ffprime=(one/mu0)**2*ffprime*twopi
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)
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))
519 . ((xs(j,1)+xs(j,nthet/2))/2-(ru+rd)/2-g2%aminor(j))/g2%aminor(j)
521 allocate(w1(g2%npsi),rxinv(g2%npsi),zxsq(g2%npsi))
523 dx=arcsur(j,2)-arcsur(j,1)
527 rxinv(j)=rxinv(j)+dx/xs(j,i)
528 zxsq(j)=zxsq(j)+dx*zs(j,i)**2
530 rxinv(j)=rxinv(j)/arcsur(j,nthet)
531 zxsq(j)=zxsq(j)/arcsur(j,nthet)
549 &((err.le.errlast) .and. (square.lt.1.) .and. (iterg.lt.100))
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
569 err=100*abs(cntrlngth_z/cntrlngth-zxsq(j))/zxsq(j)
571 g2%square(j)=squarelast
581 END SUBROUTINE mapout_nc
582 subroutine handle_err1(status)
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)
590 INTEGER,
INTENT(in) :: status
591 INTEGER :: nf_noerr=0
592 IF (status .ne. nf_noerr)
THEN
595 10
FORMAT(
'% ',i4,
"--E-- A netCDF error has occurred in")
596 END SUBROUTINE handle_err12