V3FIT
virtual_casing_mod.f90
1 !-----------------------------------------------------------------------
2 ! Module: virtual_casing_mod
3 ! Authors: S. Lazerson (lazerson@pppl.gov)
4 ! Date: 03/27/2012
5 ! Description: This subroutine calculates the plasma response in
6 ! real space coordinates through a virtual casing
7 ! principle. The code utilizes the ezspline package
8 ! which is part of the NTCC. Note if adapt_tol is
9 ! set less than zero the adaptive call will call the
10 ! non-adaptive routines.
11 ! Usage:
12 ! To load an equilibrium using virtual casing from VMEC
13 !
14 ! ALLOCATE(xm_temp(mnmax),xn_temp(mnmax)) ! Integer Arrays
15 ! ALLOCATE(rmnc_temp(mnmax,2),zmns_temp(mnmax,2)) ! DOUBLE PRECISION
16 ! ALLOCATE(bumnc_temp(mnmax,1),bvmnc_temp(mnmax,1)) ! DOUBLE PRECISION
17 ! xm_temp=xm
18 ! xn_temp=-xn/nfp
19 ! rmnc_temp(:,1)=rmnc(:,ns-1)
20 ! rmnc_temp(:,2)=rmnc(:,ns)
21 ! zmns_temp(:,1)=zmns(:,ns-1)
22 ! zmns_temp(:,2)=zmns(:,ns)
23 ! bumnc_temp(:,1) = (1.5*bsupumnc(:,ns) - 0.5*bsupumnc(:,ns-1))
24 ! bvmnc_temp(:,1) = (1.5*bsupvmnc(:,ns) - 0.5*bsupvmnc(:,ns-1))
25 ! CALL init_virtual_casing(mnmax,nu2,nv2,xm_temp,xn_temp,&
26 ! rmnc_temp,zmns_temp,nfp,&
27 ! BUMNC=bumnc_temp,BVMNC=bvmnc_temp)
28 ! -or-
29 ! CALL init_virtual_casing(mnmax,nu2,nv2,xm_temp,xn_temp,&
30 ! rmnc_temp,zmns_temp,nfp,&
31 ! RMNS=rmns_temp, ZMNC=zmnc_temp,&
32 ! BUMNC=bumnc_temp,BVMNC=bvmnc_temp,&
33 ! BUMNS=bumns_temp,BVMNS=bvmns_temp)
34 !
35 ! To load an equilibrium using current density from VMEC
36 !
37 ! ALLOCATE(xm_temp(mnmax),xn_temp(mnmax)) ! Integer Arrays
38 ! ALLOCATE(rmnc_temp(mnmax,ns),zmns_temp(mnmax,ns)) ! DOUBLE PRECISION
39 ! ALLOCATE(jumnc_temp(mnmax,ns),jvmnc_temp(mnmax,ns)) ! DOUBLE PRECISION
40 ! xm_temp = xm
41 ! xn_temp = -xn
42 ! rmnc_temp = rmnc
43 ! zmns_temp = zmns
44 ! jumnc_temp = isigng*currumnc
45 ! jvmnc_temp = isigng*currvmnc
46 ! CALL init_volint(mnmax,nu2,nv2,ns,xm_temp,xn_temp, &
47 ! rmnc_temp,zmns_temp,nfp,&
48 ! JUMNC=jumnc_temp, JVMNC=jvmnc_temp)
49 !
50 ! To evaluate the field at a point in space call
51 ! CALL bfield_vc(xp,yp,zp,bxp,byp,bzp,ier)
52 ! -or-
53 ! CALL vecpot_vc(xp,yp,zp,ax,ay,az,ier)
54 !
55 ! Note that setting adapt_tol < 0 will cause the code
56 ! to use non-adaptive integration routines. Also note
57 ! that if the code is compiled without the NTCC PSLPLINE
58 ! routines the calls will default to non-adaptive integration.
59 !
60 !
61 !
62 ! Notes: 6/15/12 - SAL
63 ! Reorganized the loops for left-hand memory
64 ! efficiency. Realspace matricies are now stored in
65 ! [nu,nv,k] order instead of [k,nu,nv]
66 ! 4/29/13 - SAL
67 ! Set minimum number of calls to equal to 0.
68 ! 7/29/13 - SAL
69 ! Minimum number of calls may be set by calling codes
70 ! 7/30/13 - SAL
71 ! Adjusted screen output to reflect non-adaptive
72 ! routine selection.
73 ! Caught bug in init_virtual_casing_dbl where part of
74 ! the nx, ny, nz arrays weren't being calculated
75 ! properly.
76 ! Added IWRK Parameter to control NAG helper array
77 ! size.
78 ! 1/9/14 - SAL
79 ! Benchmarking of code for various configurations
80 ! resulted in improved code character. Surface
81 ! normals benchmarked against VMEC ouput (surface
82 ! area). B-Field benchmarked against VMEC. These
83 ! values now appear to be correct. Although small
84 ! errors in B for stellarators and big errors in A
85 ! persist. Benchmarking on flux loops indicates
86 ! flux correctly calculated despite A discrepancy.
87 ! 9/10/14 - SAL
88 ! Fixed local deallocations.
89 ! 12/17/15 - SAL
90 ! Changed INTEGER to INTEGER(KIND=8) for gfortran and nag
91 ! Change deffinitino of norm_nag to remove factor of nfp
92 !-----------------------------------------------------------------------
93  MODULE virtual_casing_mod
94 !-----------------------------------------------------------------------
95 ! Libraries
96 !-----------------------------------------------------------------------
97  USE ezspline_obj
98  USE ezspline
99 !-----------------------------------------------------------------------
100 ! Module Variables
101 !-----------------------------------------------------------------------
102  IMPLICIT NONE
103  INTEGER :: nu_vc,nv_vc,nvp,nuvp,nr_vc,nextcur_vc,&
104  nlastcall
105  INTEGER :: MIN_CLS = 0
106  DOUBLE PRECISION :: norm, min_delta_x, x_nag, y_nag, z_nag, &
107  adapt_tol, norm_nag, norm_3d, adapt_rel, &
108  surf_area
109  DOUBLE PRECISION, ALLOCATABLE :: xsurf(:),ysurf(:),zsurf(:)
110  DOUBLE PRECISION, ALLOCATABLE :: nxsurf(:),nysurf(:),nzsurf(:)
111  DOUBLE PRECISION, ALLOCATABLE :: btopx(:),btopy(:),btopz(:),btops(:)
112  DOUBLE PRECISION, ALLOCATABLE :: jx_3d(:), jy_3d(:),jz_3d(:)
113  DOUBLE PRECISION, PARAMETER :: pi2 = 6.283185482025146d+00
114 ! INTEGER, PARAMETER :: IWRK = 67108864 ! 2^26
115  INTEGER, PARAMETER :: IWRK = 16777216 ! 2^24
116  CHARACTER(LEN=256) :: vc_type_str=''
117  DOUBLE PRECISION, PRIVATE, PARAMETER :: zero = 0.0d+0
118  DOUBLE PRECISION, PRIVATE, PARAMETER :: one = 1.0d+0
119 
120  TYPE(EZspline2_r8) :: x_spl, y_spl, z_spl
121  TYPE(EZspline2_r8) :: nx_spl, ny_spl, nz_spl
122  TYPE(EZspline2_r8) :: kx_spl, ky_spl, kz_spl, bn_spl
123  TYPE(EZspline3_r8) :: x3d_spl, y3d_spl, z3d_spl, jx3d_spl, jy3d_spl, jz3d_spl
124 
125 !-----------------------------------------------------------------------
126 ! Module Subroutines
127 ! init_virtual_casing Initializes the atop and surf arrays
128 ! init_virtual_casing_realspace Initizliaes the VC algorithm with realspace values
129 ! virtual_casing_info Outputs Information to the Screen
130 ! virtual_casing_dist Returns distance to surface
131 ! free_virtual_casing Frees Allocated Arrays
132 ! bfield_virtual_casing Calculates the field
133 ! vecpot_virtual_casing Calculates the Vector Potential (A)
134 ! bfield_virtual_casing_adapt Calculates the field using NAG adaptive integration
135 ! vecpot_virtual_casing_adapt Calculates the Vector Potential using NAG adaptive integration (A)
136 !-----------------------------------------------------------------------
137  INTERFACE init_virtual_casing
138  MODULE PROCEDURE init_virtual_casing_flt, init_virtual_casing_dbl
139  END INTERFACE
140  INTERFACE init_virtual_casing_realspace
141  MODULE PROCEDURE init_virtual_casing_realspace_flt, init_virtual_casing_realspace_dbl
142  END INTERFACE
143  INTERFACE init_volint
144  MODULE PROCEDURE init_volint_flt, init_volint_dbl
145  END INTERFACE
146  INTERFACE bfield_virtual_casing
147  MODULE PROCEDURE bfield_virtual_casing_flt, bfield_virtual_casing_dbl
148  END INTERFACE
149  INTERFACE bfield_virtual_casing_adapt
150  MODULE PROCEDURE bfield_virtual_casing_adapt_flt, bfield_virtual_casing_adapt_dbl
151  END INTERFACE
152  INTERFACE vecpot_virtual_casing
153  MODULE PROCEDURE vecpot_virtual_casing_flt, vecpot_virtual_casing_dbl
154  END INTERFACE
155  INTERFACE vecpot_virtual_casing_adapt
156  MODULE PROCEDURE vecpot_virtual_casing_adapt_flt, vecpot_virtual_casing_adapt_dbl
157  END INTERFACE
158  INTERFACE virtual_casing_dist
159  MODULE PROCEDURE virtual_casing_dist_flt, virtual_casing_dist_dbl
160  END INTERFACE
161  INTERFACE bfield_volint_adapt
162  MODULE PROCEDURE bfield_volint_adapt_flt, bfield_volint_adapt_dbl
163  END INTERFACE
164  INTERFACE vecpot_volint_adapt
165  MODULE PROCEDURE vecpot_volint_adapt_flt, vecpot_volint_adapt_dbl
166  END INTERFACE
167  INTERFACE bfield_vc
168  MODULE PROCEDURE bfield_vc_flt, bfield_vc_dbl
169  END INTERFACE
170  INTERFACE vecpot_vc
171  MODULE PROCEDURE vecpot_vc_flt, vecpot_vc_dbl
172  END INTERFACE
173  CONTAINS
174 
175  !-----------------------------------------------------------------
176  SUBROUTINE init_virtual_casing_dbl(mnmax,nu,nv,xm,xn,rmnc,zmns,nfp,bumnc,bvmnc,&
177  rmns,zmnc,bsmns,bsmnc,bumns,bvmns,dr)
178  IMPLICIT NONE
179  ! INPUT VARIABLES
180  INTEGER, INTENT(in) :: mnmax, nu, nv, nfp
181  INTEGER, INTENT(in) :: xm(1:mnmax),xn(1:mnmax)
182  DOUBLE PRECISION, INTENT(in) :: rmnc(1:mnmax,2),zmns(1:mnmax,2)
183  DOUBLE PRECISION, INTENT(in), OPTIONAL :: bumnc(1:mnmax,1),bvmnc(1:mnmax,1)
184  DOUBLE PRECISION, INTENT(in), OPTIONAL :: bsmns(1:mnmax,1),bsmnc(1:mnmax,1)
185  DOUBLE PRECISION, INTENT(in), OPTIONAL :: bumns(1:mnmax,1),bvmns(1:mnmax,1)
186  DOUBLE PRECISION, INTENT(in), OPTIONAL :: rmns(1:mnmax,2),zmnc(1:mnmax,2)
187  DOUBLE PRECISION, INTENT(in), OPTIONAL :: dr
188 
189  ! LOCAL VARIABLES
190  INTEGER :: nuv, mn, uv, i, u, v, dex1, dex2, ier, nuvm
191  INTEGER :: bcs1(2), bcs2(2)
192  DOUBLE PRECISION :: snr, snphi,snz,snx,sny,brs,bphis,bzs,bxs,bys,&
193  factor, cop, sip, dx, dy, dz, sn, signs,xt,yt,zt, dr_temp,&
194  br_vac,bphi_vac,bz_vac, xu_temp, xv_temp, yu, yv, stotal
195  DOUBLE PRECISION, ALLOCATABLE :: xu(:), xv(:)
196  DOUBLE PRECISION, ALLOCATABLE :: fmn_temp(:,:)
197  DOUBLE PRECISION, ALLOCATABLE :: xreal(:,:),yreal(:,:),zreal(:,:)
198  DOUBLE PRECISION, ALLOCATABLE :: nxreal(:,:),nyreal(:,:),nzreal(:,:)
199  DOUBLE PRECISION, ALLOCATABLE :: kxreal(:,:),kyreal(:,:),kzreal(:,:), btopreal(:,:)
200  DOUBLE PRECISION, ALLOCATABLE :: r_temp(:,:,:), z_temp(:,:,:)
201  DOUBLE PRECISION, ALLOCATABLE :: bs_temp(:,:,:), bu_temp(:,:,:), bv_temp(:,:,:)
202  DOUBLE PRECISION, ALLOCATABLE :: rs(:,:,:), ru(:,:,:), rv(:,:,:)
203  DOUBLE PRECISION, ALLOCATABLE :: zs(:,:,:), zu(:,:,:), zv(:,:,:)
204  ! BEGIN SUBROUTINE
205  ! Deallocate anything globally allocated
206  IF (ALLOCATED(xsurf)) DEALLOCATE(xsurf)
207  IF (ALLOCATED(ysurf)) DEALLOCATE(ysurf)
208  IF (ALLOCATED(zsurf)) DEALLOCATE(zsurf)
209  IF (ALLOCATED(nxsurf)) DEALLOCATE(nxsurf)
210  IF (ALLOCATED(nysurf)) DEALLOCATE(nysurf)
211  IF (ALLOCATED(nzsurf)) DEALLOCATE(nzsurf)
212  IF (ALLOCATED(btopx)) DEALLOCATE(btopx)
213  IF (ALLOCATED(btopy)) DEALLOCATE(btopy)
214  IF (ALLOCATED(btopz)) DEALLOCATE(btopz)
215  IF (ALLOCATED(btops)) DEALLOCATE(btops)
216  IF (ezspline_allocated(x_spl)) CALL ezspline_free(x_spl,ier)
217  IF (ezspline_allocated(y_spl)) CALL ezspline_free(y_spl,ier)
218  IF (ezspline_allocated(z_spl)) CALL ezspline_free(z_spl,ier)
219  IF (ezspline_allocated(nx_spl)) CALL ezspline_free(nx_spl,ier)
220  IF (ezspline_allocated(ny_spl)) CALL ezspline_free(ny_spl,ier)
221  IF (ezspline_allocated(nz_spl)) CALL ezspline_free(nz_spl,ier)
222  IF (ezspline_allocated(kx_spl)) CALL ezspline_free(kx_spl,ier)
223  IF (ezspline_allocated(ky_spl)) CALL ezspline_free(ky_spl,ier)
224  IF (ezspline_allocated(kz_spl)) CALL ezspline_free(kz_spl,ier)
225  IF (ezspline_allocated(bn_spl)) CALL ezspline_free(bn_spl,ier)
226 
227  ! Initialize varaibles
228  bcs1=(/-1,-1/)
229  bcs2=(/-1,-1/)
230  adapt_tol = 1.0d-5
231  adapt_rel = 1.0d-3
232  factor = pi2 / dble(nfp)
233  nr_vc = 1
234  nu_vc = nu
235  nv_vc = nv
236  nuv = nu * nv
237  nvp = nv * nfp
238  nuvp = nu * nv * nfp
239  norm = one/(2*pi2*nuvp)
240  norm_nag = one/(pi2*2)
241  vc_type_str='Surface Current'
242  ! Alloocations
243  ALLOCATE(xu(1:nu),xv(1:nv))
244  ALLOCATE(xsurf(1:nuvp),ysurf(1:nuvp),zsurf(1:nuvp))
245  ALLOCATE(nxsurf(1:nuvp),nysurf(1:nuvp),nzsurf(1:nuvp))
246  ALLOCATE(r_temp(1:nu,1:nv,2),z_temp(1:nu,1:nv,2))
247  ALLOCATE(bs_temp(1:nu,1:nv,1),bu_temp(1:nu,1:nv,1),bv_temp(1:nu,1:nv,1))
248  ALLOCATE(rs(1:nu,1:nv,1),ru(1:nu,1:nv,1),rv(1:nu,1:nv,1))
249  ALLOCATE(zs(1:nu,1:nv,1),zu(1:nu,1:nv,1),zv(1:nu,1:nv,1))
250  ALLOCATE(xreal(1:nu+1,1:nvp+1),yreal(1:nu+1,1:nvp+1),zreal(1:nu+1,1:nvp+1))
251  ALLOCATE(nxreal(1:nu+1,1:nvp+1),nyreal(1:nu+1,1:nvp+1),nzreal(1:nu+1,1:nvp+1))
252  ALLOCATE(kxreal(1:nu+1,1:nvp+1),kyreal(1:nu+1,1:nvp+1),kzreal(1:nu+1,1:nvp+1))
253  ALLOCATE(btopreal(1:nu+1,1:nvp+1))
254  r_temp=zero; z_temp=zero; bs_temp=zero; bu_temp=zero; bv_temp=zero
255  dr_temp=one;
256  FORALL(u=1:nu) xu(u) = dble(u-1)/dble(nu)
257  FORALL(v=1:nv) xv(v) = dble(v-1)/dble(nv)
258  CALL mntouv_local(1,2,mnmax,nu,nv,xu,xv,rmnc,xm,xn,r_temp,0,1)
259  CALL mntouv_local(1,2,mnmax,nu,nv,xu,xv,zmns,xm,xn,z_temp,1,0)
260  IF (PRESENT(bsmns)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bsmns,xm,xn,bs_temp,1,0)
261  IF (PRESENT(bumnc)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bumnc,xm,xn,bu_temp,0,0)
262  IF (PRESENT(bvmnc)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bvmnc,xm,xn,bv_temp,0,0)
263  IF (PRESENT(rmns)) CALL mntouv_local(1,2,mnmax,nu,nv,xu,xv,rmns,xm,xn,r_temp,1,0)
264  IF (PRESENT(zmnc)) CALL mntouv_local(1,2,mnmax,nu,nv,xu,xv,zmnc,xm,xn,z_temp,0,0)
265  IF (PRESENT(bsmnc)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bsmnc,xm,xn,bs_temp,0,0)
266  IF (PRESENT(bumns)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bumns,xm,xn,bu_temp,1,0)
267  IF (PRESENT(bvmns)) CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,bvmns,xm,xn,bv_temp,1,0)
268  IF (PRESENT(dr)) dr_temp = dr
269 
270  xreal = zero; yreal = zero; zreal=zero
271  uv = 1
272  DO v = 1, nv
273  DO u = 1, nu
274  xsurf(uv) = r_temp(u,v,2)*dcos(factor*xv(v))
275  ysurf(uv) = r_temp(u,v,2)*dsin(factor*xv(v))
276  zsurf(uv) = z_temp(u,v,2)
277  xreal(u,v) = xsurf(uv)
278  yreal(u,v) = ysurf(uv)
279  zreal(u,v) = zsurf(uv)
280  uv = uv + 1
281  END DO
282  END DO
283  ! Note we need to extend to nu+1
284  xreal(nu+1,:) = xreal(1,:)
285  yreal(nu+1,:) = yreal(1,:)
286  zreal(nu+1,:) = zreal(1,:)
287 
288  ! Now we calculate the edge metric elements
289  ALLOCATE(fmn_temp(1:mnmax,1))
290  rs = zero; zs = 0.0; ru = zero; zu = zero; rv = zero; zv = zero
291  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = -rmnc(mn,2)*xm(mn)
292  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,ru,1,0)
293  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = -rmnc(mn,2)*xn(mn)
294  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,rv,1,0)
295  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = zmns(mn,2)*xm(mn)
296  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,zu,0,0)
297  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = zmns(mn,2)*xn(mn)
298  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,zv,0,0)
299  IF (PRESENT(rmns)) THEN
300  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = rmns(mn,2)*xm(mn)
301  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,ru,0,0)
302  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = rmns(mn,2)*xn(mn)
303  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,rv,0,0)
304  END IF
305  IF (PRESENT(zmnc)) THEN
306  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = -zmnc(mn,2)*xm(mn)
307  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,zu,1,0)
308  FORALL(mn = 1:mnmax) fmn_temp(mn,1) = -zmnc(mn,2)*xn(mn)
309  CALL mntouv_local(1,1,mnmax,nu,nv,xu,xv,fmn_temp,xm,xn,zv,1,0)
310  END IF
311  DO u=1,nu
312  DO v=1,nv
313  rs(u,v,1) = (r_temp(u,v,2) - r_temp(u,v,1))
314  zs(u,v,1) = (z_temp(u,v,2) - z_temp(u,v,1))
315  END DO
316  END DO
317 
318  ALLOCATE(btopx(1:nuvp), btopy(1:nuvp), btopz(1:nuvp),btops(1:nuvp))
319  btopx = zero; btopy = zero; btopz = zero; btops = zero; btopreal = zero
320  kxreal = zero; kyreal = zero; kzreal = zero;
321  nxreal = zero; nyreal = zero; nzreal = zero;
322  uv = 1
323  ! Check sign of jacobian
324  signs = one
325  stotal = zero
326  IF (zreal(2,1)-zreal(1,1) < 0) signs = -one
327  ! u then v if not outputting
328  DO v = 1, nv
329  DO u = 1, nu
330  cop = dcos(factor*xv(v))
331  sip = dsin(factor*xv(v))
332  xu_temp = pi2*ru(u,v,1)*cop
333  yu = pi2*ru(u,v,1)*sip
334  xv_temp = pi2*rv(u,v,1)*cop - factor*yreal(u,v)
335  yv = pi2*rv(u,v,1)*sip + factor*xreal(u,v)
336  snx = -yu*zv(u,v,1)*pi2 + zu(u,v,1)*yv*pi2
337  sny = -xv_temp*zu(u,v,1)*pi2 + zv(u,v,1)*xu_temp*pi2
338  snz = -xu_temp*yv + yu*xv_temp
339  sn = dsqrt(snx*snx+sny*sny+snz*snz)
340  stotal = stotal + sn
341  brs = bs_temp(u,v,1)*rs(u,v,1)+bu_temp(u,v,1)*ru(u,v,1)+bv_temp(u,v,1)*rv(u,v,1)*nfp
342  bphis = bv_temp(u,v,1)*r_temp(u,v,2)
343  bzs = bs_temp(u,v,1)*zs(u,v,1)+bu_temp(u,v,1)*zu(u,v,1)+bv_temp(u,v,1)*zv(u,v,1)*nfp
344  bxs = brs*cop - bphis*sip
345  bys = brs*sip + bphis*cop
346  btopx(uv) = -( bys * snz - bzs * sny )
347  btopy(uv) = -( bzs * snx - bxs * snz )
348  btopz(uv) = -( bxs * sny - bys * snx )
349  kxreal(u,v) = btopx(uv)
350  kyreal(u,v) = btopy(uv)
351  kzreal(u,v) = btopz(uv)
352  nxsurf(uv) = snx/sn
353  nysurf(uv) = sny/sn
354  nzsurf(uv) = snz/sn
355  nxreal(u,v) = snx/sn
356  nyreal(u,v) = sny/sn
357  nzreal(u,v) = snz/sn
358  uv = uv + 1
359  END DO
360  END DO
361  surf_area = nfp*stotal/nuv
362  CALL flush(6)
363  ! NUMAX
364  kxreal(nu+1,:) = kxreal(1,:)
365  kyreal(nu+1,:) = kyreal(1,:)
366  kzreal(nu+1,:) = kzreal(1,:)
367  nxreal(nu+1,:) = nxreal(1,:)
368  nyreal(nu+1,:) = nyreal(1,:)
369  nzreal(nu+1,:) = nzreal(1,:)
370  btopreal(nu+1,:) = btopreal(1,:)
371  ! NVMAX
372  kxreal(:,nv+1) = kxreal(:,1)
373  kyreal(:,nv+1) = kyreal(:,1)
374  kzreal(:,nv+1) = kzreal(:,1)
375  nxreal(:,nv+1) = nxreal(:,1)
376  nyreal(:,nv+1) = nyreal(:,1)
377  nzreal(:,nv+1) = nzreal(:,1)
378  btopreal(:,nv+1) = btopreal(:,1)
379  ! Last point
380  kxreal(nu+1,nv+1) = kxreal(1,1)
381  kyreal(nu+1,nv+1) = kyreal(1,1)
382  kzreal(nu+1,nv+1) = kzreal(1,1)
383  nxreal(nu+1,nv+1) = nxreal(1,1)
384  nyreal(nu+1,nv+1) = nyreal(1,1)
385  nzreal(nu+1,nv+1) = nzreal(1,1)
386  btopreal(nu+1,nv+1) = btopreal(nu+1,nv+1)
387  ! Finish off
388  DO v = 2, nfp
389  cop = dcos((v-1)*factor)
390  sip = dsin((v-1)*factor)
391  dex1 = (v-1)*nuv+1
392  dex2 = v*nuv
393  btopx(dex1:dex2) = btopx(1:nuv)*cop - btopy(1:nuv)*sip
394  btopy(dex1:dex2) = btopy(1:nuv)*cop + btopx(1:nuv)*sip
395  btopz(dex1:dex2) = btopz(1:nuv)
396  btops(dex1:dex2) = btops(1:nuv)
397  xsurf(dex1:dex2) = xsurf(1:nuv)*cop - ysurf(1:nuv)*sip
398  ysurf(dex1:dex2) = ysurf(1:nuv)*cop + xsurf(1:nuv)*sip
399  zsurf(dex1:dex2) = zsurf(1:nuv)
400  nxsurf(dex1:dex2) = nxsurf(1:nuv)*cop - nysurf(1:nuv)*sip
401  nysurf(dex1:dex2) = nysurf(1:nuv)*cop + nxsurf(1:nuv)*sip
402  nzsurf(dex1:dex2) = nzsurf(1:nuv)
403  dex1 = (v-1)*nv+1
404  dex2 = v*nv
405  xreal(1:nu+1,dex1:dex2) = xreal(1:nu+1,1:nv)*cop - yreal(1:nu+1,1:nv)*sip
406  yreal(1:nu+1,dex1:dex2) = yreal(1:nu+1,1:nv)*cop + xreal(1:nu+1,1:nv)*sip
407  zreal(1:nu+1,dex1:dex2) = zreal(1:nu+1,1:nv)
408  kxreal(1:nu+1,dex1:dex2) = kxreal(1:nu+1,1:nv)*cop - kyreal(1:nu+1,1:nv)*sip
409  kyreal(1:nu+1,dex1:dex2) = kyreal(1:nu+1,1:nv)*cop + kxreal(1:nu+1,1:nv)*sip
410  kzreal(1:nu+1,dex1:dex2) = kzreal(1:nu+1,1:nv)
411  nxreal(1:nu+1,dex1:dex2) = nxreal(1:nu+1,1:nv)*cop - nyreal(1:nu+1,1:nv)*sip
412  nyreal(1:nu+1,dex1:dex2) = nyreal(1:nu+1,1:nv)*cop + nxreal(1:nu+1,1:nv)*sip
413  nzreal(1:nu+1,dex1:dex2) = nzreal(1:nu+1,1:nv)
414  btopreal(1:nu+1,dex1:dex2) = btopreal(1:nu+1,1:nv)
415  END DO
416  xreal(:,nvp+1) = xreal(:,1)
417  yreal(:,nvp+1) = yreal(:,1)
418  zreal(:,nvp+1) = zreal(:,1)
419  kxreal(:,nvp+1) = kxreal(:,1)
420  kyreal(:,nvp+1) = kyreal(:,1)
421  kzreal(:,nvp+1) = kzreal(:,1)
422  nxreal(:,nvp+1) = nxreal(:,1)
423  nyreal(:,nvp+1) = nyreal(:,1)
424  nzreal(:,nvp+1) = nzreal(:,1)
425  btopreal(:,nvp+1) = btopreal(:,1)
426  ! Calculate Return map
427  min_delta_x = 1.0d+10
428  DO v = 2, nv-1
429  DO u = 2, nu-1
430  dx = xreal(nu+1,nv) - xreal(nu-1,nv)
431  dy = yreal(nu+1,nv) - yreal(nu-1,nv)
432  dz = zreal(nu+1,nv) - zreal(nu-1,nv)
433  min_delta_x=min(min_delta_x,sqrt(dx*dx+dy*dy+dz*dz))
434  dx = xreal(nu,nv+1) - xreal(nu,nv-1)
435  dy = yreal(nu,nv+1) - yreal(nu,nv-1)
436  dz = zreal(nu,nv+1) - zreal(nu,nv-1)
437  min_delta_x=min(min_delta_x,sqrt(dx*dx+dy*dy+dz*dz))
438  END DO
439  END DO
440  ! Construct Splines
441  CALL ezspline_init(x_spl,nu+1,nvp+1,bcs1,bcs2,ier)
442  CALL ezspline_init(y_spl,nu+1,nvp+1,bcs1,bcs2,ier)
443  CALL ezspline_init(z_spl,nu+1,nvp+1,bcs1,bcs2,ier)
444  CALL ezspline_init(nx_spl,nu+1,nvp+1,bcs1,bcs2,ier)
445  CALL ezspline_init(ny_spl,nu+1,nvp+1,bcs1,bcs2,ier)
446  CALL ezspline_init(nz_spl,nu+1,nvp+1,bcs1,bcs2,ier)
447  CALL ezspline_init(kx_spl,nu+1,nvp+1,bcs1,bcs2,ier)
448  CALL ezspline_init(ky_spl,nu+1,nvp+1,bcs1,bcs2,ier)
449  CALL ezspline_init(kz_spl,nu+1,nvp+1,bcs1,bcs2,ier)
450  CALL ezspline_init(bn_spl,nu+1,nvp+1,bcs1,bcs2,ier)
451  ! Define arrays from 0 to 1 because integrand already contains
452  ! dA information.
453  DO u = 1, nu+1
454  x_spl%x1(u) = dble(u-1)/dble(nu)
455  y_spl%x1(u) = dble(u-1)/dble(nu)
456  z_spl%x1(u) = dble(u-1)/dble(nu)
457  nx_spl%x1(u) = dble(u-1)/dble(nu)
458  ny_spl%x1(u) = dble(u-1)/dble(nu)
459  nz_spl%x1(u) = dble(u-1)/dble(nu)
460  kx_spl%x1(u) = dble(u-1)/dble(nu)
461  ky_spl%x1(u) = dble(u-1)/dble(nu)
462  kz_spl%x1(u) = dble(u-1)/dble(nu)
463  bn_spl%x1(u) = dble(u-1)/dble(nu)
464  END DO
465  DO v = 1, nvp+1
466  x_spl%x2(v) = dble(v-1)/dble(nvp)
467  y_spl%x2(v) = dble(v-1)/dble(nvp)
468  z_spl%x2(v) = dble(v-1)/dble(nvp)
469  nx_spl%x2(v) = dble(v-1)/dble(nvp)
470  ny_spl%x2(v) = dble(v-1)/dble(nvp)
471  nz_spl%x2(v) = dble(v-1)/dble(nvp)
472  kx_spl%x2(v) = dble(v-1)/dble(nvp)
473  ky_spl%x2(v) = dble(v-1)/dble(nvp)
474  kz_spl%x2(v) = dble(v-1)/dble(nvp)
475  bn_spl%x2(v) = dble(v-1)/dble(nvp)
476  END DO
477  x_spl%isHermite = 1
478  y_spl%isHermite = 1
479  z_spl%isHermite = 1
480  nx_spl%isHermite = 1
481  ny_spl%isHermite = 1
482  nz_spl%isHermite = 1
483  kx_spl%isHermite = 1
484  ky_spl%isHermite = 1
485  kz_spl%isHermite = 1
486  bn_spl%isHermite = 1
487  CALL ezspline_setup(x_spl,xreal,ier)
488  CALL ezspline_setup(y_spl,yreal,ier)
489  CALL ezspline_setup(z_spl,zreal,ier)
490  CALL ezspline_setup(nx_spl,nxreal,ier)
491  CALL ezspline_setup(ny_spl,nyreal,ier)
492  CALL ezspline_setup(nz_spl,nzreal,ier)
493  CALL ezspline_setup(kx_spl,kxreal,ier)
494  CALL ezspline_setup(ky_spl,kyreal,ier)
495  CALL ezspline_setup(kz_spl,kzreal,ier)
496  CALL ezspline_setup(bn_spl,btopreal,ier)
497  DEALLOCATE(xu,xv)
498  DEALLOCATE(r_temp,z_temp)
499  DEALLOCATE(bs_temp,bu_temp,bv_temp)
500  DEALLOCATE(rs,ru,rv)
501  DEALLOCATE(zs,zu,zv)
502  DEALLOCATE(xreal,yreal,zreal)
503  DEALLOCATE(nxreal,nyreal,nzreal)
504  DEALLOCATE(kxreal,kyreal,kzreal,btopreal)
505  DEALLOCATE(fmn_temp)
506  ! END SUBROUTINE
507  END SUBROUTINE init_virtual_casing_dbl
508  !-----------------------------------------------------------------
509 
510  !-----------------------------------------------------------------
511  SUBROUTINE init_virtual_casing_realspace_dbl(nu,nv,nfp,phi,&
512  rreal,zreal,&
513  snr,snphi,snz,&
514  brreal,bphireal,bzreal)
515  IMPLICIT NONE
516  ! INPUT VARIABLES
517  INTEGER, INTENT(in) :: nu, nv, nfp
518  DOUBLE PRECISION, INTENT(in) :: phi(nv)
519  DOUBLE PRECISION, INTENT(in) :: rreal(nu,nv), zreal(nu,nv)
520  DOUBLE PRECISION, INTENT(in) :: snr(nu,nv), snphi(nu,nv), snz(nu,nv)
521  DOUBLE PRECISION, INTENT(in) :: brreal(nu,nv), bphireal(nu,nv), bzreal(nu,nv)
522  ! LOCAL VARIABLES
523  INTEGER :: nuv, mn, uv, i, u, v, dex1, dex2, ier, nuvm
524  INTEGER :: bcs1(2), bcs2(2)
525  DOUBLE PRECISION :: snx,sny,bzs,bxs,bys,&
526  factor, cop, sip, dx, dy, dz, sn, signs,xt,yt,zt, dr_temp
527  DOUBLE PRECISION, ALLOCATABLE :: fmn_temp(:,:)
528  DOUBLE PRECISION, ALLOCATABLE :: xreal(:,:),yreal(:,:),zreal2(:,:)
529  DOUBLE PRECISION, ALLOCATABLE :: nxreal(:,:),nyreal(:,:),nzreal(:,:)
530  DOUBLE PRECISION, ALLOCATABLE :: kxreal(:,:),kyreal(:,:),kzreal(:,:), btopreal(:,:)
531  DOUBLE PRECISION, ALLOCATABLE :: r_temp(:,:,:), z_temp(:,:,:)
532  ! BEGIN SUBROUTINE
533  ! Deallocate anything globally allocated
534  IF (ALLOCATED(xsurf)) DEALLOCATE(xsurf)
535  IF (ALLOCATED(ysurf)) DEALLOCATE(ysurf)
536  IF (ALLOCATED(zsurf)) DEALLOCATE(zsurf)
537  IF (ALLOCATED(nxsurf)) DEALLOCATE(nxsurf)
538  IF (ALLOCATED(nysurf)) DEALLOCATE(nysurf)
539  IF (ALLOCATED(nzsurf)) DEALLOCATE(nzsurf)
540  IF (ALLOCATED(btopx)) DEALLOCATE(btopx)
541  IF (ALLOCATED(btopy)) DEALLOCATE(btopy)
542  IF (ALLOCATED(btopz)) DEALLOCATE(btopz)
543  IF (ALLOCATED(btops)) DEALLOCATE(btops)
544  IF (ezspline_allocated(x_spl)) CALL ezspline_free(x_spl,ier)
545  IF (ezspline_allocated(y_spl)) CALL ezspline_free(y_spl,ier)
546  IF (ezspline_allocated(z_spl)) CALL ezspline_free(z_spl,ier)
547  IF (ezspline_allocated(nx_spl)) CALL ezspline_free(nx_spl,ier)
548  IF (ezspline_allocated(ny_spl)) CALL ezspline_free(ny_spl,ier)
549  IF (ezspline_allocated(nz_spl)) CALL ezspline_free(nz_spl,ier)
550  IF (ezspline_allocated(kx_spl)) CALL ezspline_free(kx_spl,ier)
551  IF (ezspline_allocated(ky_spl)) CALL ezspline_free(ky_spl,ier)
552  IF (ezspline_allocated(kz_spl)) CALL ezspline_free(kz_spl,ier)
553  IF (ezspline_allocated(bn_spl)) CALL ezspline_free(bn_spl,ier)
554  ! Initialize varaibles
555  bcs1=(/-1,-1/)
556  bcs2=(/-1,-1/)
557  adapt_tol = 1.0d-5
558  adapt_rel = 1.0d-3
559  factor = pi2 / dble(nfp)
560  nu_vc = nu
561  nv_vc = nv
562  nuv = nu * nv
563  nvp = nv * nfp
564  nuvp = nu * nv * nfp
565  norm = one/(2*pi2*nu*nv)
566  norm_nag = one/(2*pi2)
567  vc_type_str='Surface Current (REALSPACE)'
568  ! Alloocations
569  ALLOCATE(xsurf(1:nuvp),ysurf(1:nuvp),zsurf(1:nuvp))
570  ALLOCATE(nxsurf(1:nuvp),nysurf(1:nuvp),nzsurf(1:nuvp))
571  ALLOCATE(r_temp(2,1:nu,1:nv),z_temp(2,1:nu,1:nv))
572  ALLOCATE(xreal(1:nu+1,1:nvp+1),yreal(1:nu+1,1:nvp+1),zreal2(1:nu+1,1:nvp+1))
573  ALLOCATE(nxreal(1:nu+1,1:nvp+1),nyreal(1:nu+1,1:nvp+1),nzreal(1:nu+1,1:nvp+1))
574  ALLOCATE(kxreal(1:nu+1,1:nvp+1),kyreal(1:nu+1,1:nvp+1),kzreal(1:nu+1,1:nvp+1))
575  ALLOCATE(btopreal(1:nu+1,1:nvp+1))
576  dr_temp=one;
577  uv = 1
578  DO v = 1, nv
579  DO u = 1, nu
580  xsurf(uv) = rreal(u,v)*dcos(phi(nv))
581  ysurf(uv) = rreal(u,v)*dsin(phi(nv))
582  zsurf(uv) = zreal(u,v)
583  xreal(u,v) = xsurf(uv)
584  yreal(u,v) = ysurf(uv)
585  zreal2(u,v) = zsurf(uv)
586  uv = uv + 1
587  END DO
588  END DO
589  ! Note we need to extend to nu+1
590  xreal(nu+1,:) = xreal(1,:)
591  yreal(nu+1,:) = yreal(1,:)
592  zreal2(nu+1,:) = zreal2(1,:)
593 
594  ALLOCATE(btopx(1:nuvp), btopy(1:nuvp), btopz(1:nuvp),btops(1:nuvp))
595  btopx = zero; btopy = zero; btopz = zero; btops = zero; btopreal = zero
596  uv = 1
597  ! Check sign of jacobian
598  signs = -one
599  IF (zreal2(2,1) < 0) signs = 1.0
600  ! u then v if not outputting
601  DO v = 1, nv
602  DO u = 1, nu
603  cop = dcos(phi(v))
604  sip = dsin(phi(v))
605  snx = snr(u,v)*cop - snphi(u,v)*sip
606  sny = snr(u,v)*sip + snphi(u,v)*cop
607  sn = dsqrt(snx*snx+sny*sny+snz(u,v)*snz(u,v))
608  bxs = brreal(u,v)*cop - bphireal(u,v)*sip
609  bys = brreal(u,v)*sip + bphireal(u,v)*cop
610  bzs = bzreal(u,v)
611  btopx(uv) = signs * ( bys * snz(u,v) - bzs * sny )
612  btopy(uv) = signs * ( bzs * snx - bxs * snz(u,v) )
613  btopz(uv) = signs * ( bxs * sny - bys * snx )
614  kxreal(u,v) = btopx(uv)
615  kyreal(u,v) = btopy(uv)
616  kzreal(u,v) = btopz(uv)
617  nxsurf(uv) = snx/sn
618  nysurf(uv) = sny/sn
619  nzsurf(uv) = snz(u,v)/sn
620  nxreal(u,v) = snx/sn
621  nyreal(u,v) = sny/sn
622  nzreal(u,v) = snz(u,v)/sn
623  btops(uv) = bxs*snx+bys*sny+bzs*snz(u,v)
624  btopreal(u,v) = btops(uv)
625  uv = uv + 1
626  END DO
627  END DO
628  kxreal(nu+1,:) = kxreal(1,:)
629  kyreal(nu+1,:) = kyreal(1,:)
630  kzreal(nu+1,:) = kzreal(1,:)
631  nxreal(nu+1,:) = nxreal(1,:)
632  nyreal(nu+1,:) = nyreal(1,:)
633  nzreal(nu+1,:) = nzreal(1,:)
634  DO v = 2, nfp
635  cop = dcos((v-1)*factor)
636  sip = dsin((v-1)*factor)
637  dex1 = (v-1)*nuv+1
638  dex2 = v*nuv
639  btopx(dex1:dex2) = btopx(1:nuv)*cop - btopy(1:nuv)*sip
640  btopy(dex1:dex2) = btopy(1:nuv)*cop + btopx(1:nuv)*sip
641  btopz(dex1:dex2) = btopz(1:nuv)
642  btops(dex1:dex2) = btops(1:nuv)
643  xsurf(dex1:dex2) = xsurf(1:nuv)*cop - ysurf(1:nuv)*sip
644  ysurf(dex1:dex2) = ysurf(1:nuv)*cop + xsurf(1:nuv)*sip
645  zsurf(dex1:dex2) = zsurf(1:nuv)
646  nxsurf(dex1:dex2) = nxsurf(1:nuv)*cop - nysurf(1:nuv)*sip
647  nysurf(dex1:dex2) = nysurf(1:nuv)*cop + nxsurf(1:nuv)*sip
648  nzsurf(dex1:dex2) = nzsurf(1:nuv)
649  dex1 = (v-1)*nv+1
650  dex2 = v*nv
651  xreal(1:nu+1,dex1:dex2) = xreal(1:nu+1,1:nv)*cop - yreal(1:nu+1,1:nv)*sip
652  yreal(1:nu+1,dex1:dex2) = yreal(1:nu+1,1:nv)*cop + xreal(1:nu+1,1:nv)*sip
653  zreal2(1:nu+1,dex1:dex2) = zreal2(1:nu+1,1:nv)
654  kxreal(1:nu+1,dex1:dex2) = kxreal(1:nu+1,1:nv)*cop - kyreal(1:nu+1,1:nv)*sip
655  kyreal(1:nu+1,dex1:dex2) = kyreal(1:nu+1,1:nv)*cop + kxreal(1:nu+1,1:nv)*sip
656  kzreal(1:nu+1,dex1:dex2) = kzreal(1:nu+1,1:nv)
657  nxreal(1:nu+1,dex1:dex2) = nxreal(1:nu+1,1:nv)*cop - nyreal(1:nu+1,1:nv)*sip
658  nyreal(1:nu+1,dex1:dex2) = nyreal(1:nu+1,1:nv)*cop + nxreal(1:nu+1,1:nv)*sip
659  nzreal(1:nu+1,dex1:dex2) = nzreal(1:nu+1,1:nv)
660  btopreal(1:nu+1,dex1:dex2) = btopreal(1:nu+1,1:nv)
661  END DO
662  xreal(:,nvp+1) = xreal(:,1)
663  yreal(:,nvp+1) = yreal(:,1)
664  zreal2(:,nvp+1) = zreal2(:,1)
665  kxreal(:,nvp+1) = kxreal(:,1)
666  kyreal(:,nvp+1) = kyreal(:,1)
667  kzreal(:,nvp+1) = kzreal(:,1)
668  nxreal(:,nvp+1) = nxreal(:,1)
669  nyreal(:,nvp+1) = nyreal(:,1)
670  nzreal(:,nvp+1) = nzreal(:,1)
671  btopreal(:,nvp+1) = btopreal(:,1)
672  ! Calculate Return map
673  min_delta_x = 1.0d+10
674  DO v = 2, nv-1
675  DO u = 2, nu-1
676  dx = xreal(nu+1,nv) - xreal(nu-1,nv)
677  dy = yreal(nu+1,nv) - yreal(nu-1,nv)
678  dz = zreal2(nu+1,nv) - zreal2(nu-1,nv)
679  min_delta_x=min(min_delta_x,sqrt(dx*dx+dy*dy+dz*dz))
680  dx = xreal(nu,nv+1) - xreal(nu,nv-1)
681  dy = yreal(nu,nv+1) - yreal(nu,nv-1)
682  dz = zreal2(nu,nv+1) - zreal2(nu,nv-1)
683  min_delta_x=min(min_delta_x,sqrt(dx*dx+dy*dy+dz*dz))
684  END DO
685  END DO
686  ! Construct Splines
687  CALL ezspline_init(x_spl,nu+1,nvp+1,bcs1,bcs2,ier)
688  CALL ezspline_init(y_spl,nu+1,nvp+1,bcs1,bcs2,ier)
689  CALL ezspline_init(z_spl,nu+1,nvp+1,bcs1,bcs2,ier)
690  CALL ezspline_init(nx_spl,nu+1,nvp+1,bcs1,bcs2,ier)
691  CALL ezspline_init(ny_spl,nu+1,nvp+1,bcs1,bcs2,ier)
692  CALL ezspline_init(nz_spl,nu+1,nvp+1,bcs1,bcs2,ier)
693  CALL ezspline_init(kx_spl,nu+1,nvp+1,bcs1,bcs2,ier)
694  CALL ezspline_init(ky_spl,nu+1,nvp+1,bcs1,bcs2,ier)
695  CALL ezspline_init(kz_spl,nu+1,nvp+1,bcs1,bcs2,ier)
696  CALL ezspline_init(bn_spl,nu+1,nvp+1,bcs1,bcs2,ier)
697  ! Define arrays from 0 to 1 because integrand already contains
698  ! dA information.
699  DO u = 1, nu+1
700  x_spl%x1(u) = dble(u-1)/dble(nu)
701  y_spl%x1(u) = dble(u-1)/dble(nu)
702  z_spl%x1(u) = dble(u-1)/dble(nu)
703  nx_spl%x1(u) = dble(u-1)/dble(nu)
704  ny_spl%x1(u) = dble(u-1)/dble(nu)
705  nz_spl%x1(u) = dble(u-1)/dble(nu)
706  kx_spl%x1(u) = dble(u-1)/dble(nu)
707  ky_spl%x1(u) = dble(u-1)/dble(nu)
708  kz_spl%x1(u) = dble(u-1)/dble(nu)
709  bn_spl%x1(u) = dble(u-1)/dble(nu)
710  END DO
711  DO v = 1, nvp+1
712  x_spl%x2(v) = dble(v-1)/dble(nvp)
713  y_spl%x2(v) = dble(v-1)/dble(nvp)
714  z_spl%x2(v) = dble(v-1)/dble(nvp)
715  nx_spl%x2(v) = dble(v-1)/dble(nvp)
716  ny_spl%x2(v) = dble(v-1)/dble(nvp)
717  nz_spl%x2(v) = dble(v-1)/dble(nvp)
718  kx_spl%x2(v) = dble(v-1)/dble(nvp)
719  ky_spl%x2(v) = dble(v-1)/dble(nvp)
720  kz_spl%x2(v) = dble(v-1)/dble(nvp)
721  bn_spl%x2(v) = dble(v-1)/dble(nvp)
722  END DO
723  x_spl%isHermite = 1
724  y_spl%isHermite = 1
725  z_spl%isHermite = 1
726  nx_spl%isHermite = 1
727  ny_spl%isHermite = 1
728  nz_spl%isHermite = 1
729  kx_spl%isHermite = 1
730  ky_spl%isHermite = 1
731  kz_spl%isHermite = 1
732  bn_spl%isHermite = 1
733  CALL ezspline_setup(x_spl,xreal,ier)
734  CALL ezspline_setup(y_spl,yreal,ier)
735  CALL ezspline_setup(z_spl,zreal2,ier)
736  CALL ezspline_setup(kx_spl,kxreal,ier)
737  CALL ezspline_setup(ky_spl,kyreal,ier)
738  CALL ezspline_setup(kz_spl,kzreal,ier)
739  CALL ezspline_setup(nx_spl,nxreal,ier)
740  CALL ezspline_setup(ny_spl,nyreal,ier)
741  CALL ezspline_setup(nz_spl,nzreal,ier)
742  CALL ezspline_setup(bn_spl,btopreal,ier)
743 
744  DEALLOCATE(r_temp,z_temp)
745  DEALLOCATE(xreal,yreal,zreal2)
746  DEALLOCATE(kxreal,kyreal,kzreal,btopreal)
747  DEALLOCATE(nxreal,nyreal,nzreal)
748  ! END SUBROUTINE
749  END SUBROUTINE init_virtual_casing_realspace_dbl
750  !-----------------------------------------------------------------
751 
752  !-----------------------------------------------------------------
753  ! Note optional arguments must have a different name so the
754  ! Interface will pick the proper variables. Here: _FLT
755  SUBROUTINE init_virtual_casing_flt(mnmax_flt,nu_flt,nv_flt,xm_flt,&
756  xn_flt,rmnc_flt,zmns_flt,nfp_flt,&
757  bumnc_flt,bvmnc_flt,&
758  rmns_flt,zmnc_flt,bsmns_flt,&
759  bsmnc_flt,bumns_flt,bvmns_flt,&
760  dr_flt)
761  IMPLICIT NONE
762  ! INPUT VARIABLES
763  INTEGER, INTENT(in) :: mnmax_flt, nu_flt, nv_flt, nfp_flt
764  INTEGER, INTENT(in) :: xm_flt(1:mnmax_flt),xn_flt(1:mnmax_flt)
765  REAL, INTENT(in) :: rmnc_flt(1:mnmax_flt,2),zmns_flt(1:mnmax_flt,2)
766  REAL, INTENT(in), OPTIONAL :: bsmns_flt(1:mnmax_flt,1),bsmnc_flt(1:mnmax_flt,1)
767  REAL, INTENT(in), OPTIONAL :: bumnc_flt(1:mnmax_flt,1),bvmnc_flt(1:mnmax_flt,1)
768  REAL, INTENT(in), OPTIONAL :: bumns_flt(1:mnmax_flt,1),bvmns_flt(1:mnmax_flt,1)
769  REAL, INTENT(in), OPTIONAL :: rmns_flt(1:mnmax_flt,2),zmnc_flt(1:mnmax_flt,2)
770  REAL, INTENT(in), OPTIONAL :: dr_flt
771  ! LOCAL VARIABLES
772  DOUBLE PRECISION :: rmnct(1:mnmax_flt,2),zmnst(1:mnmax_flt,2)
773  DOUBLE PRECISION :: rmnst(1:mnmax_flt,2),zmnct(1:mnmax_flt,2)
774  DOUBLE PRECISION :: bumnct(1:mnmax_flt,1),bvmnct(1:mnmax_flt,1)
775  DOUBLE PRECISION :: bsmnst(1:mnmax_flt,1),bsmnct(1:mnmax_flt,1)
776  DOUBLE PRECISION :: bumnst(1:mnmax_flt,1),bvmnst(1:mnmax_flt,1)
777  DOUBLE PRECISION :: drt
778  ! BEGIN SUBROUTINE
779  rmnct = zero; rmnst = zero; zmnct = zero; zmnst = zero
780  bsmnct = zero; bsmnst = zero
781  bumnct = zero; bumnst = zero
782  bvmnct = zero; bvmnst = zero
783  drt = one
784  rmnct = rmnc_flt
785  zmnst = zmns_flt
786  IF (PRESENT(rmns_flt)) rmnst = rmns_flt
787  IF (PRESENT(zmnc_flt)) zmnct = zmnc_flt
788  IF (PRESENT(bsmnc_flt)) bsmnct = bsmnc_flt
789  IF (PRESENT(bsmns_flt)) bsmnst = bsmns_flt
790  IF (PRESENT(bumnc_flt)) bumnct = bumnc_flt
791  IF (PRESENT(bumns_flt)) bumnst = bumns_flt
792  IF (PRESENT(bvmnc_flt)) bvmnct = bvmnc_flt
793  IF (PRESENT(bvmns_flt)) bvmnst = bvmns_flt
794  IF (PRESENT(dr_flt)) drt = dr_flt
795  CALL init_virtual_casing_dbl(mnmax_flt,nu_flt,nv_flt,xm_flt,xn_flt, &
796  rmnct,zmnst,nfp_flt, &
797  rmns=rmnst,zmnc=zmnct, &
798  bsmns=bsmnst, bsmnc=bsmnct, &
799  bumnc=bumnct, bumns=bumnst, &
800  bvmnc=bvmnct, bvmns=bvmnst, &
801  dr=drt)
802  ! END SUBROUTINE
803  END SUBROUTINE init_virtual_casing_flt
804  !-----------------------------------------------------------------
805 
806  !-----------------------------------------------------------------
807  ! Note optional arguments must have a different name so the
808  ! Interface will pick the proper variables. Here: _FLT
809  SUBROUTINE init_virtual_casing_realspace_flt(nu_flt,nv_flt,nfp_flt,phi_flt,&
810  rreal_flt,zreal_flt,&
811  snr_flt,snphi_flt,snz_flt,&
812  brreal_flt,bphireal_flt,bzreal_flt)
813  IMPLICIT NONE
814  ! INPUT VARIABLES
815  INTEGER, INTENT(in) :: nu_flt, nv_flt, nfp_flt
816  REAL, INTENT(in) :: phi_flt(nv_flt)
817  REAL, INTENT(in) :: rreal_flt(nu_flt,nv_flt), zreal_flt(nu_flt,nv_flt)
818  REAL, INTENT(in) :: snr_flt(nu_flt,nv_flt), snphi_flt(nu_flt,nv_flt), snz_flt(nu_flt,nv_flt)
819  REAL, INTENT(in) :: brreal_flt(nu_flt,nv_flt), bphireal_flt(nu_flt,nv_flt), bzreal_flt(nu_flt,nv_flt)
820  ! LOCAL VARIABLES
821  DOUBLE PRECISION :: phi_dbl(nv_flt)
822  DOUBLE PRECISION :: rreal_dbl(nu_flt,nv_flt), zreal_dbl(nu_flt,nv_flt)
823  DOUBLE PRECISION :: snr_dbl(nu_flt,nv_flt), snphi_dbl(nu_flt,nv_flt), snz_dbl(nu_flt,nv_flt)
824  DOUBLE PRECISION :: brreal_dbl(nu_flt,nv_flt), bphireal_dbl(nu_flt,nv_flt), bzreal_dbl(nu_flt,nv_flt)
825  ! BEGIN SUBROUTINE
826  phi_dbl = phi_flt
827  rreal_dbl = rreal_flt
828  zreal_dbl = zreal_flt
829  snr_dbl = snr_flt
830  snphi_dbl = snphi_flt
831  snz_dbl = snz_flt
832  brreal_dbl = brreal_flt
833  bphireal_dbl = bphireal_flt
834  bzreal_dbl = bzreal_flt
835  CALL init_virtual_casing_realspace_dbl(nu_flt,nv_flt,nfp_flt,phi_dbl,&
836  rreal_dbl,zreal_dbl,&
837  snr_dbl,snphi_dbl,snz_dbl,&
838  brreal_dbl,bphireal_dbl,bzreal_dbl)
839  ! END SUBROUTINE
840  END SUBROUTINE init_virtual_casing_realspace_flt
841  !-----------------------------------------------------------------
842 
843 
844  !-----------------------------------------------------------------
845  SUBROUTINE free_virtual_casing
846  IMPLICIT NONE
847  INTEGER :: ier
848  ! BEGIN SUBROUTINE
849  nuvp = 0
850  norm = -one
851  IF (ALLOCATED(xsurf)) DEALLOCATE(xsurf, ysurf, zsurf)
852  IF (ALLOCATED(nxsurf))DEALLOCATE(nxsurf, nysurf, nzsurf)
853  IF (ALLOCATED(btopx))DEALLOCATE(btopx, btopy, btopz, btops)
854  IF (ALLOCATED(jx_3d))DEALLOCATE(jx_3d, jy_3d, jz_3d)
855 
856  IF (ezspline_allocated(x_spl)) CALL ezspline_free(x_spl,ier)
857  IF (ezspline_allocated(y_spl)) CALL ezspline_free(y_spl,ier)
858  IF (ezspline_allocated(z_spl)) CALL ezspline_free(z_spl,ier)
859  IF (ezspline_allocated(nx_spl)) CALL ezspline_free(nx_spl,ier)
860  IF (ezspline_allocated(ny_spl)) CALL ezspline_free(ny_spl,ier)
861  IF (ezspline_allocated(nz_spl)) CALL ezspline_free(nz_spl,ier)
862  IF (ezspline_allocated(kx_spl)) CALL ezspline_free(kx_spl,ier)
863  IF (ezspline_allocated(ky_spl)) CALL ezspline_free(ky_spl,ier)
864  IF (ezspline_allocated(kz_spl)) CALL ezspline_free(kz_spl,ier)
865  IF (ezspline_allocated(bn_spl)) CALL ezspline_free(bn_spl,ier)
866  IF (ezspline_allocated(x3d_spl)) CALL ezspline_free(x3d_spl,ier)
867  IF (ezspline_allocated(y3d_spl)) CALL ezspline_free(y3d_spl,ier)
868  IF (ezspline_allocated(z3d_spl)) CALL ezspline_free(z3d_spl,ier)
869  IF (ezspline_allocated(jx3d_spl)) CALL ezspline_free(jx3d_spl,ier)
870  IF (ezspline_allocated(jy3d_spl)) CALL ezspline_free(jy3d_spl,ier)
871  IF (ezspline_allocated(jz3d_spl)) CALL ezspline_free(jz3d_spl,ier)
872 
873  ! END SUBROUTINE
874  END SUBROUTINE free_virtual_casing
875  !-----------------------------------------------------------------
876 
877 
878  !-----------------------------------------------------------------
879  SUBROUTINE bfield_virtual_casing_dbl(x,y,z,bx,by,bz)
880  IMPLICIT NONE
881  ! INPUT VARIABLES
882  DOUBLE PRECISION, INTENT(in) :: x, y, z
883  DOUBLE PRECISION, INTENT(out) :: bx, by, bz
884  ! LOCAL VARIABLES
885  DOUBLE PRECISION :: gf(1:nuvp), gf3(1:nuvp)
886  ! BEGIN SUBROUTINE
887  gf(1:nuvp) = one/dsqrt( (x-xsurf(1:nuvp))**2 &
888  + (y-ysurf(1:nuvp))**2 &
889  + (z-zsurf(1:nuvp))**2)
890  gf3 = gf*gf*gf
891  bx = norm*sum( ( btopy(1:nuvp)*(z-zsurf(1:nuvp)) &
892  - btopz(1:nuvp)*(y-ysurf(1:nuvp)) ) * gf3(1:nuvp))
893  by = norm*sum( ( btopz(1:nuvp)*(x-xsurf(1:nuvp)) &
894  - btopx(1:nuvp)*(z-zsurf(1:nuvp)) ) * gf3(1:nuvp))
895  bz = norm*sum( ( btopx(1:nuvp)*(y-ysurf(1:nuvp)) &
896  - btopy(1:nuvp)*(x-xsurf(1:nuvp)) ) * gf3(1:nuvp))
897  bx = bx + norm*sum((btops(1:nuvp)*(x-xsurf(1:nuvp)))*gf3(1:nuvp))
898  by = by + norm*sum((btops(1:nuvp)*(y-ysurf(1:nuvp)))*gf3(1:nuvp))
899  bz = bz + norm*sum((btops(1:nuvp)*(z-zsurf(1:nuvp)))*gf3(1:nuvp))
900  RETURN
901  ! END SUBROUTINE
902  END SUBROUTINE bfield_virtual_casing_dbl
903  !-----------------------------------------------------------------
904 
905 
906  !-----------------------------------------------------------------
907  SUBROUTINE bfield_virtual_casing_flt(x_flt,y_flt,z_flt,bx_flt,by_flt,bz_flt)
908  IMPLICIT NONE
909  ! INPUT VARIABLES
910  REAL, INTENT(in) :: x_flt, y_flt, z_flt
911  REAL, INTENT(out) :: bx_flt, by_flt, bz_flt
912  ! LOCAL VARIABLES
913  DOUBLE PRECISION :: xt,yt,zt,bxt,byt,bzt
914  ! BEGIN SUBROUTINE
915  xt = x_flt
916  yt = y_flt
917  zt = z_flt
918  bxt = zero
919  byt = zero
920  bzt = zero
921  CALL bfield_virtual_casing_dbl(xt,yt,zt,bxt,byt,bzt)
922  bx_flt = bxt
923  by_flt = byt
924  bz_flt = bzt
925  RETURN
926  ! END SUBROUTINE
927  END SUBROUTINE bfield_virtual_casing_flt
928  !-----------------------------------------------------------------
929 
930 
931  !-----------------------------------------------------------------
932  SUBROUTINE bfield_virtual_casing_adapt_dbl(x,y,z,bx,by,bz,istat)
933  IMPLICIT NONE
934  ! INPUT VARIABLES
935  DOUBLE PRECISION, INTENT(in) :: x, y, z
936  DOUBLE PRECISION, INTENT(out) :: bx, by, bz
937  INTEGER, INTENT(inout) :: istat
938  ! LOCAL VARIABLES
939  LOGICAL :: adapt_rerun
940  INTEGER(KIND=8), PARAMETER :: ndim_nag = 2 ! theta,zeta
941  integer(kind=8), PARAMETER :: nfun_nag = 3 ! Bx, By, Bz
942  !INTEGER, PARAMETER :: lenwrk_nag = 6*ndim_nag+9*nfun_nag+(ndim_nag+nfun_nag+2)*(1+1)
943  INTEGER(KIND=8), PARAMETER :: lenwrk_nag = iwrk
944  INTEGER(KIND=8) :: maxcls_nag,mincls_nag, subs, restar, wrklen, rulcls, wrksbs, n, m, funcls
945  !DOUBLE PRECISION :: gf(1:nuvp), gf3(1:nuvp)
946  DOUBLE PRECISION :: absreq_nag, relreq_nag
947  DOUBLE PRECISION :: wrkstr_nag(lenwrk_nag)
948  DOUBLE PRECISION :: a_nag(ndim_nag), b_nag(ndim_nag), &
949  finest_nag(nfun_nag), absest_nag(nfun_nag)
950  DOUBLE PRECISION, ALLOCATABLE :: vrtwrk(:)
951 
952 #ifdef NAG
953  EXTERNAL :: d01eaf
954 
955 #else
956  EXTERNAL :: dcuhre
957 
958 #endif
959 
960  ! BEGIN SUBROUTINE
961  IF (adapt_tol < 0) THEN
962  CALL bfield_virtual_casing_dbl(x,y,z,bx,by,bz)
963  RETURN
964  END IF
965 
966  a_nag(1:2) = zero
967  b_nag(1:2) = one
968  mincls_nag = min_cls
969  maxcls_nag = iwrk
970  !absreq_nag = zero
971  absreq_nag = adapt_tol ! Talk to Stuart about proper values
972  relreq_nag = adapt_rel ! Talk to Stuart about proper values
973  finest_nag = zero
974  absest_nag = zero
975  x_nag = x
976  y_nag = y
977  z_nag = z
978  adapt_rerun = .true.
979  subs = 1
980  restar = 0
981  DO WHILE (adapt_rerun)
982 
983 #ifdef NAG
984  CALL d01eaf(ndim_nag,a_nag,b_nag,mincls_nag,maxcls_nag,nfun_nag,funsub_nag_b,absreq_nag,&
985  relreq_nag,lenwrk_nag,wrkstr_nag,finest_nag,absest_nag,istat)
986  IF (istat == 1 .or. istat == 3) THEN
987  maxcls_nag = maxcls_nag*10
988  mincls_nag = -1
989  restar = 1
990  WRITE(6,*) '!!!!! WARNING Could not reach desired tollerance !!!!!'
991  WRITE(6,*) ' BX = ',finest_nag(1),' +/- ',absest_nag(1)
992  WRITE(6,*) ' BY = ',finest_nag(2),' +/- ',absest_nag(2)
993  WRITE(6,*) ' BZ = ',finest_nag(3),' +/- ',absest_nag(3)
994  WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
995  ELSE IF (mincls_nag <= 32) THEN
996  mincls_nag = 64
997  restar = 1
998  ELSE IF (istat < 0) THEN
999  bx = zero
1000  by = zero
1001  bz = zero
1002  adapt_rerun=.false.
1003  ELSE
1004  bx = finest_nag(1)
1005  by = finest_nag(2)
1006  bz = finest_nag(3)
1007  adapt_rerun=.false.
1008  END IF
1009 
1010 #else
1011  IF (.not.ALLOCATED(vrtwrk)) THEN
1012  wrklen = ((maxcls_nag-ndim_nag)/(2*ndim_nag) + 1)*(2*ndim_nag+2*nfun_nag+2) + 17*nfun_nag + 256
1013  ALLOCATE(vrtwrk(wrklen),stat=istat)
1014  IF (istat .ne. 0) THEN
1015  WRITE(6,*) ' ALLOCATION ERROR IN: bfield_virtual_casing_adapt_dbl'
1016  WRITE(6,*) ' VARIABLE: VRTWRK, SIZE: ',wrklen
1017  WRITE(6,*) ' ISTAT: ',istat
1018  RETURN
1019  END IF
1020  END IF
1021  CALL dcuhre(ndim_nag,nfun_nag,a_nag,b_nag,mincls_nag,maxcls_nag,funsub_nag_b,absreq_nag,&
1022  relreq_nag,0,wrklen,restar,finest_nag,absest_nag,funcls,istat,vrtwrk)
1023  !PRINT *,istat,mincls_nag,funcls,maxcls_nag,finest_nag(1),finest_nag(2),finest_nag(3)
1024  IF (istat == 1) THEN
1025  maxcls_nag = maxcls_nag*10
1026  mincls_nag = funcls
1027  restar = 1
1028  istat = 0
1029  ELSE IF (istat > 1) THEN
1030  bx = zero
1031  by = zero
1032  bz = zero
1033  adapt_rerun=.false.
1034  DEALLOCATE(vrtwrk)
1035  ELSE
1036  bx = finest_nag(1)
1037  by = finest_nag(2)
1038  bz = finest_nag(3)
1039  adapt_rerun=.false.
1040  DEALLOCATE(vrtwrk)
1041  END IF
1042 
1043 #endif
1044  END DO
1045  nlastcall=mincls_nag
1046  RETURN
1047  ! END SUBROUTINE
1048  END SUBROUTINE bfield_virtual_casing_adapt_dbl
1049  !-----------------------------------------------------------------
1050 
1051 
1052  !-----------------------------------------------------------------
1053  SUBROUTINE bfield_virtual_casing_adapt_flt(x_flt,y_flt,z_flt,bx_flt,by_flt,bz_flt,istat)
1054  IMPLICIT NONE
1055  ! INPUT VARIABLES
1056  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1057  REAL, INTENT(out) :: bx_flt, by_flt, bz_flt
1058  INTEGER, INTENT(inout) :: istat
1059  ! LOCAL VARIABLES
1060  DOUBLE PRECISION :: xt,yt,zt,bxt,byt,bzt
1061  ! BEGIN SUBROUTINE
1062  xt = x_flt
1063  yt = y_flt
1064  zt = z_flt
1065  bxt = zero
1066  byt = zero
1067  bzt = zero
1068  CALL bfield_virtual_casing_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1069  bx_flt = bxt
1070  by_flt = byt
1071  bz_flt = bzt
1072  RETURN
1073  ! END SUBROUTINE
1074  END SUBROUTINE bfield_virtual_casing_adapt_flt
1075  !-----------------------------------------------------------------
1076 
1077 
1078  !-----------------------------------------------------------------
1079  SUBROUTINE funsub_nag_b(ndim, vec, nfun, f)
1080  IMPLICIT NONE
1081  ! INPUT VARIABLES
1082  INTEGER, INTENT(in) :: ndim, nfun
1083  DOUBLE PRECISION, INTENT(in) :: vec(ndim)
1084  DOUBLE PRECISION, INTENT(out) :: f(nfun)
1085  ! LOCAL VARIABLES
1086  INTEGER :: ier
1087  DOUBLE PRECISION :: bn, bx, by, bz , xs, ys, zs, gf, gf3, nx, ny, nz, ax, ay, az
1088  ! BEGIN SUBROUTINE
1089 
1090  xs = zero; ys = zero; zs = zero
1091  CALL ezspline_interp(x_spl,vec(1),vec(2),xs,ier)
1092  CALL ezspline_interp(y_spl,vec(1),vec(2),ys,ier)
1093  CALL ezspline_interp(z_spl,vec(1),vec(2),zs,ier)
1094  CALL ezspline_interp(kx_spl,vec(1),vec(2),ax,ier)
1095  CALL ezspline_interp(ky_spl,vec(1),vec(2),ay,ier)
1096  CALL ezspline_interp(kz_spl,vec(1),vec(2),az,ier)
1097 
1098  !CALL EZspline_interp(bn_spl,vec(1),vec(2),bn,ier)
1099  bn = zero
1100  gf = one/dsqrt((x_nag-xs)*(x_nag-xs)+(y_nag-ys)*(y_nag-ys)+(z_nag-zs)*(z_nag-zs))
1101  gf3 = gf*gf*gf
1102  f(1) = norm_nag*(ay*(z_nag-zs)-az*(y_nag-ys)+bn*(x_nag-xs))*gf3
1103  f(2) = norm_nag*(az*(x_nag-xs)-ax*(z_nag-zs)+bn*(y_nag-ys))*gf3
1104  f(3) = norm_nag*(ax*(y_nag-ys)-ay*(x_nag-xs)+bn*(z_nag-zs))*gf3
1105  !WRITE(427,*) vec(1),vec(2),xs,ys,zs
1106  RETURN
1107  ! END SUBROUTINE
1108  END SUBROUTINE funsub_nag_b
1109  !-----------------------------------------------------------------
1110 
1111 
1112  !-----------------------------------------------------------------
1113  SUBROUTINE vecpot_virtual_casing_dbl(x,y,z,ax,ay,az)
1114  IMPLICIT NONE
1115  ! INPUT VARIABLES
1116  DOUBLE PRECISION, INTENT(in) :: x, y, z
1117  DOUBLE PRECISION, INTENT(out) :: ax, ay,az
1118  ! LOCAL VARIABLES
1119  DOUBLE PRECISION :: gf(1:nuvp)
1120  ! BEGIN SUBROUTINE
1121  gf(1:nuvp) = one/dsqrt( (x-xsurf(1:nuvp))**2 &
1122  + (y-ysurf(1:nuvp))**2 &
1123  + (z-zsurf(1:nuvp))**2)
1124  ax = norm*sum((btopx(1:nuvp)+btops(1:nuvp)*nxsurf(1:nuvp))*gf(1:nuvp))
1125  ay = norm*sum((btopy(1:nuvp)+btops(1:nuvp)*nysurf(1:nuvp))*gf(1:nuvp))
1126  az = norm*sum((btopz(1:nuvp)+btops(1:nuvp)*nzsurf(1:nuvp))*gf(1:nuvp))
1127  ! END SUBROUTINE
1128  END SUBROUTINE vecpot_virtual_casing_dbl
1129  !-----------------------------------------------------------------
1130 
1131 
1132  !-----------------------------------------------------------------
1133  SUBROUTINE vecpot_virtual_casing_flt(x_flt,y_flt,z_flt,ax_flt,ay_flt,az_flt)
1134  IMPLICIT NONE
1135  ! INPUT VARIABLES
1136  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1137  REAL, INTENT(out) :: ax_flt, ay_flt,az_flt
1138  ! LOCAL VARIABLES
1139  DOUBLE PRECISION :: xt, yt, zt
1140  DOUBLE PRECISION :: axt, ayt,azt
1141  ! BEGIN SUBROUTINE
1142  xt = x_flt
1143  yt = y_flt
1144  zt = z_flt
1145  axt = zero
1146  ayt = zero
1147  azt = zero
1148  CALL vecpot_virtual_casing_dbl(xt,yt,zt,axt,ayt,azt)
1149  ax_flt = axt
1150  ay_flt = ayt
1151  az_flt = azt
1152  RETURN
1153  ! END SUBROUTINE
1154  END SUBROUTINE vecpot_virtual_casing_flt
1155  !-----------------------------------------------------------------
1156 
1157 
1158  !-----------------------------------------------------------------
1159  SUBROUTINE vecpot_virtual_casing_adapt_dbl(x,y,z,ax,ay,az,istat)
1160  IMPLICIT NONE
1161  ! INPUT VARIABLES
1162  DOUBLE PRECISION, INTENT(in) :: x, y, z
1163  DOUBLE PRECISION, INTENT(out) :: ax, ay, az
1164  INTEGER, INTENT(inout) :: istat
1165  ! LOCAL VARIABLES
1166  LOGICAL :: adapt_rerun
1167  INTEGER(KIND=8), PARAMETER :: ndim_nag = 2 ! theta,zeta
1168  integer(kind=8), PARAMETER :: nfun_nag = 3 ! Bx, By, Bz
1169  !INTEGER, PARAMETER :: lenwrk_nag = 6*ndim_nag+9*nfun_nag+(ndim_nag+nfun_nag+2)*(1+1)
1170  INTEGER(KIND=8), PARAMETER :: lenwrk_nag = iwrk
1171  INTEGER(KIND=8) :: maxcls_nag,mincls_nag, subs, restar, wrklen, rulcls, wrksbs, n, m, funcls
1172  !DOUBLE PRECISION :: gf(1:nuvp), gf3(1:nuvp)
1173  DOUBLE PRECISION :: absreq_nag, relreq_nag
1174  DOUBLE PRECISION :: wrkstr_nag(lenwrk_nag)
1175  DOUBLE PRECISION :: a_nag(ndim_nag), b_nag(ndim_nag), &
1176  finest_nag(nfun_nag), absest_nag(nfun_nag)
1177  DOUBLE PRECISION, ALLOCATABLE :: vrtwrk(:)
1178 
1179 #ifdef NAG
1180  EXTERNAL :: d01eaf
1181 
1182 #else
1183  EXTERNAL :: dcuhre
1184 
1185 #endif
1186  ! BEGIN SUBROUTINE
1187  IF (adapt_tol < 0) THEN
1188  CALL vecpot_virtual_casing_dbl(x,y,z,ax,ay,az)
1189  RETURN
1190  END IF
1191 
1192  a_nag(1:2) = zero
1193  b_nag(1:2) = one
1194  mincls_nag = min_cls
1195  maxcls_nag = iwrk
1196  absreq_nag = adapt_tol ! Talk to Stuart about proper values
1197  relreq_nag = adapt_rel ! Talk to Stuart about proper values
1198  finest_nag = zero
1199  absest_nag = zero
1200  x_nag = x
1201  y_nag = y
1202  z_nag = z
1203  adapt_rerun = .true.
1204  subs = 1
1205  restar = 0
1206  DO WHILE (adapt_rerun)
1207 
1208 #ifdef NAG
1209  CALL d01eaf(ndim_nag,a_nag,b_nag,mincls_nag,maxcls_nag,nfun_nag,funsub_nag_a,absreq_nag,&
1210  relreq_nag,lenwrk_nag,wrkstr_nag,finest_nag,absest_nag,istat)
1211  IF (istat == 1 .or. istat == 3) THEN
1212  maxcls_nag = maxcls_nag*10
1213  mincls_nag = -1
1214  restar = 1
1215  WRITE(6,*) '!!!!! WARNING Could not reach desired tollerance !!!!!'
1216  WRITE(6,*) ' AX = ',finest_nag(1),' +/- ',absest_nag(1)
1217  WRITE(6,*) ' AY = ',finest_nag(2),' +/- ',absest_nag(2)
1218  WRITE(6,*) ' AZ = ',finest_nag(3),' +/- ',absest_nag(3)
1219  WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
1220  ELSE IF (mincls_nag <= 32) THEN
1221  mincls_nag = 64
1222  restar = 1
1223  ELSE IF (istat < 0) THEN
1224  ax = zero
1225  ay = zero
1226  az = zero
1227  adapt_rerun=.false.
1228  print *,'error'
1229  ELSE
1230  ax = finest_nag(1)
1231  ay = finest_nag(2)
1232  az = finest_nag(3)
1233  adapt_rerun=.false.
1234  END IF
1235 
1236 #else
1237  IF (.not.ALLOCATED(vrtwrk)) THEN
1238  wrklen = ((iwrk-ndim_nag)/(2*ndim_nag) + 1)*(2*ndim_nag+2*nfun_nag+2) + 17*nfun_nag + 256
1239  ALLOCATE(vrtwrk(wrklen),stat=istat)
1240  IF (istat .ne. 0) THEN
1241  WRITE(6,*) ' ALLOCATION ERROR IN: vecpot_virtual_casing_adapt_dbl'
1242  WRITE(6,*) ' VARIABLE: VRTWRK, SIZE: ',wrklen
1243  WRITE(6,*) ' ISTAT: ',istat
1244  RETURN
1245  END IF
1246  END IF
1247  CALL dcuhre(ndim_nag,nfun_nag,a_nag,b_nag,mincls_nag,maxcls_nag,funsub_nag_a,absreq_nag,&
1248  relreq_nag,0,wrklen,restar,finest_nag,absest_nag,funcls,istat,vrtwrk)
1249  !DEALLOCATE(vrtwrk)
1250  IF (istat == 1) THEN
1251  maxcls_nag = maxcls_nag*10
1252  mincls_nag = -1
1253  restar = 1
1254  ELSE IF (istat > 1) THEN
1255  ax = zero
1256  ay = zero
1257  az = zero
1258  adapt_rerun=.false.
1259  DEALLOCATE(vrtwrk)
1260  ELSE
1261  ax = finest_nag(1)
1262  ay = finest_nag(2)
1263  az = finest_nag(3)
1264  adapt_rerun=.false.
1265  DEALLOCATE(vrtwrk)
1266  END IF
1267 
1268 #endif
1269  END DO
1270  nlastcall=mincls_nag
1271  RETURN
1272  ! END SUBROUTINE
1273  END SUBROUTINE vecpot_virtual_casing_adapt_dbl
1274  !-----------------------------------------------------------------
1275 
1276 
1277  !-----------------------------------------------------------------
1278  SUBROUTINE vecpot_virtual_casing_adapt_flt(x_flt,y_flt,z_flt,ax_flt,ay_flt,az_flt,istat)
1279  IMPLICIT NONE
1280  ! INPUT VARIABLES
1281  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1282  REAL, INTENT(out) :: ax_flt, ay_flt,az_flt
1283  INTEGER, INTENT(inout) :: istat
1284  ! LOCAL VARIABLES
1285  DOUBLE PRECISION :: xt, yt, zt
1286  DOUBLE PRECISION :: axt, ayt,azt
1287  ! BEGIN SUBROUTINE
1288  xt = x_flt
1289  yt = y_flt
1290  zt = z_flt
1291  axt = zero
1292  ayt = zero
1293  azt = zero
1294  CALL vecpot_virtual_casing_adapt_dbl(xt,yt,zt,axt,ayt,azt,istat)
1295  ax_flt = axt
1296  ay_flt = ayt
1297  az_flt = azt
1298  RETURN
1299  ! END SUBROUTINE
1300  END SUBROUTINE vecpot_virtual_casing_adapt_flt
1301  !-----------------------------------------------------------------
1302 
1303 
1304  !-----------------------------------------------------------------
1305  SUBROUTINE funsub_nag_a(ndim, vec, nfun, f)
1306  IMPLICIT NONE
1307  ! INPUT VARIABLES
1308  INTEGER, INTENT(in) :: ndim, nfun
1309  DOUBLE PRECISION, INTENT(in) :: vec(ndim)
1310  DOUBLE PRECISION, INTENT(out) :: f(nfun)
1311  ! LOCAL VARIABLES
1312  INTEGER :: ier
1313  DOUBLE PRECISION :: bn, bx, by, bz , xs, ys, zs, gf, nx, ny, nz, ax, ay, az
1314  ! BEGIN SUBROUTINE
1315 
1316  xs = zero; ys = zero; zs = zero; ax = zero; ay = zero; az = zero
1317  nx = zero; ny = zero; nz = zero; bn = zero
1318  CALL ezspline_interp(x_spl,vec(1),vec(2),xs,ier)
1319  CALL ezspline_interp(y_spl,vec(1),vec(2),ys,ier)
1320  CALL ezspline_interp(z_spl,vec(1),vec(2),zs,ier)
1321  CALL ezspline_interp(kx_spl,vec(1),vec(2),ax,ier)
1322  CALL ezspline_interp(ky_spl,vec(1),vec(2),ay,ier)
1323  CALL ezspline_interp(kz_spl,vec(1),vec(2),az,ier)
1324  CALL ezspline_interp(nx_spl,vec(1),vec(2),nx,ier)
1325  CALL ezspline_interp(ny_spl,vec(1),vec(2),ny,ier)
1326  CALL ezspline_interp(nz_spl,vec(1),vec(2),nz,ier)
1327  CALL ezspline_interp(bn_spl,vec(1),vec(2),bn,ier)
1328 
1329  gf = one/dsqrt((x_nag-xs)*(x_nag-xs)+(y_nag-ys)*(y_nag-ys)+(z_nag-zs)*(z_nag-zs))
1330  f(1) = norm_nag*(ax+bn*nx)*gf
1331  f(2) = norm_nag*(ay+bn*ny)*gf
1332  f(3) = norm_nag*(az+bn*nz)*gf
1333  RETURN
1334  ! END SUBROUTINE
1335  END SUBROUTINE funsub_nag_a
1336  !-----------------------------------------------------------------
1337 
1338 
1339  !-----------------------------------------------------------------
1340  SUBROUTINE virtual_casing_info(iunit)
1341  IMPLICIT NONE
1342  ! INPUT VARIABLES
1343  INTEGER, INTENT(in) :: iunit
1344  ! LOCAL VARIABLES
1345  ! BEGIN SUBROUTINE
1346  WRITE(iunit,'(A)') '----- Virtual Casing Information -----'
1347 #ifdef NAG
1348  WRITE(iunit,'(A,A,A)') ' INTEGRAL TYPE: ',trim(vc_type_str),' (NAG) '
1349 #else
1350  WRITE(iunit,'(A,A,A)') ' INTEGRAL TYPE: ',trim(vc_type_str),' (DCUHRE) '
1351 #endif
1352  WRITE(iunit,'(A,ES11.4)') ' MIN_GRID_DISTANCE = ',min_delta_x
1353  WRITE(iunit,'(A,ES11.4)') ' NORMAL_AREA = ',surf_area
1354  WRITE(iunit,'(A,I4,A,I4,A,I4,A,I3)') ' NR = ',nr_vc,'; NU = ',nu_vc,'; NV = ',nv_vc,'; NFP = ',nvp/nv_vc
1355  WRITE(iunit,'(A,I6)') ' NUVP = ',nuvp
1356  IF (adapt_tol > 0 .or. adapt_rel >0) THEN
1357  WRITE(iunit,'(A,ES11.4,A,ES11.4)') ' ABS_TOL = ',adapt_tol,'; REL_TOL = ',adapt_rel
1358  ELSE
1359  WRITE(iunit,'(A)') ' !!!!! Using discrete integration !!!!!'
1360  END IF
1361  WRITE(iunit,'(A,I6,A,I8,A)') ' MIN_CLS = ',min_cls,' (',iwrk,')'
1362  CALL flush(iunit)
1363  ! END SUBROUTINE
1364  END SUBROUTINE virtual_casing_info
1365  !-----------------------------------------------------------------
1366 
1367 
1368  !-----------------------------------------------------------------
1369  SUBROUTINE virtual_casing_surf_dump(iunit)
1370  IMPLICIT NONE
1371  ! INPUT VARIABLES
1372  INTEGER, INTENT(in) :: iunit
1373  ! LOCAL VARIABLES
1374  INTEGER :: i, j, k, ik
1375  ! BEGIN SUBROUTINE
1376  ik = 1
1377  WRITE(iunit,'(A)') ' n nfp nu nv x y z nx ny nz btopx btopy btopz btops'
1378  DO i = 1, nuvp/(nu_vc*nv_vc)
1379  DO k = 1, nu_vc
1380  DO j = 1, nv_vc
1381  WRITE(iunit,'(4I6,10ES23.16)') ik, i, k, j, &
1382  xsurf(ik), ysurf(ik), zsurf(ik), &
1383  nxsurf(ik), nysurf(ik), nzsurf(ik), &
1384  btopx(ik), btopy(ik), btopz(ik), &
1385  btops(ik)
1386  ik = ik + 1
1387  END DO
1388  END DO
1389  END DO
1390  CALL flush(iunit)
1391  ! END SUBROUTINE
1392  END SUBROUTINE virtual_casing_surf_dump
1393  !-----------------------------------------------------------------
1394 
1395 
1396  !-----------------------------------------------------------------
1397  DOUBLE PRECISION FUNCTION virtual_casing_dist_dbl(xp,yp,zp)
1398  IMPLICIT NONE
1399  ! INPUT VARIABLES
1400  DOUBLE PRECISION, INTENT(in) :: xp,yp,zp
1401  ! LOCAL VARIABLES
1402  ! BEGIN FUNCTION
1403  virtual_casing_dist_dbl = minval(dsqrt((xp-xsurf)**2+(yp-ysurf)**2+(zp-zsurf)**2))
1404  ! END FUNCTION
1405  END FUNCTION virtual_casing_dist_dbl
1406  !-----------------------------------------------------------------
1407 
1408 
1409  !-----------------------------------------------------------------
1410  REAL FUNCTION virtual_casing_dist_flt(xp_flt,yp_flt,zp_flt)
1411  IMPLICIT NONE
1412  ! INPUT VARIABLES
1413  REAL, INTENT(in) :: xp_flt,yp_flt,zp_flt
1414  ! LOCAL VARIABLES
1415  ! BEGIN FUNCTION
1416  virtual_casing_dist_flt = minval(dsqrt((xp_flt-xsurf)**2+(yp_flt-ysurf)**2+(zp_flt-zsurf)**2))
1417  ! END FUNCTION
1418  END FUNCTION virtual_casing_dist_flt
1419  !-----------------------------------------------------------------
1420 
1421 
1422  !-----------------------------------------------------------------
1423  SUBROUTINE init_volint_dbl(mnmax,nu,nv,ns,xm,xn,rmnc,zmns,nfp,jumnc,jvmnc,&
1424  rmns,zmnc,jumns,jvmns)
1425  IMPLICIT NONE
1426  ! INPUT VARIABLES
1427  INTEGER, INTENT(in) :: mnmax, nu, nv, nfp,ns
1428  INTEGER, INTENT(in) :: xm(1:mnmax),xn(1:mnmax)
1429  DOUBLE PRECISION, INTENT(in) :: rmnc(1:mnmax,ns),zmns(1:mnmax,ns)
1430  DOUBLE PRECISION, INTENT(in) :: jumnc(1:mnmax,ns),jvmnc(1:mnmax,ns)
1431  DOUBLE PRECISION, INTENT(in), OPTIONAL :: jumns(1:mnmax,ns),jvmns(1:mnmax,ns)
1432  DOUBLE PRECISION, INTENT(in), OPTIONAL :: rmns(1:mnmax,ns),zmnc(1:mnmax,ns)
1433  ! LOCAL VARIABLES
1434  INTEGER :: mn, uv, i, u, v, dex1, dex2, ier, nuv, k
1435  INTEGER :: bcs1(2), bcs2(2), bcs3(2)
1436  DOUBLE PRECISION :: cop, sip
1437  DOUBLE PRECISION, ALLOCATABLE :: xu(:), xv(:)
1438  DOUBLE PRECISION, ALLOCATABLE :: fmn_temp(:,:)
1439  DOUBLE PRECISION, ALLOCATABLE :: r_temp(:,:,:), z_temp(:,:,:)
1440  DOUBLE PRECISION, ALLOCATABLE :: x_temp(:,:,:), y_temp(:,:,:)
1441  DOUBLE PRECISION, ALLOCATABLE :: jv_temp(:,:,:), ju_temp(:,:,:)
1442  DOUBLE PRECISION, ALLOCATABLE :: ru(:,:,:), rv(:,:,:)
1443  DOUBLE PRECISION, ALLOCATABLE :: zu(:,:,:), zv(:,:,:)
1444  DOUBLE PRECISION, ALLOCATABLE :: jr(:,:,:), jphi(:,:,:), jz(:,:,:)
1445  DOUBLE PRECISION, ALLOCATABLE :: jx(:,:,:), jy(:,:,:)
1446  ! BEGIN SUBROUTINE
1447  ! Deallocate anything globally allocated
1448 
1449  IF (ezspline_allocated(x3d_spl)) CALL ezspline_free(x3d_spl,ier)
1450  IF (ezspline_allocated(y3d_spl)) CALL ezspline_free(y3d_spl,ier)
1451  IF (ezspline_allocated(z3d_spl)) CALL ezspline_free(z3d_spl,ier)
1452  IF (ezspline_allocated(jx3d_spl)) CALL ezspline_free(jx3d_spl,ier)
1453  IF (ezspline_allocated(jy3d_spl)) CALL ezspline_free(jy3d_spl,ier)
1454  IF (ezspline_allocated(jz3d_spl)) CALL ezspline_free(jz3d_spl,ier)
1455 
1456  ! Initialize varaibles
1457  bcs1=(/-1,-1/)
1458  bcs2=(/-1,-1/)
1459  bcs3=(/0,0/)
1460  adapt_tol = 1.0d-5
1461  adapt_rel = 1.0d-3
1462  nr_vc = ns
1463  nu_vc = nu
1464  nv_vc = nv
1465  nuv = nu * nv
1466  nvp = nv * nfp
1467  nuvp = nu * nv * nfp
1468  norm_3d = 1.0d-7
1469  vc_type_str='Volume Integral'
1470  ! Alloocations
1471  ALLOCATE(xu(nu),xv(nvp))
1472  ALLOCATE(fmn_temp(mnmax,ns))
1473  ALLOCATE(r_temp(nu,nvp,ns),z_temp(nu,nvp,ns))
1474  ALLOCATE(x_temp(nu,nvp,ns),y_temp(nu,nvp,ns))
1475  ALLOCATE(ju_temp(nu,nvp,ns),jv_temp(nu,nvp,ns))
1476  ALLOCATE(ru(nu,nvp,ns),rv(nu,nvp,ns))
1477  ALLOCATE(zu(nu,nvp,ns),zv(nu,nvp,ns))
1478  ALLOCATE(jx(nu,nvp,ns),jy(nu,nvp,ns))
1479  ALLOCATE(jr(nu,nvp,ns),jphi(nu,nvp,ns),jz(nu,nvp,ns))
1480  ALLOCATE(jx_3d(nu*nvp*ns),jy_3d(nu*nvp*ns),jz_3d(nu*nvp*ns))
1481  ALLOCATE(xsurf(nu*nvp*ns),ysurf(nu*nvp*ns),zsurf(nu*nvp*ns))
1482 
1483  ! Get the major quantities (not require use the VMEC convenction of j of j=j*jac
1484  r_temp=zero; z_temp=zero; ju_temp=zero; jv_temp=zero
1485  FORALL(u=1:nu) xu(u) = dble(u-1)/dble(nu-1)
1486  FORALL(v=1:nvp) xv(v) = dble(v-1)/dble(nvp-1)
1487  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,rmnc,xm,xn,r_temp,0,1)
1488  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,zmns,xm,xn,z_temp,1,0)
1489  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,jumnc,xm,xn,ju_temp,0,0)
1490  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,jvmnc,xm,xn,jv_temp,0,0)
1491  IF (PRESENT(rmns)) CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,rmns,xm,xn,r_temp,1,0)
1492  IF (PRESENT(zmnc)) CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,zmnc,xm,xn,z_temp,0,0)
1493  IF (PRESENT(jumns)) CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,jumns,xm,xn,ju_temp,1,0)
1494  IF (PRESENT(jvmns)) CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,jvmns,xm,xn,jv_temp,1,0)
1495  ! Now we calculate the edge metric elements
1496  ru = zero; zu = zero; rv = zero; zv = zero
1497  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = -rmnc(mn,:)*xm(mn)
1498  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,ru,1,0)
1499  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = -rmnc(mn,:)*xn(mn)
1500  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,rv,1,0)
1501  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = zmns(mn,:)*xm(mn)
1502  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,zu,0,0)
1503  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = zmns(mn,:)*xn(mn)
1504  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,zv,0,0)
1505  IF (PRESENT(rmns)) THEN
1506  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = rmns(mn,:)*xm(mn)
1507  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,ru,0,0)
1508  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = rmns(mn,:)*xn(mn)
1509  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,rv,0,0)
1510  END IF
1511  IF (PRESENT(zmnc)) THEN
1512  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = -zmnc(mn,:)*xm(mn)
1513  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,zu,1,0)
1514  FORALL(mn = 1:mnmax) fmn_temp(mn,:) = -zmnc(mn,:)*xn(mn)
1515  CALL mntouv_local(1,ns,mnmax,nu,nvp,xu,xv,fmn_temp,xm,xn,zv,1,0)
1516  END IF
1517  ! Calculate jr, jphi, jz, jx, jy
1518  jr = ju_temp*ru+jv_temp*rv
1519  jphi = r_temp * jv_temp
1520  jz = ju_temp*zu+jv_temp*zv
1521  DO u = 1, nu
1522  DO v = 1, nvp
1523  cop = dcos(pi2*xv(v))
1524  sip = dsin(pi2*xv(v))
1525  x_temp(u,v,:) = r_temp(u,v,:) * cop
1526  y_temp(u,v,:) = r_temp(u,v,:) * sip
1527  jx(u,v,:) = jr(u,v,:) * cop - jphi(u,v,:) * sip
1528  jy(u,v,:) = jr(u,v,:) * sip + jphi(u,v,:) * cop
1529  END DO
1530  END DO
1531  i=1
1532  DO u = 1, nu
1533  DO v = 1, nvp
1534  DO k = 1, ns
1535  jx_3d(i) = jx(u,v,k)
1536  jy_3d(i) = jy(u,v,k)
1537  jz_3d(i) = jz(u,v,k)
1538  xsurf(i) = x_temp(u,v,k)
1539  ysurf(i) = y_temp(u,v,k)
1540  zsurf(i) = z_temp(u,v,k)
1541  i = i + 1
1542  END DO
1543  END DO
1544  END DO
1545 
1546  ! Construct Splines
1547  CALL ezspline_init(x3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1548  CALL ezspline_init(y3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1549  CALL ezspline_init(z3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1550  CALL ezspline_init(jx3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1551  CALL ezspline_init(jy3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1552  CALL ezspline_init(jz3d_spl,nu,nvp,ns,bcs1,bcs2,bcs3,ier)
1553  x3d_spl%x1=xu*pi2
1554  y3d_spl%x1=xu*pi2
1555  z3d_spl%x1=xu*pi2
1556  jx3d_spl%x1=xu*pi2
1557  jy3d_spl%x1=xu*pi2
1558  jz3d_spl%x1=xu*pi2
1559  x3d_spl%x2=xv*pi2
1560  y3d_spl%x2=xv*pi2
1561  z3d_spl%x2=xv*pi2
1562  jx3d_spl%x2=xv*pi2
1563  jy3d_spl%x2=xv*pi2
1564  jz3d_spl%x2=xv*pi2
1565  ! default [0 1] x3
1566  x3d_spl%isHermite = 1
1567  y3d_spl%isHermite = 1
1568  z3d_spl%isHermite = 1
1569  jx3d_spl%isHermite = 1
1570  jy3d_spl%isHermite = 1
1571  jz3d_spl%isHermite = 1
1572  CALL ezspline_setup(x3d_spl,x_temp,ier)
1573  CALL ezspline_setup(y3d_spl,y_temp,ier)
1574  CALL ezspline_setup(z3d_spl,z_temp,ier)
1575  CALL ezspline_setup(jx3d_spl,jx,ier)
1576  CALL ezspline_setup(jy3d_spl,jy,ier)
1577  CALL ezspline_setup(jz3d_spl,jz,ier)
1578 
1579  ! DEALLOCATE
1580  DEALLOCATE(xu,xv)
1581  DEALLOCATE(fmn_temp)
1582  DEALLOCATE(r_temp,z_temp)
1583  DEALLOCATE(x_temp,y_temp)
1584  DEALLOCATE(ju_temp,jv_temp)
1585  DEALLOCATE(ru,rv)
1586  DEALLOCATE(zu,zv)
1587  DEALLOCATE(jx,jy,jr,jphi,jz)
1588  ! END SUBROUTINE
1589  END SUBROUTINE init_volint_dbl
1590 
1591  !-----------------------------------------------------------------
1592 
1593  !-----------------------------------------------------------------
1594  ! Note optional arguments must have a different name so the
1595  ! Interface will pick the proper variables. Here: _FLT
1596  SUBROUTINE init_volint_flt(mnmax_flt,nu_flt,nv_flt,ns_flt,xm_flt,&
1597  xn_flt,rmnc_flt,zmns_flt,nfp_flt,&
1598  jumnc_flt,jvmnc_flt,&
1599  rmns_flt,zmnc_flt,&
1600  jumns_flt,jvmns_flt)
1601  IMPLICIT NONE
1602  ! INPUT VARIABLES
1603  INTEGER, INTENT(in) :: mnmax_flt, nu_flt, nv_flt, nfp_flt, ns_flt
1604  INTEGER, INTENT(in) :: xm_flt(1:mnmax_flt),xn_flt(1:mnmax_flt)
1605  REAL, INTENT(in) :: rmnc_flt(1:mnmax_flt,ns_flt),zmns_flt(1:mnmax_flt,ns_flt)
1606  REAL, INTENT(in) :: jumnc_flt(1:mnmax_flt,ns_flt),jvmnc_flt(1:mnmax_flt,ns_flt)
1607  REAL, INTENT(in), OPTIONAL :: jumns_flt(1:mnmax_flt,ns_flt),jvmns_flt(1:mnmax_flt,ns_flt)
1608  REAL, INTENT(in), OPTIONAL :: rmns_flt(1:mnmax_flt,ns_flt),zmnc_flt(1:mnmax_flt,ns_flt)
1609  ! LOCAL VARIABLES
1610  DOUBLE PRECISION :: rmnct(1:mnmax_flt,ns_flt),zmnst(1:mnmax_flt,ns_flt)
1611  DOUBLE PRECISION :: rmnst(1:mnmax_flt,ns_flt),zmnct(1:mnmax_flt,ns_flt)
1612  DOUBLE PRECISION :: jumnct(1:mnmax_flt,ns_flt),jvmnct(1:mnmax_flt,ns_flt)
1613  DOUBLE PRECISION :: jumnst(1:mnmax_flt,ns_flt),jvmnst(1:mnmax_flt,ns_flt)
1614  ! BEGIN SUBROUTINE
1615  rmnct = zero; rmnst = zero; zmnct = zero; zmnst = zero
1616  jumnct = zero; jumnst = zero
1617  jvmnct = zero; jvmnst = zero
1618  rmnct = rmnc_flt
1619  zmnst = zmns_flt
1620  jumnct = jumnc_flt
1621  jvmnct = jvmnc_flt
1622  IF (PRESENT(rmns_flt)) rmnst = rmns_flt
1623  IF (PRESENT(zmnc_flt)) zmnct = zmnc_flt
1624  IF (PRESENT(jumns_flt)) jumnst = jumns_flt
1625  IF (PRESENT(jvmns_flt)) jvmnst = jvmns_flt
1626  CALL init_volint_dbl(mnmax_flt,nu_flt,nv_flt,ns_flt,xm_flt,xn_flt, &
1627  rmnct,zmnst,nfp_flt, &
1628  rmns=rmnst,zmnc=zmnct, &
1629  jumnc=jumnct, jumns=jumnst, &
1630  jvmnc=jvmnct, jvmns=jvmnst)
1631  !CALL init_volint_dbl(mnmax_flt,nu_flt,nv_flt,ns_flt,xm_flt,xn_flt, &
1632  ! rmnct,zmnst,nfp_flt, &
1633  ! JUMNC=jumnct, JVMNC=jvmnct)
1634  ! END SUBROUTINE
1635  END SUBROUTINE init_volint_flt
1636  !-----------------------------------------------------------------
1637 
1638 
1639  !-----------------------------------------------------------------
1640  SUBROUTINE funsub_nag_a3d(ndim, vec, nfun, f)
1641  IMPLICIT NONE
1642  ! INPUT VARIABLES
1643  INTEGER, INTENT(in) :: ndim, nfun
1644  DOUBLE PRECISION, INTENT(in) :: vec(ndim)
1645  DOUBLE PRECISION, INTENT(out) :: f(nfun)
1646  ! LOCAL VARIABLES
1647  INTEGER :: ier
1648  DOUBLE PRECISION :: ax, ay, az , xs, ys, zs, gf, nx, ny, nz
1649  ! BEGIN SUBROUTINE
1650 
1651  xs = zero; ys = zero; zs = zero
1652  ax = zero; ay = zero; az = zero
1653  CALL ezspline_interp(x3d_spl,vec(1),vec(2),vec(3),xs,ier)
1654  CALL ezspline_interp(y3d_spl,vec(1),vec(2),vec(3),ys,ier)
1655  CALL ezspline_interp(z3d_spl,vec(1),vec(2),vec(3),zs,ier)
1656  CALL ezspline_interp(jx3d_spl,vec(1),vec(2),vec(3),ax,ier)
1657  CALL ezspline_interp(jy3d_spl,vec(1),vec(2),vec(3),ay,ier)
1658  CALL ezspline_interp(jz3d_spl,vec(1),vec(2),vec(3),az,ier)
1659 
1660  gf = one/dsqrt((x_nag-xs)*(x_nag-xs)+(y_nag-ys)*(y_nag-ys)+(z_nag-zs)*(z_nag-zs))
1661  !PRINT *,xs,ys,zs,ax,ay,ax,gf
1662  f(1) = norm_3d*ax*gf
1663  f(2) = norm_3d*ay*gf
1664  f(3) = norm_3d*az*gf
1665  RETURN
1666  ! END SUBROUTINE
1667  END SUBROUTINE funsub_nag_a3d
1668  !-----------------------------------------------------------------
1669 
1670 
1671  !-----------------------------------------------------------------
1672  SUBROUTINE funsub_nag_b3d(ndim, vec, nfun, f)
1673  IMPLICIT NONE
1674  ! INPUT VARIABLES
1675  INTEGER, INTENT(in) :: ndim, nfun
1676  DOUBLE PRECISION, INTENT(in) :: vec(ndim)
1677  DOUBLE PRECISION, INTENT(out) :: f(nfun)
1678  ! LOCAL VARIABLES
1679  INTEGER :: ier
1680  DOUBLE PRECISION :: bn, ax, ay, az , xs, ys, zs, gf, gf3
1681  ! BEGIN SUBROUTINE
1682 
1683  xs = zero; ys = zero; zs = zero
1684  ax = zero; ay = zero; az = zero
1685  CALL ezspline_interp(x3d_spl,vec(1),vec(2),vec(3),xs,ier)
1686  CALL ezspline_interp(y3d_spl,vec(1),vec(2),vec(3),ys,ier)
1687  CALL ezspline_interp(z3d_spl,vec(1),vec(2),vec(3),zs,ier)
1688  CALL ezspline_interp(jx3d_spl,vec(1),vec(2),vec(3),ax,ier)
1689  CALL ezspline_interp(jy3d_spl,vec(1),vec(2),vec(3),ay,ier)
1690  CALL ezspline_interp(jz3d_spl,vec(1),vec(2),vec(3),az,ier)
1691 
1692  gf = one/dsqrt((x_nag-xs)*(x_nag-xs)+(y_nag-ys)*(y_nag-ys)+(z_nag-zs)*(z_nag-zs))
1693  gf3 = gf*gf*gf
1694  f(1) = norm_3d*(ay*(z_nag-zs)-az*(y_nag-ys))*gf3
1695  f(2) = norm_3d*(az*(x_nag-xs)-ax*(z_nag-zs))*gf3
1696  f(3) = norm_3d*(ax*(y_nag-ys)-ay*(x_nag-xs))*gf3
1697  !WRITE(427,*) xs,ys,zs
1698  RETURN
1699  ! END SUBROUTINE
1700  END SUBROUTINE funsub_nag_b3d
1701  !-----------------------------------------------------------------
1702 
1703 
1704  !-----------------------------------------------------------------
1705  SUBROUTINE vecpot_volint_adapt_flt(x_flt,y_flt,z_flt,ax_flt,ay_flt,az_flt,istat)
1706  IMPLICIT NONE
1707  ! INPUT VARIABLES
1708  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1709  REAL, INTENT(out) :: ax_flt, ay_flt,az_flt
1710  INTEGER, INTENT(inout) :: istat
1711  ! LOCAL VARIABLES
1712  DOUBLE PRECISION :: xt, yt, zt
1713  DOUBLE PRECISION :: axt, ayt,azt
1714  ! BEGIN SUBROUTINE
1715  xt = x_flt
1716  yt = y_flt
1717  zt = z_flt
1718  axt = zero
1719  ayt = zero
1720  azt = zero
1721  CALL vecpot_volint_adapt_dbl(xt,yt,zt,axt,ayt,azt,istat)
1722  ax_flt = axt
1723  ay_flt = ayt
1724  az_flt = azt
1725  RETURN
1726  ! END SUBROUTINE
1727  END SUBROUTINE vecpot_volint_adapt_flt
1728  !-----------------------------------------------------------------
1729 
1730 
1731  !-----------------------------------------------------------------
1732  SUBROUTINE vecpot_volint_dbl(x,y,z,ax,ay,az)
1733  IMPLICIT NONE
1734  ! INPUT VARIABLES
1735  DOUBLE PRECISION, INTENT(in) :: x, y, z
1736  DOUBLE PRECISION, INTENT(out) :: ax, ay, az
1737  ! LOCAL VARIABLES
1738  DOUBLE PRECISION :: gf(nuvp*nr_vc)
1739  ! BEGIN SUBROUTINE
1740  gf = one/dsqrt((x-xsurf)**2 + (y-ysurf)**2 + (z-zsurf)**2)
1741  ax = norm_3d*sum(jx_3d*gf)/nuvp
1742  ay = norm_3d*sum(jy_3d*gf)/nuvp
1743  az = norm_3d*sum(jz_3d*gf)/nuvp
1744  RETURN
1745  ! END SUBROUTINE
1746  END SUBROUTINE vecpot_volint_dbl
1747  !-----------------------------------------------------------------
1748 
1749 
1750  !-----------------------------------------------------------------
1751  SUBROUTINE bfield_volint_dbl(x,y,z,bx,by,bz)
1752  IMPLICIT NONE
1753  ! INPUT VARIABLES
1754  DOUBLE PRECISION, INTENT(in) :: x, y, z
1755  DOUBLE PRECISION, INTENT(out) :: bx, by, bz
1756  ! LOCAL VARIABLES
1757  DOUBLE PRECISION :: gf(nuvp*nr_vc), gf3(nuvp*nr_vc)
1758  ! BEGIN SUBROUTINE
1759  gf = one/dsqrt((x-xsurf)**2 + (y-ysurf)**2 + (z-zsurf)**2)
1760  gf3 = gf*gf*gf
1761  bx = norm_3d*sum((jy_3d*(z-zsurf)-jz_3d*(y-ysurf))*gf3)/nuvp
1762  by = norm_3d*sum((jz_3d*(x-xsurf)-jx_3d*(z-zsurf))*gf3)/nuvp
1763  bz = norm_3d*sum((jx_3d*(y-ysurf)-jy_3d*(x-xsurf))*gf3)/nuvp
1764  RETURN
1765  ! END SUBROUTINE
1766  END SUBROUTINE bfield_volint_dbl
1767  !-----------------------------------------------------------------
1768 
1769 
1770  !-----------------------------------------------------------------
1771  SUBROUTINE vecpot_volint_adapt_dbl(x,y,z,ax,ay,az,istat)
1772  IMPLICIT NONE
1773  ! INPUT VARIABLES
1774  DOUBLE PRECISION, INTENT(in) :: x, y, z
1775  DOUBLE PRECISION, INTENT(out) :: ax, ay, az
1776  INTEGER, INTENT(inout) :: istat
1777  ! LOCAL VARIABLES
1778  LOGICAL :: adapt_rerun
1779  INTEGER(KIND=8), PARAMETER :: ndim_nag = 3 ! theta,zeta,ns
1780  INTEGER(KIND=8), PARAMETER :: nfun_nag = 3 ! Ax, Ay, Az
1781  INTEGER(KIND=8), PARAMETER :: lenwrk_nag = iwrk
1782  INTEGER(KIND=8) :: maxcls_nag,mincls_nag, subs, restar, wrklen, rulcls, wrksbs, n, m, funcls
1783  DOUBLE PRECISION :: absreq_nag, relreq_nag
1784  DOUBLE PRECISION :: wrkstr_nag(lenwrk_nag)
1785  DOUBLE PRECISION :: a_nag(ndim_nag), b_nag(ndim_nag), &
1786  finest_nag(nfun_nag), absest_nag(nfun_nag)
1787  DOUBLE PRECISION, ALLOCATABLE :: vrtwrk(:)
1788 
1789 #ifdef NAG
1790  EXTERNAL :: d01eaf
1791 
1792 #else
1793  EXTERNAL :: dcuhre
1794 
1795 #endif
1796  ! BEGIN SUBROUTINE
1797  IF (adapt_tol < 0) THEN
1798  ! Not implmented
1799  CALL vecpot_volint_dbl(x,y,z,ax,ay,az)
1800  RETURN
1801  END IF
1802  a_nag(1:3) = zero
1803  b_nag(1:2) = one*pi2
1804  b_nag(3) = one
1805  mincls_nag = min_cls
1806  maxcls_nag = iwrk
1807  absreq_nag = adapt_tol ! Talk to Stuart about proper values
1808  relreq_nag = adapt_rel ! Talk to Stuart about proper values
1809  finest_nag = zero
1810  absest_nag = zero
1811  x_nag = x
1812  y_nag = y
1813  z_nag = z
1814  adapt_rerun = .true.
1815  subs = 1
1816  restar = 0
1817  DO WHILE (adapt_rerun)
1818 
1819 #ifdef NAG
1820  CALL d01eaf(ndim_nag,a_nag,b_nag,mincls_nag,maxcls_nag,nfun_nag,funsub_nag_a3d,absreq_nag,&
1821  relreq_nag,lenwrk_nag,wrkstr_nag,finest_nag,absest_nag,istat)
1822  IF (istat == 1 .or. istat == 3) THEN
1823  maxcls_nag = maxcls_nag*10
1824  mincls_nag = -1
1825  restar = 1
1826  WRITE(6,*) '!!!!! WARNING Could not reach desired tollerance !!!!!'
1827  WRITE(6,*) ' AX = ',finest_nag(1),' +/- ',absest_nag(1)
1828  WRITE(6,*) ' AY = ',finest_nag(2),' +/- ',absest_nag(2)
1829  WRITE(6,*) ' AZ = ',finest_nag(3),' +/- ',absest_nag(3)
1830  WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
1831  ELSE IF (istat < 0) THEN
1832  ax = zero
1833  ay = zero
1834  az = zero
1835  adapt_rerun=.false.
1836  ELSE
1837  ax = finest_nag(1)
1838  ay = finest_nag(2)
1839  az = finest_nag(3)
1840  adapt_rerun=.false.
1841  END IF
1842 
1843 #else
1844  IF (.not.ALLOCATED(vrtwrk)) THEN
1845  wrklen = ((iwrk-ndim_nag)/(2*ndim_nag) + 1)*(2*ndim_nag+2*nfun_nag+2) + 17*nfun_nag + 256
1846  ALLOCATE(vrtwrk(wrklen),stat=istat)
1847  IF (istat .ne. 0) THEN
1848  WRITE(6,*) ' ALLOCATION ERROR IN: vecpot_volint_adapt_dbl'
1849  WRITE(6,*) ' VARIABLE: VRTWRK, SIZE: ',wrklen
1850  WRITE(6,*) ' ISTAT: ',istat
1851  RETURN
1852  END IF
1853  END IF
1854  CALL dcuhre(ndim_nag,nfun_nag,a_nag,b_nag,mincls_nag,maxcls_nag,funsub_nag_a3d,absreq_nag,&
1855  relreq_nag,0,wrklen,restar,finest_nag,absest_nag,funcls,istat,vrtwrk)
1856  !DEALLOCATE(vrtwrk)
1857  IF (istat == 1) THEN
1858  maxcls_nag = maxcls_nag*10
1859  mincls_nag = -1
1860  restar = 1
1861  ELSE IF (istat > 1) THEN
1862  ax = zero
1863  ay = zero
1864  az = zero
1865  adapt_rerun=.false.
1866  DEALLOCATE(vrtwrk)
1867  ELSE
1868  ax = finest_nag(1)
1869  ay = finest_nag(2)
1870  az = finest_nag(3)
1871  adapt_rerun=.false.
1872  DEALLOCATE(vrtwrk)
1873  END IF
1874 
1875 #endif
1876  END DO
1877  nlastcall=mincls_nag
1878  RETURN
1879  ! END SUBROUTINE
1880  END SUBROUTINE vecpot_volint_adapt_dbl
1881  !-----------------------------------------------------------------
1882 
1883 
1884  !-----------------------------------------------------------------
1885  SUBROUTINE bfield_volint_adapt_flt(x_flt,y_flt,z_flt,bx_flt,by_flt,bz_flt,istat)
1886  IMPLICIT NONE
1887  ! INPUT VARIABLES
1888  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1889  REAL, INTENT(out) :: bx_flt, by_flt, bz_flt
1890  INTEGER, INTENT(inout) :: istat
1891  ! LOCAL VARIABLES
1892  DOUBLE PRECISION :: xt,yt,zt,bxt,byt,bzt
1893  ! BEGIN SUBROUTINE
1894  xt = x_flt
1895  yt = y_flt
1896  zt = z_flt
1897  bxt = zero
1898  byt = zero
1899  bzt = zero
1900  CALL bfield_volint_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1901  bx_flt = bxt
1902  by_flt = byt
1903  bz_flt = bzt
1904  RETURN
1905  ! END SUBROUTINE
1906  END SUBROUTINE bfield_volint_adapt_flt
1907  !-----------------------------------------------------------------
1908 
1909 
1910  !-----------------------------------------------------------------
1911  SUBROUTINE bfield_vc_flt(x_flt,y_flt,z_flt,bx_flt,by_flt,bz_flt,istat)
1912  IMPLICIT NONE
1913  ! INPUT VARIABLES
1914  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1915  REAL, INTENT(out) :: bx_flt, by_flt, bz_flt
1916  INTEGER, INTENT(inout) :: istat
1917  ! LOCAL VARIABLES
1918  DOUBLE PRECISION :: xt,yt,zt,bxt,byt,bzt
1919  ! BEGIN SUBROUTINE
1920  xt = x_flt
1921  yt = y_flt
1922  zt = z_flt
1923  bxt = zero
1924  byt = zero
1925  bzt = zero
1926  SELECT CASE(trim(vc_type_str))
1927  CASE('Surface Current')
1928  CALL bfield_virtual_casing_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1929  CASE('Volume Integral')
1930  CALL bfield_volint_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1931  END SELECT
1932  bx_flt = bxt
1933  by_flt = byt
1934  bz_flt = bzt
1935  RETURN
1936  ! END SUBROUTINE
1937  END SUBROUTINE bfield_vc_flt
1938  !-----------------------------------------------------------------
1939 
1940 
1941  !-----------------------------------------------------------------
1942  SUBROUTINE vecpot_vc_flt(x_flt,y_flt,z_flt,bx_flt,by_flt,bz_flt,istat)
1943  IMPLICIT NONE
1944  ! INPUT VARIABLES
1945  REAL, INTENT(in) :: x_flt, y_flt, z_flt
1946  REAL, INTENT(out) :: bx_flt, by_flt, bz_flt
1947  INTEGER, INTENT(inout) :: istat
1948  ! LOCAL VARIABLES
1949  DOUBLE PRECISION :: xt,yt,zt,bxt,byt,bzt
1950  ! BEGIN SUBROUTINE
1951  xt = x_flt
1952  yt = y_flt
1953  zt = z_flt
1954  bxt = zero
1955  byt = zero
1956  bzt = zero
1957  SELECT CASE(trim(vc_type_str))
1958  CASE('Surface Current')
1959  CALL vecpot_virtual_casing_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1960  CASE('Volume Integral')
1961  CALL vecpot_volint_adapt_dbl(xt,yt,zt,bxt,byt,bzt,istat)
1962  END SELECT
1963  bx_flt = bxt
1964  by_flt = byt
1965  bz_flt = bzt
1966  RETURN
1967  ! END SUBROUTINE
1968  END SUBROUTINE vecpot_vc_flt
1969  !-----------------------------------------------------------------
1970 
1971 
1972  !-----------------------------------------------------------------
1973  SUBROUTINE bfield_vc_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
1974  IMPLICIT NONE
1975  ! INPUT VARIABLES
1976  DOUBLE PRECISION, INTENT(in) :: x_dbl, y_dbl, z_dbl
1977  DOUBLE PRECISION, INTENT(out) :: bx_dbl, by_dbl, bz_dbl
1978  INTEGER, INTENT(inout) :: istat
1979  ! LOCAL VARIABLES
1980  ! BEGIN SUBROUTINE
1981  SELECT CASE(trim(vc_type_str))
1982  CASE('Surface Current')
1983  CALL bfield_virtual_casing_adapt_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
1984  CASE('Volume Integral')
1985  CALL bfield_volint_adapt_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
1986  END SELECT
1987  RETURN
1988  ! END SUBROUTINE
1989  END SUBROUTINE bfield_vc_dbl
1990  !-----------------------------------------------------------------
1991 
1992 
1993  !-----------------------------------------------------------------
1994  SUBROUTINE vecpot_vc_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
1995  IMPLICIT NONE
1996  ! INPUT VARIABLES
1997  DOUBLE PRECISION, INTENT(in) :: x_dbl, y_dbl, z_dbl
1998  DOUBLE PRECISION, INTENT(out) :: bx_dbl, by_dbl, bz_dbl
1999  INTEGER, INTENT(inout) :: istat
2000  ! LOCAL VARIABLES
2001  ! BEGIN SUBROUTINE
2002  SELECT CASE(trim(vc_type_str))
2003  CASE('Surface Current')
2004  CALL vecpot_virtual_casing_adapt_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
2005  CASE('Volume Integral')
2006  CALL vecpot_volint_adapt_dbl(x_dbl,y_dbl,z_dbl,bx_dbl,by_dbl,bz_dbl,istat)
2007  END SELECT
2008  RETURN
2009  ! END SUBROUTINE
2010  END SUBROUTINE vecpot_vc_dbl
2011  !-----------------------------------------------------------------
2012 
2013 
2014  !-----------------------------------------------------------------
2015  SUBROUTINE bfield_volint_adapt_dbl(x,y,z,bx,by,bz,istat)
2016  IMPLICIT NONE
2017  ! INPUT VARIABLES
2018  DOUBLE PRECISION, INTENT(in) :: x, y, z
2019  DOUBLE PRECISION, INTENT(out) :: bx, by, bz
2020  INTEGER, INTENT(inout) :: istat
2021  ! LOCAL VARIABLES
2022  LOGICAL :: adapt_rerun
2023  INTEGER(KIND=8), PARAMETER :: ndim_nag = 3 ! theta,zeta,s
2024  INTEGER(KIND=8), PARAMETER :: nfun_nag = 3 ! Bx, By, Bz
2025  INTEGER(KIND=8), PARAMETER :: lenwrk_nag = iwrk
2026  INTEGER(KIND=8) :: maxcls_nag,mincls_nag, subs, restar, wrklen, rulcls, wrksbs, n, m, funcls
2027  DOUBLE PRECISION :: absreq_nag, relreq_nag
2028  DOUBLE PRECISION :: wrkstr_nag(lenwrk_nag)
2029  DOUBLE PRECISION :: a_nag(ndim_nag), b_nag(ndim_nag), &
2030  finest_nag(nfun_nag), absest_nag(nfun_nag)
2031  DOUBLE PRECISION, ALLOCATABLE :: vrtwrk(:)
2032 
2033 #ifdef NAG
2034  EXTERNAL :: d01eaf
2035 
2036 #else
2037  EXTERNAL :: dcuhre
2038 
2039 #endif
2040  ! BEGIN SUBROUTINE
2041  IF (adapt_tol < 0) THEN
2042  ! Not implmented
2043  CALL bfield_volint_dbl(x,y,z,bx,by,bz)
2044  RETURN
2045  END IF
2046  a_nag(1:3) = zero
2047  b_nag(1:2) = one*pi2
2048  b_nag(3) = one
2049  mincls_nag = min_cls
2050  maxcls_nag = iwrk
2051  absreq_nag = adapt_tol ! Talk to Stuart about proper values
2052  relreq_nag = adapt_rel ! Talk to Stuart about proper values
2053  finest_nag = zero
2054  absest_nag = zero
2055  x_nag = x
2056  y_nag = y
2057  z_nag = z
2058  adapt_rerun = .true.
2059  subs = 1
2060  restar = 0
2061  DO WHILE (adapt_rerun)
2062 
2063 #ifdef NAG
2064  CALL d01eaf(ndim_nag,a_nag,b_nag,mincls_nag,maxcls_nag,nfun_nag,funsub_nag_b3d,absreq_nag,&
2065  relreq_nag,lenwrk_nag,wrkstr_nag,finest_nag,absest_nag,istat)
2066  IF (istat == 1 .or. istat == 3) THEN
2067  maxcls_nag = maxcls_nag*10
2068  mincls_nag = -1
2069  restar = 1
2070  WRITE(6,*) '!!!!! WARNING Could not reach desired tollerance !!!!!'
2071  WRITE(6,*) ' BX = ',finest_nag(1),' +/- ',absest_nag(1)
2072  WRITE(6,*) ' BY = ',finest_nag(2),' +/- ',absest_nag(2)
2073  WRITE(6,*) ' BZ = ',finest_nag(3),' +/- ',absest_nag(3)
2074  WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
2075  ELSE IF (istat < 0) THEN
2076  bx = zero
2077  by = zero
2078  bz = zero
2079  adapt_rerun=.false.
2080  ELSE
2081  bx = finest_nag(1)
2082  by = finest_nag(2)
2083  bz = finest_nag(3)
2084  adapt_rerun=.false.
2085  END IF
2086 
2087 #else
2088  IF (.not.ALLOCATED(vrtwrk)) THEN
2089  wrklen = ((iwrk-ndim_nag)/(2*ndim_nag) + 1)*(2*ndim_nag+2*nfun_nag+2) + 17*nfun_nag + 256
2090  ALLOCATE(vrtwrk(wrklen),stat=istat)
2091  IF (istat .ne. 0) THEN
2092  WRITE(6,*) ' ALLOCATION ERROR IN: bfield_volint_adapt_dbl'
2093  WRITE(6,*) ' VARIABLE: VRTWRK, SIZE: ',wrklen
2094  WRITE(6,*) ' ISTAT: ',istat
2095  RETURN
2096  END IF
2097  END IF
2098  CALL dcuhre(ndim_nag,nfun_nag,a_nag,b_nag,mincls_nag,maxcls_nag,funsub_nag_b3d,absreq_nag,&
2099  relreq_nag,0,wrklen,restar,finest_nag,absest_nag,funcls,istat,vrtwrk)
2100  !DEALLOCATE(vrtwrk)
2101  IF (istat == 1) THEN
2102  maxcls_nag = maxcls_nag*10
2103  mincls_nag = funcls
2104  restar = 1
2105  ELSE IF (istat > 1) THEN
2106  bx = zero
2107  by = zero
2108  bz = zero
2109  adapt_rerun=.false.
2110  DEALLOCATE(vrtwrk)
2111  ELSE
2112  bx = finest_nag(1)
2113  by = finest_nag(2)
2114  bz = finest_nag(3)
2115  adapt_rerun=.false.
2116  DEALLOCATE(vrtwrk)
2117  END IF
2118 
2119 #endif
2120  END DO
2121  nlastcall=mincls_nag
2122  RETURN
2123  ! END SUBROUTINE
2124  END SUBROUTINE bfield_volint_adapt_dbl
2125  !-----------------------------------------------------------------
2126 
2127 
2128  !-----------------------------------------------------------------
2129  SUBROUTINE mntouv_local(k1,k,mnmax,nu,nv,xu,xv,fmn,xm,xn,f,signs,calc_trig)
2130  IMPLICIT NONE
2131  ! INPUT VARIABLES
2132  INTEGER, INTENT(in) :: k1
2133  INTEGER, INTENT(in) :: k
2134  INTEGER, INTENT(in) :: mnmax
2135  INTEGER, INTENT(in) :: nu
2136  INTEGER, INTENT(in) :: nv
2137  DOUBLE PRECISION, INTENT(in) :: xu(1:nu)
2138  DOUBLE PRECISION, INTENT(in) :: xv(1:nv)
2139  DOUBLE PRECISION, INTENT(in) :: fmn(1:mnmax,k1:k)
2140  INTEGER, INTENT(in) :: xm(1:mnmax)
2141  INTEGER, INTENT(in) :: xn(1:mnmax)
2142  DOUBLE PRECISION, INTENT(inout) :: f(1:nu,1:nv,k1:k)
2143  INTEGER, INTENT(in) :: signs
2144  INTEGER, INTENT(in) :: calc_trig
2145  ! LOCAL VARIABLES
2146  INTEGER :: mn, i, ier, ik
2147  DOUBLE PRECISION :: xm_temp(1:mnmax,1)
2148  DOUBLE PRECISION :: xn_temp(1:mnmax,1)
2149  DOUBLE PRECISION :: pi2_l
2150  DOUBLE PRECISION :: mt(1:mnmax,1:nu)
2151  DOUBLE PRECISION :: nz(1:mnmax,1:nv)
2152  DOUBLE PRECISION :: fmn_temp(1:mnmax,1:nu)
2153  DOUBLE PRECISION :: xu_temp(1,1:nu)
2154  DOUBLE PRECISION :: xv_temp(1,1:nv)
2155  DOUBLE PRECISION :: fmn_help(1:mnmax)
2156  DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosmt(:,:)
2157  DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinmt(:,:)
2158  DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosnz(:,:)
2159  DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinnz(:,:)
2160  ! BEGIN SUBROUTINE
2161  pi2_l = 8.0d+0 * atan(1.)
2162  IF (calc_trig == 1) THEN
2163  IF (ALLOCATED(cosmt)) DEALLOCATE(cosmt)
2164  IF (ALLOCATED(sinmt)) DEALLOCATE(sinmt)
2165  IF (ALLOCATED(cosnz)) DEALLOCATE(cosnz)
2166  IF (ALLOCATED(sinnz)) DEALLOCATE(sinnz)
2167  ALLOCATE(cosmt(1:mnmax,1:nu),sinmt(1:mnmax,1:nu),&
2168  cosnz(1:mnmax,1:nv),sinnz(1:mnmax,1:nv),stat=ier)
2169  FORALL(i=1:mnmax) xm_temp(i,1)=dble(xm(i))
2170  FORALL(i=1:mnmax) xn_temp(i,1)=dble(xn(i))
2171  FORALL(i=1:nu) xu_temp(1,i)=xu(i)
2172  FORALL(i=1:nv) xv_temp(1,i)=xv(i)
2173  mt = matmul(xm_temp,xu_temp)
2174  nz = matmul(xn_temp,xv_temp)
2175  FORALL(mn=1:mnmax,i=1:nu) cosmt(mn,i) = dcos(pi2_l*mt(mn,i))
2176  FORALL(mn=1:mnmax,i=1:nu) sinmt(mn,i) = dsin(pi2_l*mt(mn,i))
2177  FORALL(mn=1:mnmax,i=1:nv) cosnz(mn,i) = dcos(pi2_l*nz(mn,i))
2178  FORALL(mn=1:mnmax,i=1:nv) sinnz(mn,i) = dsin(pi2_l*nz(mn,i))
2179  END IF
2180  IF (signs == 0) THEN
2181  DO ik = k1,k
2182  FORALL(mn=1:mnmax) fmn_help(mn)=fmn(mn,ik)
2183  fmn_temp=spread(fmn_help,2,nu)
2184  f(1:nu,1:nv,ik) = f(1:nu,1:nv,ik) + matmul(transpose((fmn_temp*cosmt)),cosnz) &
2185  - matmul(transpose((fmn_temp*sinmt)),sinnz)
2186  END DO
2187  ELSE IF (signs == 1) THEN
2188  DO ik = k1,k
2189  FORALL(mn=1:mnmax) fmn_help(mn)=fmn(mn,ik)
2190  fmn_temp=spread(fmn_help,2,nu)
2191  f(1:nu,1:nv,ik) = f(1:nu,1:nv,ik) + matmul(transpose((fmn_temp*sinmt)),cosnz) &
2192  + matmul(transpose((fmn_temp*cosmt)),sinnz)
2193  END DO
2194  END IF
2195  ! END SUBROUTINE
2196  END SUBROUTINE mntouv_local
2197  !-----------------------------------------------------------------
2198 
2199  !-----------------------------------------------------------------
2200  SUBROUTINE uvtomn_local(k1,k,mnmax,nu,nv,xu,xv,fmn,xm,xn,f,signs,calc_trig)
2201  IMPLICIT NONE
2202  ! INPUT VARIABLES
2203  INTEGER, INTENT(in) :: k1
2204  INTEGER, INTENT(in) :: k
2205  INTEGER, INTENT(in) :: mnmax
2206  INTEGER, INTENT(in) :: nu
2207  INTEGER, INTENT(in) :: nv
2208  DOUBLE PRECISION, INTENT(in) :: xu(1:nu)
2209  DOUBLE PRECISION, INTENT(in) :: xv(1:nv)
2210  DOUBLE PRECISION, INTENT(inout) :: fmn(1:mnmax,k1:k)
2211  INTEGER, INTENT(in) :: xm(1:mnmax)
2212  INTEGER, INTENT(in) :: xn(1:mnmax)
2213  DOUBLE PRECISION, INTENT(in) :: f(1:nu,1:nv,k1:k)
2214  INTEGER, INTENT(in) :: signs
2215  INTEGER, INTENT(in) :: calc_trig
2216  ! LOCAL VARIABLES
2217  INTEGER :: mn, i, j, ier, ik
2218  DOUBLE PRECISION :: xn_temp(1:mnmax,1)
2219  DOUBLE PRECISION :: xm_temp(1:mnmax,1)
2220  DOUBLE PRECISION :: pi2_l
2221  DOUBLE PRECISION :: mt(1:mnmax,1:nu)
2222  DOUBLE PRECISION :: nz(1:mnmax,1:nv)
2223  DOUBLE PRECISION :: xu_temp(1,1:nu)
2224  DOUBLE PRECISION :: xv_temp(1,1:nv)
2225  DOUBLE PRECISION, ALLOCATABLE, SAVE :: fnuv(:)
2226  DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosmt(:,:)
2227  DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinmt(:,:)
2228  DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosnz(:,:)
2229  DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinnz(:,:)
2230  ! BEGIN SUBROUTINE
2231  pi2_l = 8.0d+0 * atan(1.)
2232  IF (calc_trig == 1) THEN
2233  IF (ALLOCATED(cosmt)) DEALLOCATE(cosmt)
2234  IF (ALLOCATED(sinmt)) DEALLOCATE(sinmt)
2235  IF (ALLOCATED(cosnz)) DEALLOCATE(cosnz)
2236  IF (ALLOCATED(sinnz)) DEALLOCATE(sinnz)
2237  IF (ALLOCATED(fnuv)) DEALLOCATE(fnuv)
2238  ALLOCATE(fnuv(1:mnmax),stat=ier)
2239  ALLOCATE(cosmt(1:mnmax,1:nu),sinmt(1:mnmax,1:nu),&
2240  cosnz(1:mnmax,1:nv),sinnz(1:mnmax,1:nv),stat=ier)
2241  FORALL(i=1:mnmax) xm_temp(i,1)=dble(xm(i))
2242  FORALL(i=1:mnmax) xn_temp(i,1)=dble(xn(i))
2243  FORALL(i=1:mnmax) fnuv(i) = 2.0d+0/dble(nu*nv)
2244  WHERE(xm == 0) fnuv = 5.0d-1*fnuv
2245  FORALL(i=1:nu) xu_temp(1,i)=xu(i)
2246  FORALL(i=1:nv) xv_temp(1,i)=xv(i)
2247  mt = matmul(xm_temp,xu_temp)
2248  nz = matmul(xn_temp,xv_temp)
2249  FORALL(mn=1:mnmax,i=1:nu) cosmt(mn,i) = dcos(pi2_l*mt(mn,i))
2250  FORALL(mn=1:mnmax,i=1:nu) sinmt(mn,i) = dsin(pi2_l*mt(mn,i))
2251  FORALL(mn=1:mnmax,i=1:nv) cosnz(mn,i) = dcos(pi2_l*nz(mn,i))
2252  FORALL(mn=1:mnmax,i=1:nv) sinnz(mn,i) = dsin(pi2_l*nz(mn,i))
2253  END IF
2254  fmn = zero
2255  IF (signs == 0) THEN
2256  DO mn = 1, mnmax
2257  DO i = 1, nu
2258  DO j = 1, nv
2259  !PRINT *,cosmt(mn,i), cosnz(mn,j), sinmt(mn,i), sinnz(mn,j),fnuv(mn)
2260  fmn(mn,k1:k) = fmn(mn,k1:k) + f(i,j,k1:k)*(cosmt(mn,i)*cosnz(mn,j) &
2261  - sinmt(mn,i)*sinnz(mn,j))*fnuv(mn)
2262  END DO
2263  END DO
2264  END DO
2265  ELSE IF (signs == 1) THEN
2266  DO mn = 1, mnmax
2267  DO i = 1, nu
2268  DO j = 1, nv
2269  fmn(mn,k1:k) = fmn(mn,k1:k) + f(i,j,k1:k)*(sinmt(mn,i)*cosnz(mn,j) &
2270  + cosmt(mn,i)*sinnz(mn,j))*fnuv(mn)
2271  END DO
2272  END DO
2273  END DO
2274  END IF
2275  ! END SUBROUTINE
2276  END SUBROUTINE uvtomn_local
2277  !-----------------------------------------------------------------
2278 
2279 !-----------------------------------------------------------------------
2280 ! End Module
2281 !-----------------------------------------------------------------------
2282  END MODULE virtual_casing_mod
2283 
2284 
ezspline::ezspline_setup
Definition: ezspline.f90:216
ezspline_obj::ezspline_allocated
Definition: ezspline_obj.f90:19
ezspline::ezspline_init
Definition: ezspline.f90:3
ezspline::ezspline_interp
Definition: ezspline.f90:280
ezspline::ezspline_free
Definition: ezspline.f90:174
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17