V3FIT
profil1d.f
1  SUBROUTINE profil1d_par(xc, xcdot, lreset)
2  USE vmec_main
3  USE vmec_params, ONLY: signgs, lamscale, rcc, pdamp
4  USE vmec_input, ONLY: lrfp
5  USE vspline
6  USE init_geometry, ONLY: lflip
7  USE vmec_input, ONLY: nzeta
8  USE vmec_dim, ONLY: ns, ntheta3
9  USE realspace
10  USE vmec_params, ONLY: ntmax
11  USE parallel_include_module
12  IMPLICIT NONE
13 C-----------------------------------------------
14 C D u m m y A r g u m e n t s
15 C-----------------------------------------------
16  REAL(dp), DIMENSION(0:ntor,0:mpol1,ns,3*ntmax), INTENT(OUT) ::
17  1 xc, xcdot
18  LOGICAL, INTENT(IN) :: lreset
19 
20 C-----------------------------------------------
21 C L o c a l P a r a m e t e r s
22 C-----------------------------------------------
23  REAL(dp), PARAMETER :: c1p5 = 1.5_dp
24 C-----------------------------------------------
25 C L o c a l V a r i a b l e s
26 C-----------------------------------------------
27  INTEGER :: i
28  REAL(dp) :: Itor, si, tf, pedge, vpnorm, polflux_edge
29  REAL(dp) :: phipslocal, phipstotal
30  INTEGER :: j, k, l, nsmin, nsmax
31  REAL (dp) :: torflux_edge
32 C-----------------------------------------------
33 C E x t e r n a l F u n c t i o n s
34 C-----------------------------------------------
35  REAL(dp), EXTERNAL :: pcurr, pmass, piota, torflux,
36  1 torflux_deriv, polflux, polflux_deriv
37 #ifdef _ANIMEC
38  2 , photp, ptrat
39 #endif
40 C-----------------------------------------------
41 !
42 ! INDEX OF LOCAL VARIABLES
43 !
44 ! ai array of coefficients in phi-series for iota (ncurr=0)
45 ! ac array of coefficients in phi-series for the quantity d(Icurv)/ds = toroidal
46 ! current density * Vprime, so Icurv(s) = Itor(s) (used for ncurr=1)
47 ! am array of coefficients in phi-series for mass (NWT/m**2)
48 ! iotas rotational transform , on half radial mesh
49 ! Icurv (-)toroidal current inside flux surface (vanishes like s)
50 ! mass mass profile on half-grid
51 ! phiedge value of real toroidal flux at plasma edge (s=1)
52 ! phips toroidal flux (same as phip), one-dimensional array
53 ! chips poloidal flux (same as chip), one-dimensional array
54 ! presf pressure profile on full-grid, mass/phip**gamma
55 ! spres_ped value of s beyond which pressure profile is flat (pedestal)
56 
57 !
58 ! COMPUTE PHIP, IOTA PROFILES ON FULL-GRID
59 ! COMPUTE MASS PROFILE ON HALF-GRID
60 ! BY READING INPUT COEFFICIENTS. PRESSURE CONVERTED TO
61 ! INTERNAL UNITS BY MULTIPLICATION BY mu0 = 4*pi*10**-7
62 !
63 
64  IF (ncurr.EQ.1 .AND. lrfp) THEN
65  stop 'ncurr=1 inconsistent with lRFP=T!'
66  END IF
67  torflux_edge = signgs*phiedge/twopi
68  si = torflux(one)
69  IF (si .ne. zero) THEN
70  torflux_edge = torflux_edge/si
71  END IF
72  polflux_edge = torflux_edge
73  si = polflux(one)
74  IF (si .ne. zero) THEN
75  polflux_edge = polflux_edge/si
76  END IF
77  r00 = rmn_bdy(0,0,rcc)
78 
79  phips(1) = 0
80  chips(1) = 0
81  icurv(1) = 0
82 
83  nsmin = max(2, t1lglob)
84  nsmax = t1rglob
85  DO i = nsmin, nsmax
86  si = hs*(i - c1p5)
87  tf = min(one, torflux(si))
88  IF (lrfp) THEN
89  tf = si
90  END IF
91  phips(i) = torflux_edge*torflux_deriv(si)
92  chips(i) = torflux_edge*polflux_deriv(si)
93  iotas(i) = piota(tf)
94  icurv(i) = pcurr(tf)
95  END DO
96 
97 !
98 ! Compute lamscale factor for "normalizing" lambda (needed for scaling hessian)
99 ! DO IT THIS WAY (rather than ALLREDUCE) FOR PROCESSOR INDEPENDENCE
100 !
101  CALL gather1xarray(phips)
102  phipstotal = sum(phips(2:ns)**2)
103  lamscale = sqrt(hs*phipstotal)
104 
105  IF (lflip) THEN
106  iotas = -iotas
107  chips = -chips
108  END IF
109 
110  nsmin = t1lglob
111  nsmax = t1rglob
112 
113  DO i = nsmin, nsmax
114  si = hs*(i - 1)
115  tf = min(one, torflux(si))
116  IF (lrfp) THEN
117  tf = si
118  END IF
119  iotaf(i) = piota(tf)
120  phipf(i) = torflux_edge*torflux_deriv(si)
121  chipf(i) = torflux_edge*polflux_deriv(si)
122  END DO
123 !
124 ! SCALE CURRENT TO MATCH INPUT EDGE VALUE, CURTOR
125 ! FACTOR OF SIGNGS NEEDED HERE, SINCE MATCH IS MADE TO LINE
126 ! INTEGRAL OF BSUBU (IN GETIOTA) ~ SIGNGS * CURTOR
127 !
128  pedge = pcurr(one)
129  itor = 0
130  IF (abs(pedge) .gt. abs(epsilon(pedge)*curtor)) THEN
131  itor = signgs*currv/(twopi*pedge)
132  END IF
133 
134  nsmin = max(2, t1lglob)
135  nsmax = t1rglob
136  icurv(nsmin:nsmax) = itor*icurv(nsmin:nsmax)
137 
138 !
139 ! POSSIBLE PRESSURE PEDESTAL FOR S >= SPRES_PED
140 !
141  spres_ped = abs(spres_ped)
142  IF (.not.lrecon) THEN
143  nsmin = max(2,t1lglob)
144  nsmax = t1rglob
145  DO i = nsmin, nsmax
146  si = hs*(i - c1p5)
147 
148 ! NORMALIZE mass so dV/dPHI (or dV/dPSI) in pressure to mass relation
149 ! See line 195 of bcovar: pres(2:ns) = mass(2:ns)/vp(2:ns)**gamma
150 
151  IF (lrfp) THEN
152  tf=si
153  vpnorm = polflux_edge*polflux_deriv(si)
154  ELSE
155  tf = min(one, torflux(si))
156  vpnorm = torflux_edge*torflux_deriv(si)
157  END IF
158  IF (si .gt. spres_ped) THEN
159  pedge = pmass(spres_ped)
160  ELSE
161  pedge = pmass(tf)
162  END IF
163  mass(i) = pedge*(abs(vpnorm)*r00)**gamma
164 #ifdef _ANIMEC
165 ! ANISOTROPIC PRESSURE, Tper/T|| RATIOS
166  phot(i) = photp(tf)
167  tpotb(i) = ptrat(tf)
168 #endif
169  END DO
170 
171  ELSE
172  nsmin = t1lglob
173  nsmax = t1rglob
174  iotas(nsmin:nsmax) = 0
175  iotaf(nsmin:nsmax) = 0
176  mass(nsmin:nsmax) = 0
177  presf(nsmin:nsmax) = 0
178  END IF
179 
180 
181  nsmin = t1lglob
182  nsmax = min(t1rglob, ns + 1)
183  pres(nsmin:nsmax) = 0
184  xcdot(:,:,nsmin:nsmax,:) = 0
185 
186 #ifdef _ANIMEC
187  medge = pmass(one)*(abs(phips(ns))*r00)**gamma
188  phedg = photp(one)
189 #endif
190  nsmin = max(1, t1lglob - 1)
191  nsmax = min(t1rglob + 1,ns)
192  DO i = nsmin, nsmax
193  si = hs*abs(i - 1.5_dp)
194  pshalf(:,i) = sqrt(si)
195  si = hs*(i - 1)
196  psqrts(:,i) = sqrt(si)
197  bdamp(i) = 2*pdamp*(1 - si)
198  END DO
199 
200  psqrts(:,ns) = 1 !!Avoid round-off
201 
202  nsmin = max(2, t1lglob)
203  nsmax = t1rglob
204  DO i = nsmin, nsmax
205  sm(i) = pshalf(1,i)/psqrts(1,i)
206  IF (i .LT. ns) THEN
207  sp(i) = pshalf(1,i + 1)/psqrts(1,i)
208  ELSE
209  sp(i)=one/psqrts(1,i)
210  END IF
211  END DO
212 
213 
214  sm(1) = 0
215  sp(0) = 0
216  sp(1) = sm(2)
217 
218  IF (lreset) THEN
219  xc(:,:,t1lglob:t1rglob,:) = 0
220  END IF
221 
222  END SUBROUTINE profil1d_par
223 
224  SUBROUTINE profil1d(xc, xcdot, lreset)
225  USE vmec_main
226  USE vmec_params, ONLY: signgs, lamscale, rcc, pdamp
227  USE vmec_input, ONLY: lrfp
228  USE vspline
229  USE realspace, ONLY: shalf, sqrts
230  USE init_geometry, ONLY: lflip
231 
232  IMPLICIT NONE
233 C-----------------------------------------------
234 C D u m m y A r g u m e n t s
235 C-----------------------------------------------
236  REAL(dp), DIMENSION(neqs), INTENT(out) :: xc, xcdot
237  LOGICAL, INTENT(IN) :: lreset
238 C-----------------------------------------------
239 C L o c a l P a r a m e t e r s
240 C-----------------------------------------------
241  REAL(dp), PARAMETER :: c1p5 = 1.5_dp
242 C-----------------------------------------------
243 C L o c a l V a r i a b l e s
244 C-----------------------------------------------
245  INTEGER :: i
246  REAL(dp) :: Itor, si, tf, pedge, vpnorm, polflux_edge
247  REAL(dp) :: torflux_edge
248 C-----------------------------------------------
249 C E x t e r n a l F u n c t i o n s
250 C-----------------------------------------------
251  REAL(dp), EXTERNAL :: pcurr, pmass, piota, torflux,
252  1 torflux_deriv, polflux, polflux_deriv
253 #ifdef _ANIMEC
254  2 , photp, ptrat
255 #endif
256 C-----------------------------------------------
257 !
258 ! INDEX OF LOCAL VARIABLES
259 !
260 ! ai array of coefficients in phi-series for iota (ncurr=0)
261 ! ac array of coefficients in phi-series for the quantity d(Icurv)/ds = toroidal
262 ! current density * Vprime, so Icurv(s) = Itor(s) (used for ncurr=1)
263 ! am array of coefficients in phi-series for mass (NWT/m**2)
264 ! iotas rotational transform , on half radial mesh
265 ! Icurv (-)toroidal current inside flux surface (vanishes like s)
266 ! mass mass profile on half-grid
267 ! phiedge value of real toroidal flux at plasma edge (s=1)
268 ! phips toroidal flux (same as phip), one-dimensional array
269 ! chips poloidal flux (same as chip), one-dimensional array
270 ! presf pressure profile on full-grid, mass/phip**gamma
271 ! spres_ped value of s beyond which pressure profile is flat (pedestal)
272 
273 !
274 ! COMPUTE PHIP, IOTA PROFILES ON FULL-GRID
275 ! COMPUTE MASS PROFILE ON HALF-GRID
276 ! BY READING INPUT COEFFICIENTS. PRESSURE CONVERTED TO
277 ! INTERNAL UNITS BY MULTIPLICATION BY mu0 = 4*pi*10**-7
278 !
279  IF (ncurr.EQ.1 .AND. lrfp) THEN
280  stop 'ncurr=1 inconsistent with lRFP=T!'
281  END IF
282 
283  torflux_edge = signgs*phiedge/twopi
284  si = torflux(one)
285  IF (si .ne. zero) THEN
286  torflux_edge = torflux_edge/si
287  END IF
288  polflux_edge = torflux_edge
289  si = polflux(one)
290  IF (si .ne. zero) THEN
291  polflux_edge = polflux_edge/si
292  END IF
293  r00 = rmn_bdy(0,0,rcc)
294 
295  phips(1) = 0
296  chips(1) = 0
297  icurv(1) = 0
298 
299  DO i = 2, ns
300  si = hs*(i - c1p5)
301  tf = min(one, torflux(si))
302  IF (lrfp) THEN
303  tf = si
304  END IF
305  phips(i) = torflux_edge*torflux_deriv(si)
306  chips(i) = torflux_edge*polflux_deriv(si)
307  iotas(i) = piota(tf)
308  icurv(i) = pcurr(tf)
309  END DO
310 
311 !
312 ! Compute lamscale factor for "normalizing" lambda (needed for scaling hessian)
313 !
314  lamscale = sqrt(hs*sum(phips(2:ns)**2))
315 
316  IF (lflip) THEN
317  iotas = -iotas
318  chips = -chips
319  END IF
320 
321  DO i = 1,ns
322  si = hs*(i - 1)
323  tf = min(one, torflux(si))
324  IF (lrfp) THEN
325  tf = si
326  END IF
327  iotaf(i) = piota(tf)
328  phipf(i) = torflux_edge*torflux_deriv(si)
329  chipf(i) = torflux_edge*polflux_deriv(si)
330  ENDDO
331 !
332 ! SCALE CURRENT TO MATCH INPUT EDGE VALUE, CURTOR
333 ! FACTOR OF SIGNGS NEEDED HERE, SINCE MATCH IS MADE TO LINE
334 ! INTEGRAL OF BSUBU (IN GETIOTA) ~ SIGNGS * CURTOR
335 !
336  pedge = pcurr(one)
337  itor = 0
338  IF (abs(pedge) .gt. abs(epsilon(pedge)*curtor)) THEN
339  itor = signgs*currv/(twopi*pedge)
340  END IF
341  icurv(2:ns) = itor*icurv(2:ns)
342 
343 !
344 ! POSSIBLE PRESSURE PEDESTAL FOR S >= SPRES_PED
345 !
346  spres_ped = abs(spres_ped)
347  IF (.not.lrecon) THEN
348  DO i = 2, ns
349  si = hs*(i - c1p5)
350 
351 ! NORMALIZE mass so dV/dPHI (or dV/dPSI) in pressure to mass relation
352 ! See line 195 of bcovar: pres(2:ns) = mass(2:ns)/vp(2:ns)**gamma
353 
354  IF (lrfp) THEN
355  tf = si
356  vpnorm = polflux_edge*polflux_deriv(si)
357  ELSE
358  tf = min(one, torflux(si))
359  vpnorm = torflux_edge*torflux_deriv(si)
360  END IF
361  IF (si .gt. spres_ped) THEN
362  pedge = pmass(spres_ped)
363  ELSE
364  pedge = pmass(tf)
365  END IF
366  mass(i) = pedge*(abs(vpnorm)*r00)**gamma
367 #ifdef _ANIMEC
368 ! ANISOTROPIC PRESSURE, Tper/T|| RATIOS
369  phot(i) = photp(tf)
370  tpotb(i) = ptrat(tf)
371 #endif
372  END DO
373 
374  ELSE
375  iotas(:ns) = 0
376  iotaf(:ns) = 0
377  mass(:ns) = 0
378  presf(:ns) = 0
379  END IF
380 
381  pres(:ns + 1) = 0
382  xcdot(:neqs) = 0
383 
384 #ifdef _ANIMEC
385  medge = pmass(one)*(abs(phips(ns))*r00)**gamma
386  phedg = photp(one)
387 #endif
388  DO i = 1, ns
389  si = hs*abs(i - 1.5_dp)
390 ! si = torflux(si) !SPH060409: shalf = sqrt(s), NOT sqrt(phi(s))!
391  shalf(i:nrzt:ns) = sqrt(si)
392  si = hs*(i - 1)
393  sqrts(i:nrzt:ns) = sqrt(si)
394  bdamp(i) = 2*pdamp*(1 - si)
395  END DO
396 
397  sqrts(ns:nrzt:ns) = 1 !!Avoid round-off
398  shalf(nrzt + 1) = 1
399  sqrts(nrzt + 1) = 1
400 
401  DO i = 2,ns
402  sm(i) = shalf(i)/sqrts(i)
403  sp(i) = shalf(i+1)/sqrts(i)
404  END DO
405 
406  sm(1) = 0
407  sp(0) = 0
408  sp(1) = sm(2)
409 
410  IF (lreset) THEN
411  xc(:neqs) = 0
412  END IF
413 
414  END SUBROUTINE profil1d