V3FIT
funct.f
1  SUBROUTINE funct(r0c, z0c, rhoc, rhos, xpts, gr0c, gz0c, grhoc,
2  1 grhos, gpts, fsq, xin, yin, mrho_in)
3 C-----------------------------------------------
4 C M o d u l e s
5 C-----------------------------------------------
6  USE vname0
7  USE vname1
8  IMPLICIT NONE
9 C-----------------------------------------------
10 C D u m m y A r g u m e n t s
11 C-----------------------------------------------
12  INTEGER mrho_in
13  REAL(rprec) fsq, r0c, z0c, gr0c, gz0c
14  REAL(rprec), DIMENSION(ntheta) :: xpts, gpts, xin, yin
15  REAL(rprec), DIMENSION(0:mrho_in-1) ::
16  1 rhoc, rhos, grhoc, grhos
17 C-----------------------------------------------
18 C L o c a l V a r i a b l e s
19 C-----------------------------------------------
20  INTEGER :: m, mrho1
21  REAL(rprec),DIMENSION(ntheta) ::
22  1 gcon, gtt, r1, z1, rt1, zt1
23  REAL(rprec), DIMENSION(ntheta,0:mrho_in) :: cosa, sina
24  REAL(rprec) :: denom, rmc_p, zms_p, rms_p, zmc_p,
25  1 t1, t2, tnorm
26 C-----------------------------------------------
27 !
28 !
29 ! FORCES: dW/dRmn = SUM(i)[ R(i) - RIN(i)]*COS(m*u[i]) ....
30 ! dW/dZmn = SUM(i)[ Z(i) - ZIN(i)]*SIN(m*u[i]) ....
31 ! dW/du[i] = rt(i)*[R(i)-RIN(i)] + zt(i)*[Z(i) - ZIN(i)]
32 ! THE NORM ON THE ANGLE FORCE (dW/du[i]) FOLLOWS FROM NEWTON'S
33 ! LAW AND IS APPROXIMATELY GTT = RT**2 + ZT**2
34 !
35  denom = zero
36  specw = zero
37 
38  mrho1 = mrho_in-1
39  r1(:ntheta) = -xin(:ntheta)
40  z1(:ntheta) = -yin(:ntheta)
41  rt1(:ntheta) = zero
42  zt1(:ntheta) = zero
43  cosa(:,0) = one
44  cosa(:,1) = cos(xpts)
45  sina(:,0) = zero
46  sina(:,1) = sin(xpts)
47 
48 !
49 ! COMPUTE CURVE R1 = R - Rin, Z1 = Z - Zin FOR INPUT POINTS
50 ! NOTE DIMENSIONS: Rm(0:MRHO), Zm(0:MRHO), but RHO(0:MRHO-1)
51 !
52  DO m = 2, mrho_in
53  cosa(:,m) = cosa(:,m-1)*cosa(:,1) - sina(:,m-1)*sina(:,1)
54  sina(:,m) = sina(:,m-1)*cosa(:,1) + cosa(:,m-1)*sina(:,1)
55  END DO
56 
57  DO m = 0, mrho_in
58  CALL getrz(rmc_p,rms_p,zmc_p,zms_p,r0c,z0c,rhoc,rhos,
59  1 m,mrho_in)
60 !
61 ! COMPUTE SPECTRAL WIDTH OF CURVE (DIAGNOSTIC)
62 !
63  t2 = rmc_p*rmc_p + zmc_p*zmc_p + rms_p*rms_p + zms_p*zms_p
64  denom = denom + t2*xmpq(m,2)
65  specw = specw + xmpq(m,1)*t2
66 
67  gtt(:ntheta) = rmc_p*cosa(:ntheta,m) + rms_p*sina(:ntheta,m)
68  gcon(:ntheta) = zmc_p*cosa(:ntheta,m) + zms_p*sina(:ntheta,m)
69  r1(:ntheta) = r1(:ntheta) + gtt(:ntheta)
70  z1(:ntheta) = z1(:ntheta) + gcon(:ntheta)
71  rt1(:ntheta) = rt1(:ntheta) + dm1(m)*(rms_p*cosa(:ntheta,m)-
72  1 rmc_p*sina(:ntheta,m))
73  zt1(:ntheta) = zt1(:ntheta) + dm1(m)*(zms_p*cosa(:ntheta,m)-
74  1 zmc_p*sina(:ntheta,m))
75  END DO
76 
77  specw = specw/denom
78 
79  gtt(:ntheta) = rt1(:ntheta)**2 + zt1(:ntheta)**2
80  gpts(:ntheta) = r1(:ntheta)*rt1(:ntheta)+z1(:ntheta)*zt1(:ntheta)
81 
82 !
83 ! COMPUTE MEAN-SQUARE DEVIATION BETWEEN POINTS AND FIT
84 !
85 ! t1 = MAXVAL(gtt(:ntheta))
86 ! gpts(:ntheta) = gpts(:ntheta)/t1
87  gpts(:ntheta) = 0.5_dp*gpts(:ntheta)/gtt(:ntheta)
88  t1 = maxval(abs(gpts(:ntheta)))
89  t2 = 1.e-3_dp !Arbitrary threshold on angle motion
90  IF (t1.gt.t2) THEN
91 ! Careful not to make angle points move too rapidly so they DO not cross
92  t1 = t2/t1
93  gpts(:ntheta) = t1 * gpts(:ntheta)
94  END IF
95  fsq = 0.5_dp*dnorm*sum(r1(:ntheta)**2 + z1(:ntheta)**2)
96  fsq = sqrt(fsq)/r10
97 
98 !
99 ! COMPUTE CURVE FORCES
100 ! 1. Magnetic axis forces: to artificially freeze axis at R0N,Z0N values
101 ! (set in inguess), set these to zero
102 !
103  m = 0
104  gr0c = dnorm*sum(cosa(:ntheta,m)*r1(:ntheta))
105  gz0c = dnorm*sum(cosa(:ntheta,m)*z1(:ntheta))
106 
107 !
108 ! RHO FORCES
109 ! (NOTE: RHO for m=0 IS EQUIVALENT TO THE L(v)-TERM => ARCLENGTH, Eq.14 H-B paper)
110 !
111  DO m = 0, mrho1
112  IF (m .le. 0) THEN
113  t1 = dnorm/max(t1m(m+1),0.1_dp)
114  grhoc(m) = t1*sum(cosa(:ntheta,m+1)*r1(:ntheta) +
115  1 sina(:ntheta,m+1)*z1(:ntheta))
116  grhos(m) = t1*sum(sina(:ntheta,m+1)*r1(:ntheta) -
117  1 cosa(:ntheta,m+1)*z1(:ntheta))
118  ELSE
119  t1 = t1m(m+1)
120  t2 = t2m(m-1)
121  IF (t1.EQ.0 .AND. t2.EQ.0) cycle
122  tnorm = dnorm/(t1*t1 + t2*t2)
123  t1 = t1*tnorm
124  t2 = t2*tnorm
125  grhoc(m) = sum((cosa(:ntheta,m+1)*r1(:ntheta) +
126  1 sina(:ntheta,m+1)*z1(:ntheta))*t1 +
127  2 (cosa(:ntheta,m-1)*r1(:ntheta) -
128  3 sina(:ntheta,m-1)*z1(:ntheta))*t2)
129  grhos(m) = sum((sina(:ntheta,m+1)*r1(:ntheta) -
130  1 cosa(:ntheta,m+1)*z1(:ntheta))*t1 +
131  2 (sina(:ntheta,m-1)*r1(:ntheta) +
132  3 cosa(:ntheta,m-1)*z1(:ntheta))*t2)
133  ENDIF
134  END DO
135 
136  grhos(0) = zero !Enforce constraint on toroidal angle
137 
138  gnorm = sum(grhoc(0:mrho1)*grhoc(0:mrho1) +
139  1 grhos(0:mrho1)*grhos(0:mrho1)) +
140  2 gr0c**2 + gz0c**2
141  gnorm = gnorm/r10**2
142 
143  gnorm = gnorm + dnorm*sum(gpts(:ntheta)*gpts(:ntheta))
144 
145  END SUBROUTINE funct