V3FIT
equalarc.f
1  SUBROUTINE equalarc(psivin,ipsi)
2  USE mapg_mod
3 c-----------------------------------------------------------------------
4 c stripped DOwn version of chali from mbc
5 c rlm 6/22/94
6 c this routine finds psi contour psivin=psiv(ipsi) and RETURN through
7 c common blocks the equal arc LENgth values of xs(ipsi,j) and zs(ipsi,j)
8 c arcsur(ipsi,j), the arc LENgth along the contour, is also RETURNed
9 c-----------------------------------------------------------------------
10  IMPLICIT NONE
11  REAL(rprec) :: psivin
12  INTEGER :: ipsi, icount=0
13 
14  INTEGER is,js,npc,ntw,npco,npcm1
15  REAL(rprec) :: delxz
16  REAL(rprec) :: xp(nxzd),zp(nxzd),bpsqp(nxzd),tp(nxzd)
17 c REAL(rprec) :: csx(3,nxzd),csz(3,nxzd),csbp(3,nxzd)
18  REAL(rprec) :: csx(nxzd),csz(nxzd),csbp(nxzd)
19  REAL(rprec) :: arc(nxzd)
20  REAL(rprec) :: dl,ds,bpsqval
21  INTEGER j,nlow
22 c REAL(rprec) :: sterpl,f2s,bb(4)
23 c INTEGER ier
24 c-----------------------------------------------------------------------
25 c the following common is shared with routine cntour
26 c-----------------------------------------------------------------------
27  REAL(rprec) :: xemin,xemax !output form cntour
28  REAL(rprec) :: yemin,yemax !output form cntour
29  REAL(rprec) :: yxmin,yxmax !output form cntour
30  REAL(rprec) :: xymin,xymax !output form cntour
31  REAL(rprec) :: xaxd,yaxd !input to cntour
32  REAL(rprec) :: dang,arcl,bperr !input to cntour
33  REAL(rprec) :: xmin,xmax !input to cntour
34  REAL(rprec) :: ymin,ymax !input to cntour
35  common/cntd/xaxd,yaxd,
36  $ xemin,xemax,yemin,yemax,yxmin,yxmax,xymin,xymax,
37  $ dang,arcl,bperr,xmin,xmax,ymin,ymax
38 c-----------------------------------------------------------------------
39 c attempt to find contour using furpl
40 c is and js are starting values for search for flux surfaces
41 c-----------------------------------------------------------------------
42  npc=1
43  ntw=nxzd
44  is=nx/2
45  js=nz/2
46  IF(ipsi .eq. 2222) THEN ! set to 2 for PRINTing
47  print*," shape(psixz) ",shape(psixz)
48  print*," shape(bpsq) ",shape(bpsq)
49  print*," shape(xgrid) ",shape(xgrid)
50  print*," shape(zgrid) ",shape(zgrid)
51  print*," shape(psivin) ",shape(psivin)
52  print*," shape(xp) ",shape(xp)
53  print*," shape(zp) ",shape(zp)
54  print*," shape(bpsqp) ",shape(bpsqp)
55  print*,"npc,nx,nz,ntw,nxd,nzd,is,js=",
56  & npc,nx,nz,ntw,nxd,nzd,is,js
57  ENDIF
58  CALL furplm(psixz,bpsq,xgrid,zgrid,psivin,xp,zp,bpsqp,npc,
59  $ nx,nz,ntw,nxd,nzd,is,js)
60  IF(ipsi .eq. 2222) THEN ! set to 2 for PRINTing
61  print*,npc,psivin, maxval(psixz),minval(psixz)
62  print*,(psivin > maxval(psixz)).and.(psivin > minval(psixz))
63  ENDIF
64 
65 c npc=0
66 c npfit=400
67  npfit=40
68  IF(npc.le.npfit) THEN
69 c WRITE(6,'("used contour",2i5)') ipsi,npc
70  xaxd=xaxis
71  yaxd=zaxis
72  npco=npc
73  xmin=xgrid(1)
74  xmax=xgrid(nx)
75  ymin=zgrid(1)
76  ymax=zgrid(nz)
77  bperr=0.1
78  IF(ipsi.eq.2) THEN
79  arcl=max(dx,dz)*0.1
80  ELSE
81  arcl=arcsur(ipsi-1,nthet)/npfit
82  ENDIF
83  npfit=400
84  dang=1./npfit
85  CALL cntourp(xgrid,nx,zgrid,nz,csplpsi,xp,zp,bpsqp,
86  $ npc,dx,dz,nxzd,psivin,nxd)
87  CALL sortr(xp,zp,bpsqp,npc,xaxis,zaxis,xaxis,zaxis,0,1)
88  ELSE
89  IF(xp(1).ne.xp(npc).or.zp(1).ne.zp(npc)) THEN
90  WRITE(6,'("error in equalarc finding contour ",
91  $ i3," of ",i3)') ipsi,npsi
92  WRITE(6,'("IF this is final contour,"
93  $ ," try reducing percenflux")')
94  stop
95  ENDIF
96  CALL sortr(xp,zp,bpsqp,npc,xaxis,zaxis,xaxis,zaxis,0,1)
97  delxz=0.20*min(xgrid(2)-xgrid(1),zgrid(2)-zgrid(1))
98  CALL packk(bpsqp,xp,zp,npc,delxz)
99  ENDIF
100  IF(npc.le.5) THEN
101  WRITE(6,'("error in equalarc--npc too small=")') npc
102  stop
103  ENDIF
104 c-----------------------------------------------------------------------
105 c set up arclengths on surface
106 c-----------------------------------------------------------------------
107 ! icount=icount+1; PRINT*,'HERE',icount
108  tp(1)=0.
109  npcm1=npc-1
110  DO j=2,npc
111  tp(j)=tp(j-1)+sqrt((xp(j)-xp(j-1))**2+(zp(j)-zp(j-1))**2)
112  END DO
113 c-----------------------------------------------------------------------
114 c spline xp wrt t
115 c-----------------------------------------------------------------------
116  CALL spline(tp,xp,npc,-1.e31_dbl,-1.e31_dbl,csx)
117 c-----------------------------------------------------------------------
118 c spline zp wrt t
119 c-----------------------------------------------------------------------
120  CALL spline(tp,zp,npc,-1.e31_dbl,-1.e31_dbl,csz)
121 c-----------------------------------------------------------------------
122 c calculate arclength to each point
123 c-----------------------------------------------------------------------
124  CALL garc(tp,xp,zp,csx,csz,nxzd,arc,npc)
125 c-----------------------------------------------------------------------
126 c spline xp wrt arc
127 c-----------------------------------------------------------------------
128  CALL spline(arc,xp,npc,-1.e31_dbl,-1.e31_dbl,csx)
129 c-----------------------------------------------------------------------
130 c spline zp wrt arc
131 c-----------------------------------------------------------------------
132  CALL spline(arc,zp,npc,-1.e31_dbl,-1.e31_dbl,csz)
133 c-----------------------------------------------------------------------
134 c spline bpsqp wrt arc
135 c-----------------------------------------------------------------------
136  CALL spline(arc,bpsqp,npc,-1.e31_dbl,-1.e31_dbl,csbp)
137 c-----------------------------------------------------------------------
138 c convert to equal arc LENgth along contour
139 c-----------------------------------------------------------------------
140  ds=arc(npc)/(nthet-1.)
141  DO j=1,nthet-1
142  dl=(j-1)*ds
143  nlow=0
144  CALL splint(arc,xp,csx,npc,dl,xs(ipsi,j),nlow)
145  CALL splint(arc,zp,csz,npc,dl,zs(ipsi,j),nlow)
146  CALL splint(arc,bpsqp,csbp,npc,dl,bpsqval,nlow)
147  IF(bpsqval.gt.0.) THEN
148  bps(ipsi,j)=sqrt(bpsqval)
149  ELSE
150  bps(ipsi,j)=0.
151  ENDIF
152  arcsur(ipsi,j)=dl
153  END DO
154  xs(ipsi,nthet)=xs(ipsi,1)
155  zs(ipsi,nthet)=zs(ipsi,1)
156  bps(ipsi,nthet)=bps(ipsi,1)
157  arcsur(ipsi,nthet)=arc(npc)
158  END SUBROUTINE equalarc
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17