V3FIT
mercier.f
1  SUBROUTINE mercier(gsqrt, bsq, bdotj, iotas, wint,
2  1 r1, rt, rz, zt, zz, bsubu, vp, phips, pres, ns, nznt)
3  USE safe_open_mod
4  USE vmercier
5  USE vmec_input, ONLY: input_extension
6  USE vparams, ONLY: one, zero, twopi, nmercier0
7  IMPLICIT NONE
8 C-----------------------------------------------
9 C D u m m y A r g u m e n t s
10 C-----------------------------------------------
11  INTEGER, INTENT(in) :: ns, nznt
12  REAL(rprec), DIMENSION(ns,nznt), INTENT(in) ::
13  1 gsqrt, bsq
14  REAL(rprec), DIMENSION(ns,nznt), INTENT(inout) :: bdotj
15  REAL(rprec), DIMENSION(ns*nznt), INTENT(in) :: wint, bsubu
16  REAL(rprec), DIMENSION(ns,nznt,0:1), INTENT(in) ::
17  1 r1, rt, rz, zt, zz
18  REAL(rprec), DIMENSION(ns), INTENT(in) ::
19  1 iotas, vp, phips, pres
20 C-----------------------------------------------
21 C L o c a l P a r a m e t e r s
22 C-----------------------------------------------
23  REAL(rprec), PARAMETER :: p5 = 0.5_dp, two = 2
24 C-----------------------------------------------
25 C L o c a l V a r i a b l e s
26 C-----------------------------------------------
27  INTEGER :: ns1, i, imercier0, nmerc = nmercier0, nrzt
28  REAL(rprec) :: sign_jac, hs, sqs, denom
29  REAL(rprec), DIMENSION(:,:), ALLOCATABLE ::
30  1 gpp, gsqrt_full, b2
31  REAL(rprec), DIMENSION(nznt) :: gtt, rtf, ztf,
32  1 rzf, zzf, r1f, jdotb, ob2, b2i
33  REAL(rprec), DIMENSION(ns) :: vp_real, phip_real,
34  1 shear, vpp, presp, torcur, ip, sj, tpp, tjb, tbb, tjj
35  CHARACTER(LEN=120) mercier_file
36 C-----------------------------------------------
37  nrzt = ns*nznt
38  mercier_file = 'mercier.'//trim(input_extension)
39  CALL safe_open (nmerc, imercier0, mercier_file, 'replace',
40  1 'formatted')
41  IF (imercier0 .ne. 0) RETURN
42 
43 
44  ALLOCATE (gpp(ns,nznt), gsqrt_full(ns,nznt), b2(ns,nznt), stat=i)
45  IF (i .ne. 0) stop 'allocation error in Mercier'
46 
47 !
48 ! SCALE VP, PHIPS TO REAL UNITS (VOLUME, TOROIDAL FLUX DERIVATIVES)
49 ! AND PUT GSQRT IN ABS UNITS (SIGNGS MAY BE NEGATIVE)
50 ! NOTE: VP has (coming into this routine) the sign of the jacobian multiplied out
51 ! i.e., vp = signgs*<gsqrt>
52 ! THE SHEAR TERM MUST BE MULTIPLIED BY THE SIGN OF THE JACOBIAN
53 ! (OR A BETTER SOLUTION IS TO RETAIN THE JACOBIAN SIGN IN ALL TERMS, INCLUDING
54 ! VP, THAT DEPEND EXPLICITLY ON THE JACOBIAN. WE CHOOSE THIS LATTER METHOD...)
55 !
56 ! COMING INTO THIS ROUTINE, THE JACOBIAN(gsqrt) = 1./(grad-s . grad-theta X grad-zeta)
57 ! WE CONVERT THIS FROM grad-s to grad-phi DEPENDENCE BY DIVIDING gsqrt by PHIP_real
58 !
59 ! NOTE: WE ARE USING 0 < s < 1 AS THE FLUX VARIABLE, BEING CAREFUL
60 ! TO KEEP d(phi)/ds == PHIP_real FACTORS WHERE REQUIRED
61 ! THE V'' TERM IS d2V/d(PHI)**2, PHI IS REAL TOROIDAL FLUX
62 !
63 ! SHEAR = d(iota)/d(phi) : FULL MESH
64 ! VPP = d(vp)/d(phi) : FULL MESH
65 ! PRESP = d(pres)/d(phi) : FULL MESH (PRES IS REAL PRES*mu0)
66 ! IP = d(Itor)/d(phi) : FULL MESH
67 !
68 ! ON ENTRY, BDOTJ = Jacobian * J*B ON THE FULL RADIAL GRID
69 ! BSQ = 0.5*|B**2| + p IS ON THE HALF RADIAL GRID
70 !
71 
72  ns1 = ns - 1
73  IF (ns1 .le. 0) RETURN
74  hs = one/ns1
75  sign_jac = zero
76  IF (gsqrt(ns,1) .ne. zero)
77  1 sign_jac = abs(gsqrt(ns,1))/gsqrt(ns,1)
78 
79  IF (sign_jac .eq. zero) RETURN
80  phip_real = twopi * phips * sign_jac
81 !
82 ! NOTE: phip_real should be > 0 to get the correct physical sign of REAL-space gradients
83 ! For example, grad-p, grad-Ip, etc. However, with phip_real defined this way,
84 ! Mercier will be correct
85 !
86  vp_real(2:ns) = sign_jac*(twopi*twopi)*vp(2:ns)/phip_real(2:ns) !!dV/d(PHI) on half mesh
87 
88 !
89 ! COMPUTE INTEGRATED TOROIDAL CURRENT
90 !
91  DO i = 2,ns
92  torcur(i)=sign_jac*twopi*sum(bsubu(i:nrzt:ns)*wint(i:nrzt:ns))
93  END DO
94 
95 !
96 ! COMPUTE SURFACE AVERAGE VARIABLES ON FULL RADIAL MESH
97 !
98  DO i = 2,ns1
99  phip_real(i) = p5*(phip_real(i+1) + phip_real(i))
100  denom = one/(hs*phip_real(i))
101  shear(i) = (iotas(i+1) - iotas(i))*denom !!d(iota)/d(PHI)
102  vpp(i) = (vp_real(i+1) - vp_real(i))*denom !!d(VP)/d(PHI)
103  presp(i) = (pres(i+1) - pres(i))*denom !!d(p)/d(PHI)
104  ip(i) = (torcur(i+1) - torcur(i))*denom !!d(Itor)/d(PHI)
105  END DO
106 
107 !
108 ! COMPUTE GPP == |grad-phi|**2 = PHIP**2*|grad-s|**2 (on full mesh)
109 ! Gsqrt_FULL = JACOBIAN/PHIP == jacobian based on flux (on full mesh)
110 !
111 
112  DO i = 2, ns1
113  gsqrt_full(i,:) = p5*(gsqrt(i,:) + gsqrt(i+1,:))
114  bdotj(i,:) = bdotj(i,:)/gsqrt_full(i,:)
115  gsqrt_full(i,:) = gsqrt_full(i,:)/phip_real(i)
116  sj(i) = hs*(i-1)
117  sqs = sqrt(sj(i))
118  rtf(:) = rt(i,:,0) + sqs*rt(i,:,1)
119  ztf(:) = zt(i,:,0) + sqs*zt(i,:,1)
120  gtt(:) = rtf(:)*rtf(:) + ztf(:)*ztf(:)
121  rzf(:) = rz(i,:,0) + sqs*rz(i,:,1)
122  zzf(:) = zz(i,:,0) + sqs*zz(i,:,1)
123  r1f(:) = r1(i,:,0) + sqs*r1(i,:,1)
124  gpp(i,:) = gsqrt_full(i,:)**2/(gtt(:)*r1f(:)**2 +
125  1 (rtf(:)*zzf(:) - rzf(:)*ztf(:))**2) !!1/gpp
126  END DO
127 
128 !
129 ! COMPUTE SURFACE AVERAGES OVER dS/|grad-PHI|**3 => |Jac| du dv / |grad-PHI|**2
130 ! WHERE Jac = gsqrt/phip_real
131 !
132  DO i = 2,ns
133  b2(i,:) = two*(bsq(i,:) - pres(i))
134  END DO
135 
136  DO i = 2,ns1
137  b2i(:) = p5*(b2(i+1,:) + b2(i,:))
138  ob2(:) = gsqrt_full(i,:)/b2i(:)
139  tpp(i) = sum(ob2(:)*wint(i:nrzt:ns))
140  ob2(:) = b2i(:) * gsqrt_full(i,:) * gpp(i,:)
141  tbb(i) = sum(ob2(:)*wint(i:nrzt:ns))
142  jdotb(:) = bdotj(i,:) * gpp(i,:) * gsqrt_full(i,:)
143  tjb(i) = sum(jdotb(:)*wint(i:nrzt:ns))
144  jdotb(:) = jdotb(:) * bdotj(i,:) / b2i(:)
145  tjj(i) = sum(jdotb(:)*wint(i:nrzt:ns))
146  END DO
147 
148  DEALLOCATE (gpp, gsqrt_full, b2, stat=i)
149 
150 !
151 ! REFERENCE: BAUER, BETANCOURT, GARABEDIAN, MHD Equilibrium and Stability of Stellarators
152 ! We break up the Omega-subs into a positive shear term (Dshear) and a net current term, Dcurr
153 ! Omega_subw == Dwell and Omega-subd == Dgeod (geodesic curvature, Pfirsch-Schluter term)
154 !
155 ! Include (eventually) Suydam for reference (cylindrical limit)
156 !
157 
158  WRITE(nmerc,90)
159  90 FORMAT(6x,'S',10x,'PHI',9x,'IOTA',8x,'SHEAR',7x,' VP ',8x,'WELL',
160  1 8x,'ITOR',7x,'ITOR''',7x,'PRES',7x,'PRES''',/,120('-'))
161 
162  DO i = 2,ns1
163  sqs = p5*(vp_real(i) + vp_real(i+1))*sign_jac
164  IF (sqs .eq. zero) cycle
165  WRITE(nmerc,100) sj(i), hs*sum(phip_real(2:i)),
166  1 p5*(iotas(i+1)+iotas(i)), shear(i)/sqs,
167  2 sqs, -vpp(i)*sign_jac,
168  3 p5*(torcur(i) + torcur(i+1)), ip(i)/sqs,
169  4 p5*(pres(i) + pres(i+1)), presp(i)/sqs
170  END DO
171 
172  100 FORMAT(1p,10e12.4)
173 
174  WRITE(nmerc,190)
175  190 FORMAT(/,6x,'S',8x,'DMerc',8x,'DShear',7x,'DCurr',7x,'DWell',
176  1 7x,'Dgeod',/,100('-'))
177 
178  DO i = 2,ns1
179  tpp(i) = (twopi*twopi)*tpp(i)
180  tjb(i) = (twopi*twopi)*tjb(i)
181  tbb(i) = (twopi*twopi)*tbb(i)
182  tjj(i) = (twopi*twopi)*tjj(i)
183  dshear(i) = shear(i) * shear(i)/4
184  dcurr(i) =-shear(i) * (tjb(i) - ip(i) *tbb(i))
185  dwell(i) = presp(i) * (vpp(i) - presp(i) *tpp(i))*tbb(i)
186  dgeod(i) = tjb(i) *tjb(i) - tbb(i) *tjj(i)
187  dmerc(i) = dshear(i) + dcurr(i) + dwell(i) + dgeod(i)
188  WRITE(nmerc,100) sj(i), dmerc(i), dshear(i),
189  1 dcurr(i), dwell(i), dgeod(i)
190  END DO
191 
192  CLOSE (nmerc)
193 
194  END SUBROUTINE mercier
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17