1 SUBROUTINE mercier(gsqrt, bsq, bdotj, iotas, wint,
2 1 r1, rt, rz, zt, zz, bsubu, vp, phips, pres, ns, nznt)
5 USE vmec_input,
ONLY: input_extension
6 USE vparams,
ONLY: one, zero, twopi, nmercier0
11 INTEGER,
INTENT(in) :: ns, nznt
12 REAL(rprec),
DIMENSION(ns,nznt),
INTENT(in) ::
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) ::
18 REAL(rprec),
DIMENSION(ns),
INTENT(in) ::
19 1 iotas, vp, phips, pres
23 REAL(rprec),
PARAMETER :: p5 = 0.5_dp, two = 2
27 INTEGER :: ns1, i, imercier0, nmerc = nmercier0, nrzt
28 REAL(rprec) :: sign_jac, hs, sqs, denom
29 REAL(rprec),
DIMENSION(:,:),
ALLOCATABLE ::
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
38 mercier_file =
'mercier.'//trim(input_extension)
39 CALL safe_open (nmerc, imercier0, mercier_file,
'replace',
41 IF (imercier0 .ne. 0)
RETURN
44 ALLOCATE (gpp(ns,nznt), gsqrt_full(ns,nznt), b2(ns,nznt), stat=i)
45 IF (i .ne. 0) stop
'allocation error in Mercier'
73 IF (ns1 .le. 0)
RETURN
76 IF (gsqrt(ns,1) .ne. zero)
77 1 sign_jac = abs(gsqrt(ns,1))/gsqrt(ns,1)
79 IF (sign_jac .eq. zero)
RETURN
80 phip_real = twopi * phips * sign_jac
86 vp_real(2:ns) = sign_jac*(twopi*twopi)*vp(2:ns)/phip_real(2:ns)
92 torcur(i)=sign_jac*twopi*sum(bsubu(i:nrzt:ns)*wint(i:nrzt:ns))
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
102 vpp(i) = (vp_real(i+1) - vp_real(i))*denom
103 presp(i) = (pres(i+1) - pres(i))*denom
104 ip(i) = (torcur(i+1) - torcur(i))*denom
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)
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)
133 b2(i,:) = two*(bsq(i,:) - pres(i))
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))
148 DEALLOCATE (gpp, gsqrt_full, b2, stat=i)
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(
'-'))
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
172 100
FORMAT(1p,10e12.4)
175 190
FORMAT(/,6x,
'S',8x,
'DMerc',8x,
'DShear',7x,
'DCurr',7x,
'DWell',
176 1 7x,
'Dgeod',/,100(
'-'))
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)
194 END SUBROUTINE mercier