text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine resultsmech_u1(co,kon,ipkon,lakon,ne,v,
     &  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,eme,iperturb,fn,iout,qa,vold,nmethod,
     &  veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc,
     &  reltime,calcul_fn,calcul_qa,calcul_cauchy,nener,
     &  ikin,nal,ne0,thicke,emeini,i,ielprop,prop)
!
!     calculates nal,qa,fn,xstiff,ener,eme,eei,stx for user element 1
!
!     This is a beam type element. Reference:
!     Yunhua Luo, An Efficient 3D Timoshenko Beam Element with
!     Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1,
!     2008, no. 3, 95-106
!
!
!     special case for which the beam axis goes through the
!     center of gravity of the cross section and the 1-direction
!     corresponds with a principal axis
!
!     note that the strain components are Lagrange strain components,
!     the stress components are Piola-Kirchhoff components of the
!     second kind. For linear geometric calculations the strains
!     reduce to the infinitesimal strains
!
!
!     INPUT:
!
!     co(1..3,i)         coordinates of node i
!     kon(*)             contains the topology of all elements. The
!                        topology of element i starts at kon(ipkon(i)+1)
!                        and continues until all nodes are covered. The
!                        number of nodes depends on the element label
!     ipkon(i)           points to the location in field kon preceding
!                        the topology of element i
!     lakon(i)           label of element i (character*8)
!     ne                 highest element number in the mesh
!     v(0..mi(2),i)      value of the variables in node i at the end
!                        of the present iteration
!     elcon              elastic constants (cf. List of variables and
!                        their meaning in the User's Manual)
!     nelcon             integers describing the elastic constant fields
!                        (cf. User's Manual)
!     rhcon              density constants (cf. User's Manual)
!     nrhcon             integers describing the density constant fields
!                        (cf. User's Manual)
!     alcon              thermal expansion constants (cf. User's Manual)
!     nalcon             integers describing the thermal expansion
!                        constants (cf. User's Manual)
!     alzero             thermal expansion reference values (cf. User's Manual)
!     ielmat(i)          material for element i
!     ielorien(i)        orientation for element i
!     norien             number of orientations
!     orab(7,*)          description of all local coordinate systems.
!                        (cf. List of variables and their meaning in the
!                        User's manual)
!     ntmat_             maximum number of material temperature points
!     t0(i)              temperature in node i at start of calculation
!     t1(i)              temperature in node i at the end of the current
!                        increment
!     ithermal(1..2)     cf. List of variables and
!                        their meaning in the User's Manual
!     prestr(i,j,k)      residual stress component i in integration point j
!                        of element k
!     iprestr            if 0: no residual stresses
!                        else: residual stresses
!     iperturb(*)        describes the kind of nonlinearity of the
!                        calculation, cf. User's Manual
!     iout               if -2: v is assumed to be known and is used to
!                               calculate strains, stresses..., no result output
!                               corresponds to iout=-1 with in addition the
!                               calculation of the internal energy density
!                        if -1: v is assumed to be known and is used to
!                               calculate strains, stresses..., no result
!                               output; is used to take changes in SPC's and
!                               MPC's at the start of a new increment or
!                               iteration into account
!                        if 0: v is calculated from the system solution
!                              and strains, stresses.. are calculated,
!                              no result output
!                        if 1: v is calculated from the system solution and
!                              strains, stresses.. are calculated, requested
!                              results output
!                        if 2: v is assumed to be known and is used to
!                              calculate strains, stresses..., requested
!                              results output
!     vold(0..mi(2),i)   value of the variables in node i at the end
!                        of the previous iteration
!     nmethod            procedure:
!                        1: static analysis
!                        2: frequency analysis
!                        3: buckling analysis
!                        4: (linear or nonlinear) dynamic analysis
!                        5: steady state dynamics analysis
!                        6: Coriolis frequency calculation
!                        7: flutter frequency calculation
!                        8:  magnetostatics
!                        9:  magnetodynamics
!                        10: electromagnetic eigenvalue problems
!                        11: superelement creation or Green function
!                            calculation
!                        12: sensitivity analysis
!     veold(j,i)         time rate of variable j in node i at the end
!                        of the previous iteration
!     dtime              length of present time increment
!     time               step time at the end of the present increment
!     ttime              total time at the start of the present increment
!     plicon,nplicon     fields describing isotropic hardening of
!                        a plastic material or spring constants of
!                        a nonlinear spring (cf. User's Manual)
!     plkcon,nplkcon     fields describing kinematic hardening of
!                        a plastic material or gap conductance
!                        constants (cf. User's Manual)
!     xstateini(i,j,k)   state variable component i in integration point j
!                        of element k at the start of the present increment
!     xstate(i,j,k)      state variable component i in integration point j
!                        of element k at the end of the present increment
!     npmat_             maximum number of plastic constants
!     matname(i)         name of material i
!     mi(1)              max # of integration points per element (max
!                        over all elements)
!     mi(2)              max degree of freedom per node (max over all
!                        nodes) in fields like v(0:mi(2))...
!     mi(3)              max number of layers in the structure
!     ielas              0: no elastic iteration: irreversible effects
!                        are allowed
!                        1: elastic iteration, i.e. no irreversible
!                           deformation allowed
!     icmd               not equal to 3: calculate stress and stiffness
!                        3: calculate only stress
!     ncmat_             max number of elastic constants
!     nstate_            max number of state variables in any integration
!                        point
!     stiini(i,j,k)      stress component i in integration point j
!                        of element k at the start of the present
!                        increment (= end of last increment)
!     vini(0..mi(2),i)   value of the variables in node i at the start
!                        of the present increment
!     enerini(j,k)       internal energy density at integration point j
!                        of element k at the start of the present increment
!     istep              current step number
!     iinc               current increment number within the actual step
!     reltime            relative step time (between 0 and 1)
!     calcul_fn          if 0: no nodal forces have to be calculated
!                        else: nodal forces are required on output
!     calcul_qa          if 0: no mean forces have to be calculated
!                        else: mean forces are required on output
!     calcul_cauchy      if 0: no Cauchy stresses are required
!                        else: Cauchy stresses are required on output: have
!                              to be calculated from the PK2 stresses
!     nener              if 0: internal energy calculation is not required
!                        else: internal energy is required on output
!     ikin               if 0: kinetic energy calculation is not requred
!                        else: kinetic energy is required on output
!     ne0                largest element number without contact elements (are
!                        stored after all other elements)
!     thicke(j,i)        layer thickness for layer j in element i
!     emeini(i,j,k)      mechanical strain component i in integration point j
!                        of element k at the start of the present increment
!     i                  actual element at stake
!     ielprop(i)         points to the location in field prop preceding
!                        the properties of element i
!     prop(*)            contains the properties and some beam
!                        elements (cf. User's Manual)
!
!
!     OUTPUT:
!
!     stx(i,j,k)         stress component i in integration point j
!                        of element k at the end of the present
!                        iteration
!     eme(i,j,k)         mechanical strain component i in integration point j
!                        of element k at the end of the present iteration
!     fn(j,i)            nodal force component j in node i
!     qa(1..4)           qa(1): average force
!                        qa(2): average flux
!                        ... cf. User's Manual
!     xstiff(i,j,k)      stiffness (i=1...21) and conductivity
!                        (i=22..27) constants in integration point j
!                        of element k
!     ener(j,k)          internal energy density at integration point j
!                        of element k at the end of the present increment
!     eei(i,j,k)         total strain component i in integration point j
!                        of element k at the end of the present iteration
!     nal                number of nodal force contributions
!
      implicit none
!
      character*8 lakon(*)
      character*80 amat,matname(*)
!
      integer kon(*),konl(26),mi(*),
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
     &  ielorien(mi(3),*),ntmat_,ipkon(*),ne0,
     &  istep,iinc,ne,mattyp,ithermal(2),iprestr,i,j,k,m1,m2,jj,
     &  i1,kk,nener,indexe,nope,norien,iperturb(*),iout,
     &  nal,icmd,ihyper,nmethod,kode,imat,iorien,ielas,
     &  istiff,ncmat_,nstate_,ikin,ielprop(*),
     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn,
     &  calcul_cauchy,calcul_qa,index,node
!
      real*8 co(3,*),v(0:mi(2),*),stiini(6,mi(1),*),
     &  stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),
     &  alcon(0:6,ntmat_,*),vini(0:mi(2),*),
     &  alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*),
     &  q(0:mi(2),26),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*),
     &  vold(0:mi(2),*),eloc(9),elconloc(21),eth(6),coords(3),
     &  ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*),
     &  veold(0:mi(2),*),e,un,um,tt,dl,qa(*),t0l,t1l,dtime,time,ttime,
     &  plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
     &  xstateini(nstate_,mi(1),*),tm(3,3),a,reltime,
     &  thicke(mi(3),*),emeini(6,mi(1),*),aly,alz,bey,bez,xi(2),
     &  vlp(6,2),xi11,xi12,xi22,xk,offset1,offset2,e1(3),e2(3),e3(3)
!
      intent(in) co,kon,ipkon,lakon,ne,v,
     &  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,iperturb,iout,vold,nmethod,
     &  veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstate,npmat_,matname,mi,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,enerini,istep,iinc,
     &  reltime,calcul_fn,calcul_qa,calcul_cauchy,nener,
     &  ikin,ne0,thicke,emeini,i,ielprop,prop
!
      intent(inout) nal,qa,fn,xstiff,ener,eme,eei,stx
!
      nope=2
!
      indexe=ipkon(i)
!
!     material and orientation
!
      imat=ielmat(1,i)
      amat=matname(imat)
      if(norien.gt.0) then
         iorien=ielorien(1,i)
      else
         iorien=0
      endif
      if(iorien.gt.0) then
         write(*,*) '*ERROR in resultsmech_u1: no orientation'
         write(*,*) '       calculation for this type of element'
         call exit(201)
      endif
!
      if(nelcon(1,imat).lt.0) then
         ihyper=1
      else
         ihyper=0
      endif
!
!     properties of the cross section
!
      index=ielprop(i)
      a=prop(index+1)
      xi11=prop(index+2)
      xi12=prop(index+3)
      xi22=prop(index+4)
      xk=prop(index+5)
      e2(1)=prop(index+6)
      e2(2)=prop(index+7)
      e2(3)=prop(index+8)
      offset1=prop(index+9)
      offset2=prop(index+10)
!
      if(dabs(xi12).gt.0.d0) then
         write(*,*) '*ERROR in resultsmech_u1: no nonzero cross moment'
         write(*,*) '       of inertia for this type of element'
         call exit(201)
      endif
!
      if(dabs(offset1).gt.0.d0) then
         write(*,*) '*ERROR in resultsmech_u1: no offset in 1-direction'
         write(*,*) '       for this type of element'
         call exit(201)
      endif
!
      if(dabs(offset2).gt.0.d0) then
         write(*,*) '*ERROR in resultsmech_u1: no offset in 2-direction'
         write(*,*) '       for this type of element'
         call exit(201)
      endif
!
!     computation of the coordinates and displacements of the local nodes
!
      do k=1,nope
         konl(k)=kon(indexe+k)
         do j=1,3
            xl(j,k)=co(j,konl(k))
            vl(j,k)=v(j,konl(k))
         enddo
      enddo
!
!     q contains the nodal forces per element; initialization of q
!
      if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1)))
     &     then
         do m1=1,nope
            do m2=0,mi(2)
               q(m2,m1)=fn(m2,konl(m1))
            enddo
         enddo
      endif
!
!     local axes e1-e2-e3  (e2 is given by the user (in the input deck
!     this is called e1), e1 is parallel to the beam axis)
!
      do j=1,3
         e1(j)=xl(j,2)-xl(j,1)
      enddo
!
      dl=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3))
!
      do j=1,3
         e1(j)=e1(j)/dl
      enddo
!
!     e3 = e1 x e2
!
      e3(1)=e1(2)*e2(3)-e1(3)*e2(2)
      e3(2)=e1(3)*e2(1)-e1(1)*e2(3)
      e3(3)=e1(1)*e2(2)-e1(2)*e2(1)
!
!     transformation matrix from the global into the local system
!
      do j=1,3
         tm(1,j)=e1(j)
         tm(2,j)=e2(j)
         tm(3,j)=e3(j)
      enddo
!
!     calculating the temperature in the integration
!     point
!
      t0l=0.d0
      t1l=0.d0
      if(ithermal(1).eq.1) then
         do i1=1,nope
            t0l=t0l+t0(konl(i1))/2.d0
            t1l=t1l+t1(konl(i1))/2.d0
         enddo
      elseif(ithermal(1).ge.2) then
         write(*,*) '*ERROR in resultsmech_u1: no thermal'
         write(*,*) '       calculation for this type of element'
         call exit(201)
      endif
      tt=t1l-t0l
!
      kode=nelcon(1,imat)
!
      if(kode.eq.2) then
         mattyp=1
      else
         mattyp=0
      endif
!
!     material data and local stiffness matrix
!
      istiff=0
      call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
     &     imat,amat,iorien,coords,orab,ntmat_,elas,rho,
     &     i,ithermal,alzero,mattyp,t0l,t1l,
     &     ihyper,istiff,elconloc,eth,kode,plicon,
     &     nplicon,plkcon,nplkcon,npmat_,
     &     plconloc,mi(1),dtime,kk,
     &     xstiff,ncmat_)
!
!     correcting the thermal strains
!
      do k=2,6
         eth(k)=0.d0
      enddo
!
      if(mattyp.eq.1) then
         e=elconloc(1)
         un=elconloc(2)
         um=e/(1.d0+un)
         um=um/2.d0
      else
         write(*,*) '*ERROR in resultsmech_u1: no anisotropic material'
         write(*,*) '       calculation for this type of element'
         call exit(201)
      endif
!
!     calculating the local displacement and rotation values
!
!     values 1..6 of vl are u,v,w,phi,psi,theta from Luo
!
      do k=1,nope
         node=kon(indexe+k)
         do j=1,3
            vl(j,k)=tm(j,1)*v(1,node)
     &             +tm(j,2)*v(2,node)
     &             +tm(j,3)*v(3,node)
            vl(j+3,k)=tm(j,1)*v(4,node)
     &               +tm(j,2)*v(5,node)
     &               +tm(j,3)*v(6,node)
         enddo
      enddo
!
      xi(1)=0.d0
      xi(2)=1.d0
!
      aly=12.d0*e*xi11/(xk*um*a*dl*dl)
      alz=12.d0*e*xi22/(xk*um*a*dl*dl)
      bey=1.d0/(1.d0-aly)
      bez=1.d0/(1.d0-alz)
!
      do jj=1,nope
!
!        calculating the derivative of the local displacement and
!        rotation values (from Eqn. (19) and (20) in Luo)
!
         vlp(1,jj)=(-vl(1,1)+vl(1,2))/dl
         vlp(2,jj)=
     &     bey*(6.d0*xi(jj)**2-6.d0*xi(jj)+aly*xi(jj))/dl*vl(2,1)+
     &    bey*(3.d0*xi(jj)**2+(aly-4.d0)*xi(jj)+(1.d0-aly/2.d0))*vl(6,1)
     &    +bey*(-6.d0*xi(jj)**2+6.d0*xi(jj)-aly)/dl*vl(2,2)
     &    +bey*(3.d0*xi(jj)**2-(2.d0+aly)*xi(jj)+aly/2.d0)*vl(6,2)
         vlp(3,jj)=bez*(6.d0*xi(jj)*xi(jj)-6.d0*xi(jj)+alz)/dl*vl(3,1)+
     &    bez*(3.d0*xi(jj)**2+(alz-4.d0)*xi(jj)+(1.d0-alz/2.d0))*vl(5,1)
     &    +bez*(-6.d0*xi(jj)**2+6.d0*xi(jj)-alz)/dl*vl(3,2)
     &    +bez*(3.d0*xi(jj)**2-(2.d0+alz)*xi(jj)+alz/2.d0)*vl(5,2)
         vlp(4,jj)=(vl(4,2)-vl(4,1))/dl
         vlp(5,jj)=6.d0*bez*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(3,1)
     &           +bez*(6.d0*xi(jj)+(alz-4.d0))/dl*vl(5,1)
     &           +6.d0*bez*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(3,2)
     &           +bez*(6.d0*xi(jj)-(alz+2))/dl*vl(5,2)
         vlp(6,jj)=6.d0*bey*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(2,1)
     &           +bey*(6.d0*xi(jj)+(aly-4.d0))/dl*vl(6,1)
     &           +6.d0*bey*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(2,2)
     &           +bey*(6.d0*xi(jj)-(aly+2.d0))/dl*vl(6,2)
!
!        calculation of the strains (Eqn. (8) in Luo)
!
         eloc(1)=vlp(1,jj)
         eloc(2)=-vlp(6,jj)
         eloc(3)=vlp(5,jj)
         eloc(4)=vlp(2,jj)-vl(6,jj)
         eloc(5)=vlp(3,jj)+vl(5,jj)
         eloc(6)=vlp(4,jj)
!
!        determining the mechanical strain
!
         if(ithermal(1).ne.0) then
            do m1=2,6
               emec(m1)=eloc(m1)-eth(m1)
            enddo
         else
            do m1=1,6
               emec(m1)=eloc(m1)
            enddo
         endif
!
!        subtracting initial strains
!
         if(iprestr.ne.0) then
            write(*,*) '*ERROR in resultsmech_u1:'
            write(*,*) '       no initial strains allowed for'
            write(*,*) '       this user element'
            call exit(201)
         endif
!
!        calculating the section forces
!        simplified version of Eqn. (11) in Luo (symmetric case
!        for which Ay=Az=J=0)
!
         stre(1)=e*a*emec(1)
         stre(2)=e*xi11*emec(2)
         stre(3)=e*xi22*emec(3)
         stre(4)=xk*um*a*emec(4)
         stre(5)=xk*um*a*emec(5)
         stre(6)=xk*um*(xi11+xi22)*emec(6)
!
!        updating the internal energy and mechanical strain
!
         if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100).or.
     &        ((nmethod.eq.4).and.(iperturb(1).gt.1).and.
     &        (ithermal(1).le.1))) then
            if(ithermal(1).eq.0) then
               do m1=1,6
                  eth(m1)=0.d0
               enddo
            endif
            if(nener.eq.1) then
               ener(jj,i)=enerini(jj,i)+
     &              ((eloc(1)-eth(1)-emeini(1,jj,i))*
     &              (stre(1)+stiini(1,jj,i))+
     &              (eloc(2)-eth(2)-emeini(2,jj,i))*
     &              (stre(2)+stiini(2,jj,i))+
     &              (eloc(3)-eth(3)-emeini(3,jj,i))*
     &              (stre(3)+stiini(3,jj,i)))/2.d0+
     &         (eloc(4)-eth(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
     &         (eloc(5)-eth(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
     &         (eloc(6)-eth(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i))
            endif
            eme(1,jj,i)=eloc(1)-eth(1)
            eme(2,jj,i)=eloc(2)-eth(2)
            eme(3,jj,i)=eloc(3)-eth(3)
            eme(4,jj,i)=eloc(4)-eth(4)
            eme(5,jj,i)=eloc(5)-eth(5)
            eme(6,jj,i)=eloc(6)-eth(6)
         endif
!
         if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100)) then
!
            eei(1,jj,i)=eloc(1)
            eei(2,jj,i)=eloc(2)
            eei(3,jj,i)=eloc(3)
            eei(4,jj,i)=eloc(4)
            eei(5,jj,i)=eloc(5)
            eei(6,jj,i)=eloc(6)
         endif
!
!        updating the kinetic energy
!
         if(ikin.eq.1) then
            write(*,*) '*ERROR in resultsmech_u1:'
            write(*,*) '       no initial strains allowed for'
            write(*,*) '       this user element'
            call exit(201)
         endif
!
         stx(1,jj,i)=stre(1)
         stx(2,jj,i)=stre(2)
         stx(3,jj,i)=stre(3)
         stx(4,jj,i)=stre(4)
         stx(5,jj,i)=stre(5)
         stx(6,jj,i)=stre(6)
!
!     calculation of the global nodal forces
!
         if(calcul_fn.eq.1)then
            node=kon(indexe+jj)
            do j=1,3
               fn(j,node)=fn(j,node)+tm(1,j)*stre(1)
     &              +tm(2,j)*stre(2)
     &              +tm(3,j)*stre(3)
               fn(j+3,node)=fn(j+3,node)+tm(1,j)*stre(4)
     &              +tm(2,j)*stre(5)
     &              +tm(3,j)*stre(6)
               enddo
         endif
      enddo
!
!     q contains the contributions to the nodal force in the nodes
!     belonging to the element at stake from other elements (elements
!     already treated). These contributions have to be
!     subtracted to get the contributions attributable to the element
!     at stake only
!
      if(calcul_qa.eq.1) then
         do m1=1,nope
            do m2=1,3
               qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1))
            enddo
         enddo
         nal=nal+3*nope
      endif
!
      return
      end


text:
C=====================================================================72
      subroutine gen_vecfield_file(d,nx,ny,xmin,xmax,ymin,ymax,ix1,iy1)
      implicit none
      integer nx,ny, ix1,iy1, ix,iy, ic, i, j, nstrokes, ndel,ixoff,iix
      real d(nx,ny,2),xmin,xmax,ymin,ymax,scale, stroke(4,10000)
      real eps, dx, dy, x0, y0, x, y, signx
      real percentile_value   ! (d,n,p)
      real v2, vv, v80
      allocatable v2(:)

      x0 = xmin
      y0 = ymin
      if(ix1 .lt. 1) x0 = 0.0
      if(iy1 .lt. 1) y0 = 0.0

      dx = (xmax-x0)/float(nx)
      dy = (ymax-y0)/float(ny)
      ndel = max(2, int(max(nx,ny)/30))

      ! Calculate scales & sizes
      allocate (v2(nx*ny))
      ic = 0
      do iy = 1, ny
      do ix = 1, nx
        vv = d(ix,iy,1)**2 + d(ix,iy,2)**2
        if(vv .gt. 0.0) then
          ic = ic + 1
          v2(ic) = vv
        endif
      enddo
      enddo
      v80 = sqrt(percentile_value(v2,ic,80.0))
      deallocate (v2)

      eps = v80 / 1000000.0
      scale = 2.0 * float(ndel) * min(dx,dy) / v80

      ! Generate & mirror vectors as needed
      nstrokes = 0
      ixoff = 1
      do iy=ndel/2,ny,ndel
      ixoff = ixoff + max(1, ndel/5)
      if(ixoff .gt. ndel) ixoff = 1

      do iix=ix1+ixoff,nx,ndel
        ix = iix
        signx = scale
        if(iix .lt. 1) then
          ix = 1-iix
          signx  = -1.0 * scale
        endif
        if(sqrt(d(ix,iy,1)**2+d(ix,iy,2)**2) .ge. eps) then
          nstrokes = nstrokes + 1
          stroke(1,nstrokes) = x0 + dx*(0.5+float(iix-1))
          stroke(2,nstrokes) = y0 + dy*(0.5+float(iy -1))
          stroke(3,nstrokes) = signx * d(ix,iy,1)
          stroke(4,nstrokes) = scale * d(ix,iy,2)
        endif
      enddo
      enddo
c      if(ix1 .lt. 1) call mirror_strokes(nstrokes,stroke,-1.0,+1.0)
c      if(iy1 .lt. 1) call mirror_strokes(nstrokes,stroke,+1.0,-1.0)

      ! Write vectors to file
      open(unit=23,file="vecfield.vec",form="formatted")
      do i = 1, nstrokes
        write (23,999) (stroke(j,i), j=1,4)
999     format(1p4e15.6)
      enddo
      close(23)

      return
      end

C=====================================================================72
      ! Returns the value of the p-th percentile in d(1:n)
      real function percentile_value(d,n,p)
      implicit none
      real d(n), p, dbin, val, aval, pval
      integer n, i, mbin, ibin, icp, ic, ibinp
      parameter(mbin=10000)
      integer idist(-mbin:mbin)

      dbin = 0.01
      do ibin=-mbin,mbin
        idist(ibin) = 0
      enddo
      ic = 0

      do i = 1, n
        val = d(i)
        if(val .eq. 0) then
          idist(0) = idist(0) + 1
        else
          aval = abs(val)
          ibin = int(float(mbin)+alog(abs(aval))/dbin) - mbin
          ibin = max(-mbin/2,min(mbin/2,ibin))
          if(val .gt. 0.0) ibin = ibin + mbin/2
          if(val .lt. 0.0) ibin = ibin - mbin/2
          idist(ibin) = idist(ibin) + 1
        endif
        ic = ic + 1
      enddo

      icp = int(0.01 * p * float(ic))
      ibinp = -mbin
      ic = 0
      do ibin=-mbin,mbin
       if(ic .lt. icp) ibinp = ibin
       ic = ic + idist(ibin)
      enddo

      if(ibinp .eq. 0) pval = 0.0
      if(ibinp .gt. 0) pval = exp( dbin*float(ibinp - mbin/2))
      if(ibinp .lt. 0) pval = exp(-dbin*float(ibinp + mbin/2))

      percentile_value = pval
      return
      end

C=====================================================================72
      subroutine vecplot(d,nx,ny,xmin,xmax,ymin,ymax,coutfile,
     1                   cx,cy, vmin,vmax,ix1,ix2,iy1,iy2)
      implicit none
      character*(*) coutfile, cx, cy
      character*64  cpltfile
      integer nx,ny, ix1,ix2,iy1,iy2
      real d(nx,ny,2),xmin,xmax,ymin,ymax, vmin, vmax, v1, v2,scale,ex
      integer ix,iy, i,j,n,m, imap, nx1, ny1, ixx, iyy
      real  xtic, auto_tic, x1, x2, y1, y2, fmap, vx1,vx2,vy1,vy2
      integer nstrokes
      real stroke(4,10000)

      vx1 = vmin
      vx2 = vmax
      call autoscale(d(1,1,1),nx,ny,vx1,vx2,imap,scale)
      vy1 = vmin
      vy2 = vmax
      call autoscale(d(2,1,1),nx,ny,vy1,vy2,imap,scale)
      v1 = 0.0
      v2 = max(vx2,vy2)

      ! Generate gnuplot script
      cpltfile( 1:32) = "                                "
      cpltfile(33:64) = "                                "
      n = len(trim(coutfile))
      cpltfile(1:n) = trim(coutfile)
      cpltfile(n+1:n+6) = ".plt  "
      open(unit=23,file=cpltfile,form="formatted")
999   format(a, a, a, a, a, a)
998   format(a, " [", e15.6,":",e15.6,"]")
      write (23,996) "set size ratio", (y2-y1)/(x2-x1)
#ifndef ISGFORTRAN
996   format(a, f, f)
#else
996   format(a, f6.2, f6.2)
#endif
      if((y2-y1)/(x2-x1) .gt. 2.0) then
        xtic = auto_tic(x1, x2)
        write (23,996) "set xtic", xtic, xtic
      endif

      write (23,999) "set xlabel '", trim(cx), "'"
      write (23,999) "set ylabel '", trim(cy), "'"
      write (23,998) "set xrange", x1, x2
      write (23,998) "set yrange", y1, y2

      call get_vec_geom(nstrokes,stroke)
      if(nstrokes .le. 0) then
        write (23,999) "splot 'vecfield.vec' using 1:2:3:4"
      else
        if(ix1 .lt. 1) call mirror_strokes(nstrokes,stroke,-1.0,+1.0)
        if(iy1 .lt. 1) call mirror_strokes(nstrokes,stroke,+1.0,-1.0)
        call myclip(nstrokes,stroke,x1,x2,y1,y2)
        open(unit=37,file="geometry.4vec",form="formatted")
        do i=1,nstrokes
          write (37,*) (stroke(j,i),j=1,4)
        enddo
        close(37)
        write (23,999) "set clip two"
        write (23,999) "set style arrow 7 nohead ls 1"
#ifndef ISGFORTRAN
        write (23,999) "splot 'vecfield.vec using' 1:2:3:4"
     1                 ", 'geometry.4vec' using 1:2:3:4",
     1                 " notitle with vectors arrowstyle 7"
#else
        write (23,999) "splot 'vecfield.vec using' 1:2:3:4",
     1                 ", 'geometry.4vec' using 1:2:3:4",
     1                 " notitle with vectors arrowstyle 7"
#endif

      endif
      close(23)
      return
      end

C=====================================================================72

ext: f

max_stars_repo_name: nannau/astrohubs

max_issues_repo_name: nannau/astrohubs

max_forks_repo_name: nannau/astrohubs


!***********************************/
!	Name:数値実験のまとめ
!	Alias:PreNumericalEx
!	Description:
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:
!	Date:
!	Update:
!	Other:
!***********************************/
subroutine UOutputExact(UConf,UG,NE,iStep)
    use StructVar_Mod
    use LoopVar_Mod
    use ExactSol
    implicit none
    type(Configulation), intent(in) :: UConf
    type(UnstructuredGrid), intent(in) :: UG
    type(NumericalExperiment), intent(in) :: NE
    integer, intent(in) :: iStep
    character(len=64) :: cDirectory,cFileName, cCaseName
    character(len=32) :: cStep
    double precision :: Width
    integer :: iDir

    Width = 1.0d0/UG%GI%AllCells
    iDir = UConf%UseResume - 1

    write(cStep,*) iStep
    cDirectory = "EXP/" !UConf%SaveDirectiry
    cFileName = trim(adjustl(cDirectory))//"ExactSolution"//trim(adjustl(cStep))//"th_Result.vtk"
    cCaseName = "UnstructuredShockTube" !UConf%CaseName
    open(unit = 1, file =trim(adjustl(cFileName)), status = 'unknown')
        write(1,"('# vtk DataFile Version 3.0')")
        write(1,*) cCaseName
        write(1,"('ASCII')")
        write(1,"('DATASET UNSTRUCTURED_GRID')")
        write(1,"('POINTS ',(1x,i7),' double')") UG%GI%Points

        do iPoint=1, UG%GI%Points
            !write(1,"(3(1x,f22.17))") (UG%CD%Point(iPoint,iLoop),iLoop=1,3) !Adjust Point Number
            write(1,"(3(1x,f22.17))") UG%CD%Point(iPoint,1),UG%CD%Point(iPoint,2),UG%CD%Point(iPoint,3)+1 !for OverSet
        end do
        write(1,*) ""

        write(1,"('CELLS ',2(1x,i7))") UG%GI%RealCells,UG%GI%RealCells*4
        do iCell=1, UG%GI%RealCells
                write(1,"(4(1x,i7))") 3,(UG%Tri%Point(iCell,iLoop)-1,iLoop=1,3) !Adjusted Point Number
        end do
        write(1,*) ""

        write(1,"('CELL_TYPES ',(1x,i7))") UG%GI%RealCells
        do iCell=1, UG%GI%RealCells
            write(1,"(1x,i1)") 5
        end do

        write(1,"('CELL_DATA ',i7)") UG%GI%RealCells
!密度プロット
        write(1,"('SCALARS ExDensity float')")
        write(1,"('LOOKUP_TABLE default')")
        do iCell=1, UG%GI%RealCells
            iPoint = nint(UG%CD%Cell(iCell,iDir)/Width)
            write(1,"((2x,f22.14))") NE%STT%PrimitiveVariable(1,iPoint)
        end do

        write(1,"('VECTORS ExVelocity float')")

            do iCell=1, UG%GI%RealCells
                iPoint = nint(UG%CD%Cell(iCell,iDir)/Width)
                write(1,"((2x,f22.14))") (NE%STT%PrimitiveVariable(iLoop,iPoint),iLoop=2,3),0.0d0
            end do

        write(1,"('SCALARS ExPressure float')")
        write(1,"('LOOKUP_TABLE default')")
        do iCell=1, UG%GI%RealCells
            iPoint = nint(UG%CD%Cell(iCell,iDir)/Width)
            write(1,"((2x,f22.14))") NE%STT%PrimitiveVariable(4,iPoint)
        end do

    close(1)

    return
end subroutine UOutputExact


text:
! vmec_to_gs2_geometry_interface.f90
! Written by Matt Landreman, University of Maryland
! Initial code written August 2017.
! Skip down ~25 lines for detailed description of the input and output parameters.

module vmec_to_gs2_geometry_interface_mod

  implicit none

  private

  public :: vmec_to_gs2_geometry_interface

  real :: theta_pest_target, zeta0
  real, dimension(2) :: vmec_radial_weight_full, vmec_radial_weight_half
  integer, dimension(2) :: vmec_radial_index_full, vmec_radial_index_half

contains

  subroutine vmec_to_gs2_geometry_interface(vmec_filename, nalpha, nzgrid, zeta_center, number_of_field_periods_to_include, &
       desired_normalized_toroidal_flux, vmec_surface_option, verbose, &
       normalized_toroidal_flux_used, safety_factor_q, shat, L_reference, B_reference, &
       alpha, zeta, bmag, gradpar, gds2, gds21, gds22, gbdrift, gbdrift0, cvdrift, cvdrift0)

    use read_wout_mod, nzgrid_vmec => nzgrid  ! VMEC has a variable nzgrid which conflicts with our nzgrid, so rename vmec's version.

    implicit none

    !*********************************************************************
    ! Input parameters
    !*********************************************************************

    ! vmec_filename is the vmec wout_* file that will be read.
    character(*), intent(in) :: vmec_filename

    ! nalpha is the number of grid points in the alpha coordinate:
    integer, intent(in) :: nalpha

    ! The zeta grid has nzgrid*2+1 points, including the "repeated" point at index -nzgrid and +nzgrid.
    integer, intent(in) :: nzgrid

    ! The zeta domain is centered at zeta_center. Setting zeta_center = 2*pi*N/nfp for any integer N should
    ! yield identical results to setting zeta_center = 0, where nfp is the number of field periods (as in VMEC).
    real, intent(in) :: zeta_center

    ! If number_of_field_periods_to_include is > 0, then this parameter does what you think:
    ! the extent of the toroidal in zeta will be 2*pi*number_of_field_periods_to_include/nfp.
    ! If number_of_field_periods_to_include is <= 0, the entire 2*pi toroidal domain will be included.
    real, intent(in) :: number_of_field_periods_to_include

    ! The parameter desired_normalized_toroidal_flux determines which flux surface from the VMEC file will be used
    ! for the computation. This parameter should lie in the interval [0,1].
    real, intent(in) :: desired_normalized_toroidal_flux

    ! If vmec_surface_option = 0, the magnetic surface specified by desired_normalized_toroidal_flux will be used,
    ! by interpolating between the surfaces available in the vmec file.
    ! If vmec_surface_option = 1, the magnetic surface on vmec's HALF radial mesh will be used that is closest to desired_normalized_toroidal_flux.
    ! If vmec_surface_option = 2, the magnetic surface on vmec's FULL radial mesh will be used that is closest to desired_normalized_toroidal_flux.
    ! Other values of vmec_surface_option will cause the program to abort with an error.
    integer, intent(in) :: vmec_surface_option

    ! If verbose is .true., lots of diagnostic information is printed.
    logical, intent(in) :: verbose

    !*********************************************************************
    ! Output quantities
    !*********************************************************************

    ! On exit, normalized_toroidal_flux_used holds the flux surface that was actually used for the geometry,
    ! as measured by psi_toroidal / psi_{toroidal,edge}
    real, intent(out) :: normalized_toroidal_flux_used

    ! Safety factor q = 1/iota
    real, intent(out) :: safety_factor_q

    ! Magnetic shear shat = (x/q) * (d q / d x) where x = Aminor_p * sqrt(psi_toroidal / psi_{toroidal,edge})
    ! and Aminor_p is the minor radius calculated by VMEC.
    real, intent(out) :: shat

    ! L_reference is the reference length used for gs2's normalization, in meters.
    real, intent(out) :: L_reference

    ! B_reference is the reference magnetic field strength used for gs2's normalization, in Tesla.
    real, intent(out) :: B_reference

    ! On exit, alpha holds the grid points in alpha = theta_p - iota * zeta, where theta_p is the PEST toroidal angle
    real, dimension(nalpha), intent(out) :: alpha

    ! On exit, zeta holds the grid points in the toroidal angle zeta
    real, dimension(-nzgrid:nzgrid), intent(out) :: zeta

    real, dimension(nalpha, -nzgrid:nzgrid), intent(out) :: bmag, gradpar, gds2, gds21, gds22, gbdrift, gbdrift0, cvdrift, cvdrift0

    !*********************************************************************
    ! Variables used internally by this subroutine
    !*********************************************************************

    real, parameter :: pi = 3.1415926535897932d+0
    real, parameter :: zero = 0.0d+0
    real, parameter :: one = 1.0d+0
    real, parameter :: mu_0 = 4*pi*(1.0d-7)

    integer :: j, index, izeta, ialpha, which_surface, isurf, m, n, imn, imn_nyq
    real :: angle, sin_angle, cos_angle, temp, edge_toroidal_flux_over_2pi
    real, dimension(:,:), allocatable :: theta_vmec
    integer :: ierr, iopen, fzero_flag
    real :: number_of_field_periods_to_include_final
    real :: dphi, iota, min_dr2, ds, d_pressure_d_s, d_iota_d_s, scale_factor
    real :: theta_vmec_min, theta_vmec_max, sqrt_s
    real, dimension(:), allocatable :: dr2, normalized_toroidal_flux_full_grid, normalized_toroidal_flux_half_grid
    real, dimension(:), allocatable :: d_pressure_d_s_on_half_grid, d_iota_d_s_on_half_grid
    real :: root_solve_absolute_tolerance, root_solve_relative_tolerance
    logical :: non_Nyquist_mode_available, found_imn
    real, dimension(:,:), allocatable :: B, sqrt_g, R, B_dot_grad_theta_pest_over_B_dot_grad_zeta, temp2D
    real, dimension(:,:), allocatable :: d_B_d_theta_vmec, d_B_d_zeta, d_B_d_s
    real, dimension(:,:), allocatable :: d_R_d_theta_vmec, d_R_d_zeta, d_R_d_s
    real, dimension(:,:), allocatable :: d_Z_d_theta_vmec, d_Z_d_zeta, d_Z_d_s
    real, dimension(:,:), allocatable :: d_X_d_theta_vmec, d_X_d_zeta, d_X_d_s
    real, dimension(:,:), allocatable :: d_Y_d_theta_vmec, d_Y_d_zeta, d_Y_d_s
    real, dimension(:,:), allocatable :: d_Lambda_d_theta_vmec, d_Lambda_d_zeta, d_Lambda_d_s
    real, dimension(:,:), allocatable :: B_sub_s, B_sub_theta_vmec, B_sub_zeta
    real, dimension(:,:), allocatable :: B_sup_theta_vmec, B_sup_zeta
    real, dimension(:), allocatable :: d_B_d_s_mnc, d_B_d_s_mns
    real, dimension(:), allocatable :: d_R_d_s_mnc, d_R_d_s_mns
    real, dimension(:), allocatable :: d_Z_d_s_mnc, d_Z_d_s_mns
    real, dimension(:), allocatable :: d_Lambda_d_s_mnc, d_Lambda_d_s_mns
    real, dimension(:,:), allocatable :: grad_s_X, grad_s_Y, grad_s_Z
    real, dimension(:,:), allocatable :: grad_theta_vmec_X, grad_theta_vmec_Y, grad_theta_vmec_Z
    real, dimension(:,:), allocatable :: grad_zeta_X, grad_zeta_Y, grad_zeta_Z
    real, dimension(:,:), allocatable :: grad_psi_X, grad_psi_Y, grad_psi_Z
    real, dimension(:,:), allocatable :: grad_alpha_X, grad_alpha_Y, grad_alpha_Z
    real, dimension(:,:), allocatable :: B_cross_grad_B_dot_grad_alpha, B_cross_grad_B_dot_grad_alpha_alternate
    real, dimension(:,:), allocatable :: B_cross_grad_s_dot_grad_alpha, B_cross_grad_s_dot_grad_alpha_alternate
    real, dimension(:,:), allocatable :: grad_B_X, grad_B_Y, grad_B_Z
    real, dimension(:,:), allocatable :: B_X, B_Y, B_Z

    !*********************************************************************
    ! VMEC variables of interest:
    ! ns = number of flux surfaces used by VMEC
    ! nfp = number of field periods, e.g. 5 for W7-X, 4 for HSX
    ! iotas = rotational transform (1/q) on the half grid.
    ! iotaf = rotational transform on the full grid.
    ! presf = pressure on the full grid.
    !
    ! All VMEC quantities (B, pressure, etc) are in SI units.
    !
    ! In VMEC, quantities on the half grid have the same number of array elements (ns) as quantities on the full grid,
    ! but the first array element is 0.
    !
    !*********************************************************************

    !*********************************************************************
    ! Beginning of executable statements.
    !*********************************************************************

    if (verbose) print *,"Entering subroutine vmec_to_gs2_geometry_interface."

    !*********************************************************************
    ! Do some validation.
    !*********************************************************************

    if (nalpha<1) then
       print *,"Error! nalpha must be >= 1. Instead it is",nalpha
       stop
    end if

    if (nzgrid<1) then
       print *,"Error! nzgrid must be >= 1. Instead it is",nzgrid
       stop
    end if

    if (desired_normalized_toroidal_flux <= 0) then
       print *,"Error! desired_normalized_toroidal_flux must be >0. Instead it is",desired_normalized_toroidal_flux
       stop
    end if

    if (desired_normalized_toroidal_flux > 1) then
       print *,"Error! desired_normalized_toroidal_flux must be <= 1. Instead it is",desired_normalized_toroidal_flux
       stop
    end if

    !*********************************************************************
    ! Read in everything from the vmec wout file using libstell.
    !*********************************************************************

    if (verbose) print *,"  About to read VMEC wout file ",trim(vmec_filename)
    call read_wout_file(vmec_filename, ierr, iopen)
    if (iopen .ne. 0) stop 'error opening wout file'
    if (ierr .ne. 0) stop 'error reading wout file'
    if (verbose) print *,"  Successfully read VMEC data from ",trim(vmec_filename)

    if (verbose) print *,"  Number of field periods (nfp):",nfp
    if (verbose) print *,"  Stellarator-asymmetric? (lasym):",lasym

    ! There is a bug in libstell read_wout_file for ASCII-format wout files, in which the xm_nyq and xn_nyq arrays are sometimes
    ! not populated. The next few lines here provide a workaround:
    if (maxval(abs(xm_nyq)) < 1 .and. maxval(abs(xn_nyq)) < 1) then
       if (mnmax_nyq == mnmax) then
          if (verbose) print *,"xm_nyq and xn_nyq arrays are not populated in the wout file. Using xm and xn instead."
          xm_nyq = xm
          xn_nyq = xn
       else
          print *,"Error! xm_nyq and xn_nyq arrays are not populated in the wout file, and mnmax_nyq != mnmax."
          stop
       end if
    end if

    edge_toroidal_flux_over_2pi = phi(ns) / (2*pi) * isigng ! isigns is called signgs in the wout*.nc file. Why is this signgs here?

    ! Set reference length and magnetic field for GS2's normalization,
    ! using the choices made by Pavlos Xanthopoulos in GIST:
    L_reference = Aminor ! Note that 'Aminor' in read_wout_mod is called 'Aminor_p' in the wout*.nc file.
    B_reference = 2 * abs(edge_toroidal_flux_over_2pi) / (L_reference * L_reference)
    if (verbose) then
       print *,"  Reference length for GS2 normalization:",L_reference," meters."
       print *,"  Reference magnetic field strength for GS2 normalization:",B_reference," Tesla."
    end if

    ! --------------------------------------------------------------------------------
    ! Do some sanity checking to ensure the VMEC arrays have some expected properties.
    ! --------------------------------------------------------------------------------

    ! 'phi' is vmec's array of the toroidal flux (not divided by 2pi!) on vmec's radial grid.
    if (abs(phi(1)) > 1d-14) then
       print *,"Error! VMEC phi array does not begin with 0."
       print *,"phi:",phi
       stop
    end if

    dphi = phi(2) - phi(1)
    do j=3,ns
       if (abs(phi(j)-phi(j-1)-dphi) > 1d-11) then
          print *,"Error! VMEC phi array is not uniformly spaced."
          print *,"phi:",phi
          stop
       end if
    end do

    ! The variable called 'phips' in the wout file is called just 'phip' in read_wout_mod.F.
    ! phips is on the half-mesh, so skip first point.
    do j=2,ns
       if (abs(phip(j)+phi(ns)/(2*pi)) > 1d-11) then
          print *,"Error! VMEC phips array is not constant and equal to -phi(ns)/(2*pi)."
          print *,"phip(s):",phip
          stop
       end if
    end do

    ! The first mode in the m and n arrays should be m=n=0:
    if (xm(1) .ne. 0) stop "First element of xm in the wout file should be 0."
    if (xn(1) .ne. 0) stop "First element of xn in the wout file should be 0."
    if (xm_nyq(1) .ne. 0) stop "First element of xm_nyq in the wout file should be 0."
    if (xn_nyq(1) .ne. 0) stop "First element of xn_nyq in the wout file should be 0."

    ! Lambda should be on the half mesh, so its value at radial index 1 should be 0 for all (m,n)
    if (maxval(abs(lmns(:,1))) > 0) then
       print *,"Error! Expected lmns to be on the half mesh, but its value at radial index 1 is nonzero."
       print *,"Here comes lmns(:,1):", lmns(:,1)
       stop
    end if
    if (lasym) then
       if (maxval(abs(lmnc(:,1))) > 0) then
          print *,"Error! Expected lmnc to be on the half mesh, but its value at radial index 1 is nonzero."
          print *,"Here comes lmnc(:,1):", lmnc(:,1)
          stop
       end if
    end if

    ! --------------------------------------------------------------------------------
    ! End of sanity checks.
    ! --------------------------------------------------------------------------------

    allocate(normalized_toroidal_flux_full_grid(ns))
    normalized_toroidal_flux_full_grid = [( real(j-1)/(ns-1), j=1,ns )]

    ! Build an array of the half grid points:
    allocate(normalized_toroidal_flux_half_grid(ns-1))
    do j = 1,ns-1
       normalized_toroidal_flux_half_grid(j) = (normalized_toroidal_flux_full_grid(j) + normalized_toroidal_flux_full_grid(j+1))*(0.5d+0)
    end do

    !*********************************************************************
    ! Determine which flux surface to use, based on
    ! desired_normalized_toroidal_flux and vmec_surface_option.
    !*********************************************************************

    ! Possible values of vmec_surface_option:
    ! 0 = Use the exact radius requested.
    ! 1 = Use the nearest value of the VMEC half grid.
    ! 2 = Use the nearest value of the VMEC full grid.

    select case (vmec_surface_option)
    case (0)
       ! Use exact radius requested.
       normalized_toroidal_flux_used = desired_normalized_toroidal_flux

    case (1)
       ! Use nearest value of the VMEC half grid

       ! Compute differences
       allocate(dr2(ns-1))
       dr2 = (normalized_toroidal_flux_half_grid - desired_normalized_toroidal_flux) ** 2

       index = 1
       min_dr2 = dr2(1)
       ! Find the index of minimum error:
       do j=2,ns-1
          if (dr2(j)<min_dr2) then
             index = j
             min_dr2 = dr2(j)
          end if
       end do

       normalized_toroidal_flux_used = normalized_toroidal_flux_half_grid(index)
       deallocate(dr2)

    case (2)
       ! Use nearest value of the VMEC full grid

       ! Compute differences
       allocate(dr2(ns))
       dr2 = (normalized_toroidal_flux_full_grid - desired_normalized_toroidal_flux) ** 2

       index = 1
       min_dr2 = dr2(1)
       ! Find the index of minimum error:
       do j=2,ns
          if (dr2(j)<min_dr2) then
             index = j
             min_dr2 = dr2(j)
          end if
       end do

       normalized_toroidal_flux_used = normalized_toroidal_flux_full_grid(index)
       deallocate(dr2)

    case default
       print *,"Error! vmec_surface_option must be 0, 1, or 2. It is instead ",vmec_surface_option
       stop
    end select

    ! --------------------------------------------------------------------------------
    ! Done choosing the actual radius to use.
    ! --------------------------------------------------------------------------------

    ! In general, we get quantities for gs2 by linear interpolation, taking a weighted average of the quantity from
    ! 2 surfaces in the VMEC file. Sometimes the weights are 0 and 1, i.e. no interpolation is needed.

    ! For any VMEC quantity Q on the full grid, the value used in GS2 will be
    !  Q_gs2 = Q(vmec_radial_index_full(1))*vmec_radial_weight_full(1) + Q(vmec_radial_index_full(2))*vmec_radial_weight_full(2)

    ! For any VMEC quantity Q on the half grid, the value used in GS2 will be
    !  Q_gs2 = Q(vmec_radial_index_half(1))*vmec_radial_weight_half(1) + Q(vmec_radial_index_half(2))*vmec_radial_weight_half(2)

    ! Handle quantities for the full grid
    if (normalized_toroidal_flux_used>1) then
       stop "Error! normalized_toroidal_flux_used cannot be >1"
    elseif (normalized_toroidal_flux_used<0) then
       stop "Error! normalized_toroidal_flux_used cannot be <0"
    elseif (normalized_toroidal_flux_used==1) then
       vmec_radial_index_full(1) = ns-1
       vmec_radial_index_full(2) = ns
       vmec_radial_weight_full(1) = zero
    else
       ! normalized_toroidal_flux_used is >= 0 and <1
       ! This is the most common case.
       vmec_radial_index_full(1) = floor(normalized_toroidal_flux_used*(ns-1))+1
       vmec_radial_index_full(2) = vmec_radial_index_full(1) + 1
       vmec_radial_weight_full(1) = vmec_radial_index_full(1) - normalized_toroidal_flux_used*(ns-one)
    end if
    vmec_radial_weight_full(2) = one - vmec_radial_weight_full(1)

    ! Handle quantities for the half grid
    if (normalized_toroidal_flux_used < normalized_toroidal_flux_half_grid(1)) then
       print *,"Warning: extrapolating beyond the end of VMEC's half grid."
       print *,"(Extrapolating towards the magnetic axis.) Results are likely to be inaccurate."

       ! We start at element 2 since element 1 is always 0 for quantities on the half grid.
       vmec_radial_index_half(1) = 2
       vmec_radial_index_half(2) = 3
       vmec_radial_weight_half(1) = (normalized_toroidal_flux_half_grid(2) - normalized_toroidal_flux_used) / (normalized_toroidal_flux_half_grid(2) - normalized_toroidal_flux_half_grid(1))

    elseif (normalized_toroidal_flux_used > normalized_toroidal_flux_half_grid(ns-1)) then
       print *,"Warning: extrapolating beyond the end of VMEC's half grid."
       print *,"(Extrapolating towards the last closed flux surface.) Results may be inaccurate."
       vmec_radial_index_half(1) = ns-1
       vmec_radial_index_half(2) = ns
       vmec_radial_weight_half(1) = (normalized_toroidal_flux_half_grid(ns-1) - normalized_toroidal_flux_used) &
            / (normalized_toroidal_flux_half_grid(ns-1) - normalized_toroidal_flux_half_grid(ns-2))

    elseif (normalized_toroidal_flux_used == normalized_toroidal_flux_half_grid(ns-1)) then
       ! We are exactly at the last point of the half grid
       vmec_radial_index_half(1) = ns-1
       vmec_radial_index_half(2) = ns
       vmec_radial_weight_half(1) = zero
    else
       ! normalized_toroidal_flux_used is inside the half grid.
       ! This is the most common case.
       vmec_radial_index_half(1) = floor(normalized_toroidal_flux_used*(ns-1) + 0.5d+0)+1
       if (vmec_radial_index_half(1) < 2) then
          ! This can occur sometimes due to roundoff error.
          vmec_radial_index_half(1) = 2
       end if
       vmec_radial_index_half(2) = vmec_radial_index_half(1) + 1
       vmec_radial_weight_half(1) = vmec_radial_index_half(1) - normalized_toroidal_flux_used*(ns-one) - (0.5d+0)
    end if
    vmec_radial_weight_half(2) = one-vmec_radial_weight_half(1)

    if (verbose) then
       if (abs(vmec_radial_weight_half(1)) < 1e-14) then
          print "(a,i4,a,i4,a)","   Using radial index ",vmec_radial_index_half(2)," of ",ns," from vmec's half mesh."
       elseif (abs(vmec_radial_weight_half(2)) < 1e-14) then
          print "(a,i4,a,i4,a)","   Using radial index ",vmec_radial_index_half(1)," of ",ns," from vmec's half mesh."
       else
          print "(a,i4,a,i4,a,i4,a)", "   Interpolating using radial indices ",vmec_radial_index_half(1)," and ",vmec_radial_index_half(2),&
               " of ",ns," from vmec's half mesh."
          print "(a,f17.14,a,f17.14)", "   Weights for half mesh = ",vmec_radial_weight_half(1)," and ",vmec_radial_weight_half(2)
          print "(a,i4,a,i4,a,i4,a)", "   Interpolating using radial indices ",vmec_radial_index_full(1)," and ",vmec_radial_index_full(2),&
               " of ",ns," from vmec's full mesh."
          print "(a,f17.14,a,f17.14)", "   Weights for full mesh = ",vmec_radial_weight_full(1)," and ",vmec_radial_weight_full(2)
       end if
    end if

    !*********************************************************************
    ! Evaluate several radial-profile functions at the flux surface
    ! we ended up choosing.
    !*********************************************************************

    iota = iotas(vmec_radial_index_half(1)) * vmec_radial_weight_half(1) &
         + iotas(vmec_radial_index_half(2)) * vmec_radial_weight_half(2)
    if (verbose) print *,"  iota =",iota
    safety_factor_q = 1/iota

    allocate(d_iota_d_s_on_half_grid(ns))
    d_iota_d_s_on_half_grid = 0
    ds = normalized_toroidal_flux_full_grid(2) - normalized_toroidal_flux_full_grid(1)
    if (verbose) print *,"  ds =",ds
    d_iota_d_s_on_half_grid(2:ns) = (iotaf(2:ns) - iotaf(1:ns-1)) / ds
    d_iota_d_s =  &
         d_iota_d_s_on_half_grid(vmec_radial_index_half(1)) * vmec_radial_weight_half(1) &
         + d_iota_d_s_on_half_grid(vmec_radial_index_half(2)) * vmec_radial_weight_half(2)
    deallocate(d_iota_d_s_on_half_grid)
    if (verbose) print *,"  d iota / d s =",d_iota_d_s
    ! shat = (r/q)(dq/dr) where r = a sqrt(s).
    !      = - (r/iota) (d iota / d r) = -2 (s/iota) (d iota / d s)
    shat = (-2 * normalized_toroidal_flux_used / iota) * d_iota_d_s

    allocate(d_pressure_d_s_on_half_grid(ns))
    d_pressure_d_s_on_half_grid = 0
    ds = normalized_toroidal_flux_full_grid(2) - normalized_toroidal_flux_full_grid(1)
    d_pressure_d_s_on_half_grid(2:ns) = (presf(2:ns) - presf(1:ns-1)) / ds
    d_pressure_d_s =  &
         d_pressure_d_s_on_half_grid(vmec_radial_index_half(1)) * vmec_radial_weight_half(1) &
         + d_pressure_d_s_on_half_grid(vmec_radial_index_half(2)) * vmec_radial_weight_half(2)
    deallocate(d_pressure_d_s_on_half_grid)
    if (verbose) print *,"  d pressure / d s =",d_pressure_d_s

    !*********************************************************************
    ! Set up the coordinate grids.
    !*********************************************************************

    alpha = [( ((j-1)*2*pi) / nalpha, j=1, nalpha )]

!!$    if (number_of_field_periods_to_include > nfp) then
!!$       print *,"Error! number_of_field_periods_to_include > nfp"
!!$       print *,"  number_of_field_periods_to_include =",number_of_field_periods_to_include
!!$       print *,"  nfp =",nfp
!!$       stop
!!$    end if
    number_of_field_periods_to_include_final = number_of_field_periods_to_include
    if (number_of_field_periods_to_include <= 0) then
       number_of_field_periods_to_include_final = nfp
       if (verbose) print *,"  Since number_of_field_periods_to_include was <= 0, it is being reset to nfp =",nfp
    end if

    zeta = [( zeta_center + (pi*j*number_of_field_periods_to_include_final)/(nfp*nzgrid), j=-nzgrid, nzgrid )]

    !*********************************************************************
    ! We know theta_pest = alpha + iota * zeta, but we need to determine
    ! theta_vmec = theta_pest - Lambda.
    !*********************************************************************

    allocate(theta_vmec(nalpha, -nzgrid:nzgrid))

    if (verbose) print *,"  Beginning root solves to determine theta_vmec."
    root_solve_absolute_tolerance = 1.0d-10
    root_solve_relative_tolerance = 1.0d-10
    do izeta = -nzgrid, nzgrid
       zeta0 = zeta(izeta)
       do ialpha = 1,nalpha
          theta_pest_target = alpha(ialpha) + iota * zeta0
          ! Guess that theta_vmec will be within 0.3 radians of theta_pest:
          theta_vmec_min = theta_pest_target - 1.0
          theta_vmec_max = theta_pest_target + 1.0

          ! In the 4th argument, we are telling the root-finder (fzero) to use theta_pest as the initial guess for theta_vmec.
          call fzero(fzero_residual, theta_vmec_min, theta_vmec_max, theta_pest_target, &
               root_solve_relative_tolerance, root_solve_absolute_tolerance, fzero_flag)
          ! Note: fzero returns its answer in theta_vmec_min.
          theta_vmec(ialpha,izeta) = theta_vmec_min
          if (fzero_flag == 4) then
             stop "ERROR: fzero returned error 4: no sign change in residual"
          else if (fzero_flag > 2) then
             print *,"WARNING: fzero returned an error code:",fzero_flag
          end if
       end do
    end do
    if (verbose) then
       print *,"  Done with root solves. Here comes theta_vmec:"
       do j = 1, nalpha
          !print *,theta_vmec(j,:)
       end do
    end if

    !*********************************************************************
    ! Initialize geometry arrays
    !*********************************************************************

    bmag = 0
    gradpar = 0
    gds2 = 0
    gds21 = 0
    gds22 = 0
    gbdrift = 0
    gbdrift0 = 0
    cvdrift = 0
    cvdrift0 = 0

    allocate(B(nalpha,-nzgrid:nzgrid))
    allocate(temp2D(nalpha,-nzgrid:nzgrid))
    allocate(sqrt_g(nalpha,-nzgrid:nzgrid))
    allocate(R(nalpha,-nzgrid:nzgrid))
    allocate(d_B_d_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(d_B_d_zeta(nalpha,-nzgrid:nzgrid))
    allocate(d_B_d_s(nalpha,-nzgrid:nzgrid))
    allocate(d_R_d_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(d_R_d_zeta(nalpha,-nzgrid:nzgrid))
    allocate(d_R_d_s(nalpha,-nzgrid:nzgrid))
    allocate(d_Z_d_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(d_Z_d_zeta(nalpha,-nzgrid:nzgrid))
    allocate(d_Z_d_s(nalpha,-nzgrid:nzgrid))
    allocate(d_Lambda_d_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(d_Lambda_d_zeta(nalpha,-nzgrid:nzgrid))
    allocate(d_Lambda_d_s(nalpha,-nzgrid:nzgrid))
    allocate(B_sub_s(nalpha,-nzgrid:nzgrid))
    allocate(B_sub_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(B_sub_zeta(nalpha,-nzgrid:nzgrid))
    allocate(B_sup_theta_vmec(nalpha,-nzgrid:nzgrid))
    allocate(B_sup_zeta(nalpha,-nzgrid:nzgrid))

    allocate(d_B_d_s_mnc(ns))
    allocate(d_B_d_s_mns(ns))
    allocate(d_R_d_s_mnc(ns))
    allocate(d_R_d_s_mns(ns))
    allocate(d_Z_d_s_mnc(ns))
    allocate(d_Z_d_s_mns(ns))
    allocate(d_Lambda_d_s_mnc(ns))
    allocate(d_Lambda_d_s_mns(ns))

    allocate(d_X_d_s(nalpha, -nzgrid:nzgrid))
    allocate(d_X_d_theta_vmec(nalpha, -nzgrid:nzgrid))
    allocate(d_X_d_zeta(nalpha, -nzgrid:nzgrid))
    allocate(d_Y_d_s(nalpha, -nzgrid:nzgrid))
    allocate(d_Y_d_theta_vmec(nalpha, -nzgrid:nzgrid))
    allocate(d_Y_d_zeta(nalpha, -nzgrid:nzgrid))

    allocate(grad_s_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_s_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_s_Z(nalpha, -nzgrid:nzgrid))
    allocate(grad_theta_vmec_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_theta_vmec_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_theta_vmec_Z(nalpha, -nzgrid:nzgrid))
    allocate(grad_zeta_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_zeta_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_zeta_Z(nalpha, -nzgrid:nzgrid))
    allocate(grad_psi_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_psi_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_psi_Z(nalpha, -nzgrid:nzgrid))
    allocate(grad_alpha_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_alpha_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_alpha_Z(nalpha, -nzgrid:nzgrid))

    allocate(B_X(nalpha, -nzgrid:nzgrid))
    allocate(B_Y(nalpha, -nzgrid:nzgrid))
    allocate(B_Z(nalpha, -nzgrid:nzgrid))
    allocate(grad_B_X(nalpha, -nzgrid:nzgrid))
    allocate(grad_B_Y(nalpha, -nzgrid:nzgrid))
    allocate(grad_B_Z(nalpha, -nzgrid:nzgrid))
    allocate(B_cross_grad_B_dot_grad_alpha(nalpha, -nzgrid:nzgrid))
    allocate(B_cross_grad_B_dot_grad_alpha_alternate(nalpha, -nzgrid:nzgrid))
    allocate(B_cross_grad_s_dot_grad_alpha(nalpha, -nzgrid:nzgrid))
    allocate(B_cross_grad_s_dot_grad_alpha_alternate(nalpha, -nzgrid:nzgrid))

    B = 0
    sqrt_g = 0
    R = 0
    d_B_d_theta_vmec = 0
    d_B_d_zeta = 0
    d_B_d_s = 0
    d_R_d_theta_vmec = 0
    d_R_d_zeta = 0
    d_R_d_s = 0
    d_Z_d_theta_vmec = 0
    d_Z_d_zeta = 0
    d_Z_d_s = 0
    d_Lambda_d_theta_vmec = 0
    d_Lambda_d_zeta = 0
    d_Lambda_d_s = 0
    B_sub_s = 0
    B_sub_theta_vmec = 0
    B_sub_zeta = 0
    B_sup_theta_vmec = 0
    B_sup_zeta = 0

    !*********************************************************************
    ! Now that we know the grid points in theta_vmec, we can evaluate
    ! all the geometric quantities on the grid points.
    !*********************************************************************

    do imn_nyq = 1, mnmax_nyq ! All the quantities we need except R, Z, and Lambda use the _nyq mode numbers.
       m = xm_nyq(imn_nyq)
       n = xn_nyq(imn_nyq)/nfp

       if (abs(m) >= mpol .or. abs(n) > ntor) then
          non_Nyquist_mode_available = .false.
       else
          non_Nyquist_mode_available = .true.
          ! Find the imn in the non-Nyquist arrays that corresponds to the same m and n.
          found_imn = .false.
          do imn = 1,mnmax
             if (xm(imn)==m .and. xn(imn)==n*nfp) then
                found_imn = .true.
                exit
             end if
          end do
          if ((xm(imn) .ne. m) .or. (xn(imn) .ne. n*nfp)) stop "Something went wrong!"
          if (.not. found_imn) stop "Error! imn could not be found matching the given imn_nyq."
       end if

       ! All quantities are multiplied by a variable scale_factor which can in principle depend on m and n.
       ! For now we just set scale_factor = 1. In the future, scale_factor could be used to lower the
       ! symmetry-breaking Fourier components, or filter out certain Fourier components in some way.
       scale_factor = 1

       ! -----------------------------------------------------
       ! First, consider just the stellarator-symmetric terms:
       ! -----------------------------------------------------

       ! Evaluate the radial derivatives we will need:

       ! B and Lambda are on the half mesh, so their radial derivatives are on the full mesh.
       ! R and Z are on the full mesh, so their radial derivatives are on the half mesh.

       d_B_d_s_mnc(2:ns-1) = (bmnc(imn_nyq,3:ns) - bmnc(imn_nyq,2:ns-1)) / ds
       ! Simplistic extrapolation at the endpoints:
       d_B_d_s_mnc(1) = d_B_d_s_mnc(2)
       d_B_d_s_mnc(ns) = d_B_d_s_mnc(ns-1)

       if (non_Nyquist_mode_available) then
          ! R is on the full mesh:
          d_R_d_s_mnc(2:ns) = (rmnc(imn,2:ns) - rmnc(imn,1:ns-1)) / ds
          d_R_d_s_mnc(1) = 0

          ! Z is on the full mesh:
          d_Z_d_s_mns(2:ns) = (zmns(imn,2:ns) - zmns(imn,1:ns-1)) / ds
          d_Z_d_s_mns(1) = 0

          ! Lambda is on the half mesh:
          d_Lambda_d_s_mns(2:ns-1) = (lmns(imn,3:ns) - lmns(imn,2:ns-1)) / ds
          ! Simplistic extrapolation at the endpoints:
          d_Lambda_d_s_mns(1) = d_Lambda_d_s_mns(2)
          d_Lambda_d_s_mns(ns) = d_Lambda_d_s_mns(ns-1)
       else
          d_R_d_s_mnc = 0
          d_Z_d_s_mns = 0
          d_Lambda_d_s_mns = 0
       end if

       ! End of evaluating radial derivatives.

       do ialpha = 1,nalpha
          do izeta = -nzgrid, nzgrid
             angle = m * theta_vmec(ialpha,izeta) - n * nfp * zeta(izeta)
             cos_angle = cos(angle)
             sin_angle = sin(angle)

             do isurf = 1,2

                ! Handle |B|:
                temp = bmnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                B(ialpha,izeta) = B(ialpha,izeta) + temp * cos_angle
                d_B_d_theta_vmec(ialpha,izeta) = d_B_d_theta_vmec(ialpha,izeta) - m * temp * sin_angle
                d_B_d_zeta(ialpha,izeta) = d_B_d_zeta(ialpha,izeta) + n * nfp * temp * sin_angle

                ! Handle Jacobian:
                temp = gmnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                sqrt_g(ialpha,izeta) = sqrt_g(ialpha,izeta) + temp * cos_angle

                ! Handle B sup theta:
                temp = bsupumnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                B_sup_theta_vmec(ialpha,izeta) = B_sup_theta_vmec(ialpha,izeta) + temp * cos_angle

                ! Handle B sup zeta:
                temp = bsupvmnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                B_sup_zeta(ialpha,izeta) = B_sup_zeta(ialpha,izeta) + temp * cos_angle

                ! Handle B sub theta:
                temp = bsubumnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                B_sub_theta_vmec(ialpha,izeta) = B_sub_theta_vmec(ialpha,izeta) + temp * cos_angle

                ! Handle B sub zeta:
                temp = bsubvmnc(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                temp = temp*scale_factor
                B_sub_zeta(ialpha,izeta) = B_sub_zeta(ialpha,izeta) + temp * cos_angle

                ! Handle B sub psi.
                ! Unlike the other components of B, this one is on the full mesh.
                temp = bsubsmns(imn_nyq,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                temp = temp*scale_factor
                B_sub_s(ialpha,izeta) = B_sub_s(ialpha,izeta) + temp * sin_angle

                ! Handle d B / d s
                ! Since bmnc is on the half mesh, its radial derivative is on the full mesh.
                temp = d_B_d_s_mnc(vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                temp = temp*scale_factor
                d_B_d_s(ialpha,izeta) = d_B_d_s(ialpha,izeta) + temp * cos_angle

                ! Handle arrays that use xm and xn instead of xm_nyq and xn_nyq.
                if (non_Nyquist_mode_available) then

                   ! Handle R, which is on the full mesh
                   temp = rmnc(imn,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                   temp = temp*scale_factor
                   R(ialpha,izeta) = R(ialpha,izeta) + temp * cos_angle
                   d_R_d_theta_vmec(ialpha,izeta) = d_R_d_theta_vmec(ialpha,izeta) - temp * m * sin_angle
                   d_R_d_zeta(ialpha,izeta)  = d_R_d_zeta(ialpha,izeta)  + temp * n * nfp * sin_angle

                   ! Handle Z, which is on the full mesh
                   temp = zmns(imn,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                   temp = temp*scale_factor
                   !Z(ialpha,izeta) = Z(ialpha,izeta) + temp * sin_angle  ! We don't actually need Z itself, only derivatives of Z.
                   d_Z_d_theta_vmec(ialpha,izeta) = d_Z_d_theta_vmec(ialpha,izeta) + temp * m * cos_angle
                   d_Z_d_zeta(ialpha,izeta)  = d_Z_d_zeta(ialpha,izeta)  - temp * n * nfp * cos_angle

                   ! Handle Lambda:
                   temp = lmns(imn,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   ! We don't need Lambda itself, just its derivatives.
                   d_Lambda_d_theta_vmec(ialpha,izeta) = d_Lambda_d_theta_vmec(ialpha,izeta) + m * temp * cos_angle
                   d_Lambda_d_zeta(ialpha,izeta) = d_Lambda_d_zeta(ialpha,izeta) - n * nfp * temp * cos_angle

                   ! Handle d R / d s
                   ! Since R is on the full mesh, its radial derivative is on the half mesh.
                   temp = d_R_d_s_mnc(vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   d_R_d_s(ialpha,izeta) = d_R_d_s(ialpha,izeta) + temp * cos_angle

                   ! Handle d Z / d s
                   ! Since Z is on the full mesh, its radial derivative is on the half mesh.
                   temp = d_Z_d_s_mns(vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   d_Z_d_s(ialpha,izeta) = d_Z_d_s(ialpha,izeta) + temp * sin_angle

                   ! Handle d Lambda / d s
                   ! Since Lambda is on the half mesh, its radial derivative is on the full mesh.
                   temp = d_Lambda_d_s_mns(vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                   temp = temp*scale_factor
                   d_Lambda_d_s(ialpha,izeta) = d_Lambda_d_s(ialpha,izeta) + temp * sin_angle

                end if
             end do
          end do
       end do

       ! -----------------------------------------------------
       ! Now consider the stellarator-asymmetric terms.
       ! -----------------------------------------------------

       if (lasym) then

          ! Evaluate the radial derivatives we will need:

          ! B and Lambda are on the half mesh, so their radial derivatives are on the full mesh.
          ! R and Z are on the full mesh, so their radial derivatives are on the half mesh.

          d_B_d_s_mns(2:ns-1) = (bmns(imn_nyq,3:ns) - bmns(imn_nyq,2:ns-1)) / ds
          ! Simplistic extrapolation at the endpoints:
          d_B_d_s_mns(1) = d_B_d_s_mns(2)
          d_B_d_s_mns(ns) = d_B_d_s_mns(ns-1)

          if (non_Nyquist_mode_available) then
             ! R is on the full mesh:
             d_R_d_s_mns(2:ns) = (rmns(imn,2:ns) - rmns(imn,1:ns-1)) / ds
             d_R_d_s_mns(1) = 0

             ! Z is on the full mesh:
             d_Z_d_s_mnc(2:ns) = (zmnc(imn,2:ns) - zmnc(imn,1:ns-1)) / ds
             d_Z_d_s_mnc(1) = 0

             ! Lambda is on the half mesh:
             d_Lambda_d_s_mnc(2:ns-1) = (lmnc(imn_nyq,3:ns) - lmnc(imn_nyq,2:ns-1)) / ds
             ! Simplistic extrapolation at the endpoints:
             d_Lambda_d_s_mnc(1) = d_Lambda_d_s_mnc(2)
             d_Lambda_d_s_mnc(ns) = d_Lambda_d_s_mnc(ns-1)
          else
             d_R_d_s_mns = 0
             d_Z_d_s_mnc = 0
             d_Lambda_d_s_mnc = 0
          end if

          ! End of evaluating radial derivatives.

          do ialpha = 1,nalpha
             do izeta = -nzgrid, nzgrid
                angle = m * theta_vmec(ialpha,izeta) - n * nfp * zeta(izeta)
                cos_angle = cos(angle)
                sin_angle = sin(angle)

                do isurf = 1,2

                   ! Handle |B|:
                   temp = bmns(imn_nyq,vmec_radial_index_half(1)) * vmec_radial_weight_half(1)
                   temp = temp * scale_factor
                   B(ialpha,izeta) = B(ialpha,izeta) + temp * sin_angle
                   d_B_d_theta_vmec(ialpha,izeta) = d_B_d_theta_vmec(ialpha,izeta) + m * temp * cos_angle
                   d_B_d_zeta(ialpha,izeta) = d_B_d_zeta(ialpha,izeta) - n * nfp * temp * cos_angle

                   ! Handle Jacobian:
                   temp = gmns(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   sqrt_g(ialpha,izeta) = sqrt_g(ialpha,izeta) + temp * sin_angle

                   ! Handle B sup theta:
                   temp = bsupumns(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   B_sup_theta_vmec(ialpha,izeta) = B_sup_theta_vmec(ialpha,izeta) + temp * sin_angle

                   ! Handle B sup zeta:
                   temp = bsupvmns(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   B_sup_zeta(ialpha,izeta) = B_sup_zeta(ialpha,izeta) + temp * sin_angle

                   ! Handle B sub theta:
                   temp = bsubumns(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   B_sub_theta_vmec(ialpha,izeta) = B_sub_theta_vmec(ialpha,izeta) + temp * sin_angle

                   ! Handle B sub zeta:
                   temp = bsubvmns(imn_nyq,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                   temp = temp*scale_factor
                   B_sub_zeta(ialpha,izeta) = B_sub_zeta(ialpha,izeta) + temp * sin_angle

                   ! Handle B sub psi.
                   ! Unlike the other components of B, this one is on the full mesh.
                   temp = bsubsmnc(imn_nyq,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                   temp = temp*scale_factor
                   B_sub_s(ialpha,izeta) = B_sub_s(ialpha,izeta) + temp * cos_angle

                   ! Handle d B / d s.
                   ! Since bmns is on the half mesh, its radial derivative is on the full mesh.
                   temp = d_B_d_s_mns(vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                   temp = temp*scale_factor
                   d_B_d_s(ialpha,izeta) = d_B_d_s(ialpha,izeta) + temp * sin_angle

                   ! Handle arrays that use xm and xn instead of xm_nyq and xn_nyq.
                   if (non_Nyquist_mode_available) then

                      ! Handle R, which is on the full mesh
                      temp = rmns(imn,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                      temp = temp*scale_factor
                      R(ialpha,izeta) = R(ialpha,izeta) + temp * sin_angle
                      d_R_d_theta_vmec(ialpha,izeta) = d_R_d_theta_vmec(ialpha,izeta) + temp * m * cos_angle
                      d_R_d_zeta(ialpha,izeta)  = d_R_d_zeta(ialpha,izeta)  - temp * n * nfp * cos_angle

                      ! Handle Z, which is on the full mesh
                      temp = zmnc(imn,vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                      temp = temp*scale_factor
                      ! Z(ialpha,izeta) = Z(ialpha,izeta) + temp * cos_angle   ! We don't actually need Z itself, only derivatives of Z.
                      d_Z_d_theta_vmec(ialpha,izeta) = d_Z_d_theta_vmec(ialpha,izeta) - temp * m * sin_angle
                      d_Z_d_zeta(ialpha,izeta)  = d_Z_d_zeta(ialpha,izeta)  + temp * n * nfp * sin_angle

                      ! Handle Lambda, which is on the half mesh
                      temp = lmnc(imn,vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                      temp = temp*scale_factor
                      ! We don't actually need Lambda itself, only derivatives of Lambda.
                      d_Lambda_d_theta_vmec(ialpha,izeta) = d_Lambda_d_theta_vmec(ialpha,izeta) - temp * m * sin_angle
                      d_Lambda_d_zeta(ialpha,izeta)  = d_Lambda_d_zeta(ialpha,izeta)  + temp * n * nfp * sin_angle

                      ! Handle d R / d s.
                      ! Since R is on the full mesh, its radial derivative is on the half mesh.
                      temp = d_R_d_s_mns(vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                      temp = temp*scale_factor
                      d_R_d_s(ialpha,izeta) = d_R_d_s(ialpha,izeta) + temp * sin_angle

                      ! Handle d Z / d s.
                      ! Since Z is on the full mesh, its radial derivative is on the half mesh.
                      temp = d_Z_d_s_mnc(vmec_radial_index_half(isurf)) * vmec_radial_weight_half(isurf)
                      temp = temp*scale_factor
                      d_Z_d_s(ialpha,izeta) = d_Z_d_s(ialpha,izeta) + temp * cos_angle

                      ! Handle d Lambda / d s.
                      ! Since Lambda is on the half mesh, its radial derivative is on the full mesh.
                      temp = d_Lambda_d_s_mnc(vmec_radial_index_full(isurf)) * vmec_radial_weight_full(isurf)
                      temp = temp*scale_factor
                      d_Lambda_d_s(ialpha,izeta) = d_Lambda_d_s(ialpha,izeta) + temp * cos_angle
                   end if
                end do
             end do
          end do
       end if
    end do

    !*********************************************************************
    ! Sanity check: If the conversion to theta_pest has been done
    ! correctly, we should find that
    ! (B dot grad theta_pest) / (B dot grad zeta) = iota.
    ! Let's verify this:
    !*********************************************************************

    allocate(B_dot_grad_theta_pest_over_B_dot_grad_zeta(nalpha, -nzgrid:nzgrid))
    ! Compute (B dot grad theta_pest) / (B dot grad zeta):
    B_dot_grad_theta_pest_over_B_dot_grad_zeta = (B_sup_theta_vmec * (1 + d_Lambda_d_theta_vmec) + B_sup_zeta * d_Lambda_d_zeta) / B_sup_zeta
    temp2D = iota
    call test_arrays(B_dot_grad_theta_pest_over_B_dot_grad_zeta, temp2D, .false., 0.035, 'iota')
    deallocate(B_dot_grad_theta_pest_over_B_dot_grad_zeta)

    !*********************************************************************
    ! Using R(theta,zeta) and Z(theta,zeta), compute the Cartesian
    ! components of the gradient basis vectors using the dual relations:
    !*********************************************************************

    do izeta = -nzgrid, nzgrid
       cos_angle = cos(zeta(izeta))
       sin_angle = sin(zeta(izeta))

       ! X = R * cos(zeta)
       d_X_d_theta_vmec(:,izeta) = d_R_d_theta_vmec(:,izeta) * cos_angle
       d_X_d_zeta(:,izeta) = d_R_d_zeta(:,izeta) * cos_angle - R(:,izeta) * sin_angle
       d_X_d_s(:,izeta) = d_R_d_s(:,izeta) * cos_angle

       ! Y = R * sin(zeta)
       d_Y_d_theta_vmec(:,izeta) = d_R_d_theta_vmec(:,izeta) * sin_angle
       d_Y_d_zeta(:,izeta) = d_R_d_zeta(:,izeta) * sin_angle + R(:,izeta) * cos_angle
       d_Y_d_s(:,izeta) = d_R_d_s(:,izeta) * sin_angle

!!$       ! Y = -R * sin(zeta)
!!$       d_Y_d_theta_vmec(:,izeta) = -d_R_d_theta_vmec(:,izeta) * sin_angle
!!$       d_Y_d_zeta(:,izeta) = -(d_R_d_zeta(:,izeta) * sin_angle + R(:,izeta) * cos_angle)
!!$       d_Y_d_s(:,izeta) = -d_R_d_s(:,izeta) * sin_angle

    end do

    ! Use the dual relations to get the Cartesian components of grad s, grad theta_vmec, and grad zeta:
    grad_s_X = (d_Y_d_theta_vmec * d_Z_d_zeta - d_Z_d_theta_vmec * d_Y_d_zeta) / sqrt_g
    grad_s_Y = (d_Z_d_theta_vmec * d_X_d_zeta - d_X_d_theta_vmec * d_Z_d_zeta) / sqrt_g
    grad_s_Z = (d_X_d_theta_vmec * d_Y_d_zeta - d_Y_d_theta_vmec * d_X_d_zeta) / sqrt_g

    grad_theta_vmec_X = (d_Y_d_zeta * d_Z_d_s - d_Z_d_zeta * d_Y_d_s) / sqrt_g
    grad_theta_vmec_Y = (d_Z_d_zeta * d_X_d_s - d_X_d_zeta * d_Z_d_s) / sqrt_g
    grad_theta_vmec_Z = (d_X_d_zeta * d_Y_d_s - d_Y_d_zeta * d_X_d_s) / sqrt_g

    grad_zeta_X = (d_Y_d_s * d_Z_d_theta_vmec - d_Z_d_s * d_Y_d_theta_vmec) / sqrt_g
    grad_zeta_Y = (d_Z_d_s * d_X_d_theta_vmec - d_X_d_s * d_Z_d_theta_vmec) / sqrt_g
    grad_zeta_Z = (d_X_d_s * d_Y_d_theta_vmec - d_Y_d_s * d_X_d_theta_vmec) / sqrt_g
    ! End of the dual relations.

    ! Sanity check: grad_zeta_X should be -sin(zeta) / R:
    do izeta = -nzgrid,nzgrid
       temp2D(:,izeta) = -sin(zeta(izeta)) / R(:,izeta)
    end do
    call test_arrays(grad_zeta_X, temp2D, .false., 3.0e-2, 'grad_zeta_X')
    grad_zeta_X = temp2D ! We might as well use the exact value, which is in temp2D.

    ! Sanity check: grad_zeta_Y should be cos(zeta) / R:
    do izeta = -nzgrid,nzgrid
       temp2D(:,izeta) = cos(zeta(izeta)) / R(:,izeta)
    end do
    call test_arrays(grad_zeta_Y, temp2D, .false., 3.0e-2, 'grad_zeta_Y')
    grad_zeta_Y = temp2D ! We might as well use the exact value, which is in temp2D.

    ! grad_zeta_Z should be 0:
    call test_arrays(grad_zeta_Z, temp2D, .true., 1.0e-12, 'grad_zeta_Z')
    grad_zeta_Z = 0

    !*********************************************************************
    ! Compute the Cartesian components of other quantities we need:
    !*********************************************************************

    grad_psi_X = grad_s_X * edge_toroidal_flux_over_2pi
    grad_psi_Y = grad_s_Y * edge_toroidal_flux_over_2pi
    grad_psi_Z = grad_s_Z * edge_toroidal_flux_over_2pi

    ! Form grad alpha = grad (theta_vmec + Lambda - iota * zeta)
    do izeta = -nzgrid,nzgrid
       grad_alpha_X(:,izeta) = (d_Lambda_d_s(:,izeta) - zeta(izeta) * d_iota_d_s) * grad_s_X(:,izeta)
       grad_alpha_Y(:,izeta) = (d_Lambda_d_s(:,izeta) - zeta(izeta) * d_iota_d_s) * grad_s_Y(:,izeta)
       grad_alpha_Z(:,izeta) = (d_Lambda_d_s(:,izeta) - zeta(izeta) * d_iota_d_s) * grad_s_Z(:,izeta)
    end do
    grad_alpha_X = grad_alpha_X + (1 + d_Lambda_d_theta_vmec) * grad_theta_vmec_X + (-iota + d_Lambda_d_zeta) * grad_zeta_X
    grad_alpha_Y = grad_alpha_Y + (1 + d_Lambda_d_theta_vmec) * grad_theta_vmec_Y + (-iota + d_Lambda_d_zeta) * grad_zeta_Y
    grad_alpha_Z = grad_alpha_Z + (1 + d_Lambda_d_theta_vmec) * grad_theta_vmec_Z + (-iota + d_Lambda_d_zeta) * grad_zeta_Z

    grad_B_X = d_B_d_s * grad_s_X + d_B_d_theta_vmec * grad_theta_vmec_X + d_B_d_zeta * grad_zeta_X
    grad_B_Y = d_B_d_s * grad_s_Y + d_B_d_theta_vmec * grad_theta_vmec_Y + d_B_d_zeta * grad_zeta_Y
    grad_B_Z = d_B_d_s * grad_s_Z + d_B_d_theta_vmec * grad_theta_vmec_Z + d_B_d_zeta * grad_zeta_Z

    !temp2D = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec) * d_R_d_zeta + (iota - d_Lambda_d_zeta) * d_R_d_theta_vmec) / sqrt_g
!!$    do izeta = -nzgrid,nzgrid
!!$       !B_X(:,izeta) = temp2D(:,izeta) * cos(zeta(izeta))
!!$       !B_Y(:,izeta) = temp2D(:,izeta) * sin(zeta(izeta))
!!$       sin_angle = sin(zeta(izeta))
!!$       cos_angle = cos(zeta(izeta))
!!$       B_X(:,izeta) = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec(:,izeta)) * (d_R_d_zeta(:,izeta)*cos_angle - R(:,izeta)*sin_angle) &
!!$            + (iota - d_Lambda_d_zeta(:,izeta)) * d_R_d_theta_vmec(:,izeta)*cos_angle) / sqrt_g(:,izeta)
!!$       B_Y(:,izeta) = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec(:,izeta)) * (d_R_d_zeta(:,izeta)*sin_angle + R(:,izeta)*cos_angle) &
!!$            + (iota - d_Lambda_d_zeta(:,izeta)) * d_R_d_theta_vmec(:,izeta)*sin_angle) / sqrt_g(:,izeta)
!!$    end do
    B_X = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec) * d_X_d_zeta + (iota - d_Lambda_d_zeta) * d_X_d_theta_vmec) / sqrt_g
    B_Y = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec) * d_Y_d_zeta + (iota - d_Lambda_d_zeta) * d_Y_d_theta_vmec) / sqrt_g
    B_Z = edge_toroidal_flux_over_2pi * ((1 + d_Lambda_d_theta_vmec) * d_Z_d_zeta + (iota - d_Lambda_d_zeta) * d_Z_d_theta_vmec) / sqrt_g

    sqrt_s = sqrt(normalized_toroidal_flux_used)

    !*********************************************************************
    ! Sanity tests: Verify that the Jacobian equals the appropriate
    ! cross product of the basis vectors.
    !*********************************************************************

    temp2D = 0 &
         + d_X_d_s * d_Y_d_theta_vmec * d_Z_d_zeta &
         + d_Y_d_s * d_Z_d_theta_vmec * d_X_d_zeta &
         + d_Z_d_s * d_X_d_theta_vmec * d_Y_d_zeta &
         - d_Z_d_s * d_Y_d_theta_vmec * d_X_d_zeta &
         - d_X_d_s * d_Z_d_theta_vmec * d_Y_d_zeta &
         - d_Y_d_s * d_X_d_theta_vmec * d_Z_d_zeta
    call test_arrays(sqrt_g, temp2D, .false., 3.0e-2, 'sqrt_g')

    temp2D = 0 &
         + grad_s_X * grad_theta_vmec_Y * grad_zeta_Z &
         + grad_s_Y * grad_theta_vmec_Z * grad_zeta_X &
         + grad_s_Z * grad_theta_vmec_X * grad_zeta_Y &
         - grad_s_Z * grad_theta_vmec_Y * grad_zeta_X &
         - grad_s_X * grad_theta_vmec_Z * grad_zeta_Y &
         - grad_s_Y * grad_theta_vmec_X * grad_zeta_Z
    call test_arrays(1/sqrt_g, temp2D, .false., 3.0e-2, '1/sqrt_g')

    !*********************************************************************
    ! Sanity tests: Verify that
    ! \vec{B} dot (each of the covariant and contravariant basis vectors)
    ! matches the corresponding term from VMEC.
    !*********************************************************************

    call test_arrays(B_X * d_X_d_theta_vmec + B_Y * d_Y_d_theta_vmec + B_Z * d_Z_d_theta_vmec, B_sub_theta_vmec, .false., 3.0e-2, 'B_sub_theta_vmec')
    call test_arrays(B_X * d_X_d_s          + B_Y * d_Y_d_s          + B_Z * d_Z_d_s,          B_sub_s,          .false., 3.0e-2, 'B_sub_s')
    call test_arrays(B_X * d_X_d_zeta       + B_Y * d_Y_d_zeta       + B_Z * d_Z_d_zeta,       B_sub_zeta,       .false., 3.0e-2, 'B_sub_zeta')

    call test_arrays(B_X *          grad_s_X + B_Y *          grad_s_Y + B_Z *          grad_s_Z,           temp2D,  .true., 3.0e-2, 'B_sup_s')
    call test_arrays(B_X *       grad_zeta_X + B_Y *       grad_zeta_Y + B_Z *       grad_zeta_Z,       B_sup_zeta, .false., 3.0e-2, 'B_sup_zeta')
    call test_arrays(B_X * grad_theta_vmec_X + B_Y * grad_theta_vmec_Y + B_Z * grad_theta_vmec_Z, B_sup_theta_vmec, .false., 3.0e-2, 'B_sup_theta_vmec')

    !*********************************************************************
    ! For gbdrift, we need \vect{B} cross grad |B| dot grad alpha.
    ! For cvdrift, we also need \vect{B} cross grad s dot grad alpha.
    ! Let us compute both of these quantities 2 ways, and make sure the two
    ! approaches give the same answer (within some tolerance).
    !*********************************************************************

    B_cross_grad_s_dot_grad_alpha = (B_sub_zeta * (1 + d_Lambda_d_theta_vmec) &
         - B_sub_theta_vmec * (d_Lambda_d_zeta - iota) ) / sqrt_g

    B_cross_grad_s_dot_grad_alpha_alternate = 0 &
         + B_X * grad_s_Y * grad_alpha_Z &
         + B_Y * grad_s_Z * grad_alpha_X &
         + B_Z * grad_s_X * grad_alpha_Y &
         - B_Z * grad_s_Y * grad_alpha_X &
         - B_X * grad_s_Z * grad_alpha_Y &
         - B_Y * grad_s_X * grad_alpha_Z

    call test_arrays(B_cross_grad_s_dot_grad_alpha, B_cross_grad_s_dot_grad_alpha_alternate, &
         .false., 3.0e-2, 'B_cross_grad_s_dot_grad_alpha')

    do izeta = -nzgrid,nzgrid
       B_cross_grad_B_dot_grad_alpha(:,izeta) = 0 &
            + (B_sub_s(:,izeta) * d_B_d_theta_vmec(:,izeta) * (d_Lambda_d_zeta(:,izeta) - iota) &
            + B_sub_theta_vmec(:,izeta) * d_B_d_zeta(:,izeta) * (d_Lambda_d_s(:,izeta) - zeta(izeta) * d_iota_d_s) &
            + B_sub_zeta(:,izeta) * d_B_d_s(:,izeta) * (1 + d_Lambda_d_theta_vmec(:,izeta)) &
            - B_sub_zeta(:,izeta) * d_B_d_theta_vmec(:,izeta) * (d_Lambda_d_s(:,izeta) - zeta(izeta) * d_iota_d_s) &
            - B_sub_theta_vmec(:,izeta) * d_B_d_s(:,izeta) * (d_Lambda_d_zeta(:,izeta) - iota) &
            - B_sub_s(:,izeta) * d_B_d_zeta(:,izeta) * (1 + d_Lambda_d_theta_vmec(:,izeta))) / sqrt_g(:,izeta)
    end do

    B_cross_grad_B_dot_grad_alpha_alternate = 0 &
         + B_X * grad_B_Y * grad_alpha_Z &
         + B_Y * grad_B_Z * grad_alpha_X &
         + B_Z * grad_B_X * grad_alpha_Y &
         - B_Z * grad_B_Y * grad_alpha_X &
         - B_X * grad_B_Z * grad_alpha_Y &
         - B_Y * grad_B_X * grad_alpha_Z

    call test_arrays(B_cross_grad_B_dot_grad_alpha, B_cross_grad_B_dot_grad_alpha_alternate, &
         .false., 3.0e-2, 'B_cross_grad_B_dot_grad_alpha')

    !*********************************************************************
    ! Finally, assemble the quantities needed for gs2.
    !*********************************************************************

    ! See the latex note gs2_full_surface_stellarator_geometry in the "doc" directory for a derivation of the formulae that follow.

    bmag = B / B_reference

    gradpar = L_reference * B_sup_zeta / B

    gds2 = (grad_alpha_X * grad_alpha_X + grad_alpha_Y * grad_alpha_Y + grad_alpha_Z * grad_alpha_Z) &
         * L_reference * L_reference * normalized_toroidal_flux_used

    gds21 = (grad_alpha_X * grad_psi_X + grad_alpha_Y * grad_psi_Y + grad_alpha_Z * grad_psi_Z) &
         * shat / B_reference

    gds22 = (grad_psi_X * grad_psi_X + grad_psi_Y * grad_psi_Y + grad_psi_Z * grad_psi_Z) &
         * shat * shat / (L_reference * L_reference * B_reference * B_reference * normalized_toroidal_flux_used)

    gbdrift = 2 * B_reference * L_reference * L_reference * sqrt_s * B_cross_grad_B_dot_grad_alpha &
         / (B * B * B)

    gbdrift0 = (B_sub_theta_vmec * d_B_d_zeta - B_sub_zeta * d_B_d_theta_vmec) / sqrt_g * edge_toroidal_flux_over_2pi &
         * 2 * shat / (B * B * B * sqrt_s)
    ! In the above 2-line expression for gbdrift0, the first line is \vec{B} \times \nabla B \cdot \nabla \psi.

    cvdrift = gbdrift + 2 * B_reference * L_reference * L_reference * sqrt_s * mu_0 * d_pressure_d_s &
         * B_cross_grad_s_dot_grad_alpha / (B * B * B * B)

    cvdrift0 = gbdrift0

    !*********************************************************************
    ! Free all arrays that were allocated.
    !*********************************************************************

    deallocate(B)
    deallocate(temp2D)
    deallocate(sqrt_g)
    deallocate(R)
    deallocate(d_B_d_theta_vmec)
    deallocate(d_B_d_zeta)
    deallocate(d_B_d_s)
    deallocate(d_R_d_theta_vmec)
    deallocate(d_R_d_zeta)
    deallocate(d_R_d_s)
    deallocate(d_Z_d_theta_vmec)
    deallocate(d_Z_d_zeta)
    deallocate(d_Z_d_s)
    deallocate(d_Lambda_d_theta_vmec)
    deallocate(d_Lambda_d_zeta)
    deallocate(d_Lambda_d_s)
    deallocate(B_sub_s)
    deallocate(B_sub_theta_vmec)
    deallocate(B_sub_zeta)
    deallocate(B_sup_theta_vmec)
    deallocate(B_sup_zeta)

    deallocate(d_B_d_s_mnc)
    deallocate(d_B_d_s_mns)
    deallocate(d_R_d_s_mnc)
    deallocate(d_R_d_s_mns)
    deallocate(d_Z_d_s_mnc)
    deallocate(d_Z_d_s_mns)
    deallocate(d_Lambda_d_s_mnc)
    deallocate(d_Lambda_d_s_mns)

    deallocate(d_X_d_s)
    deallocate(d_X_d_theta_vmec)
    deallocate(d_X_d_zeta)
    deallocate(d_Y_d_s)
    deallocate(d_Y_d_theta_vmec)
    deallocate(d_Y_d_zeta)

    deallocate(grad_s_X)
    deallocate(grad_s_Y)
    deallocate(grad_s_Z)
    deallocate(grad_theta_vmec_X)
    deallocate(grad_theta_vmec_Y)
    deallocate(grad_theta_vmec_Z)
    deallocate(grad_zeta_X)
    deallocate(grad_zeta_Y)
    deallocate(grad_zeta_Z)
    deallocate(grad_psi_X)
    deallocate(grad_psi_Y)
    deallocate(grad_psi_Z)
    deallocate(grad_alpha_X)
    deallocate(grad_alpha_Y)
    deallocate(grad_alpha_Z)

    deallocate(B_X)
    deallocate(B_Y)
    deallocate(B_Z)
    deallocate(grad_B_X)
    deallocate(grad_B_Y)
    deallocate(grad_B_Z)
    deallocate(B_cross_grad_B_dot_grad_alpha)
    deallocate(B_cross_grad_B_dot_grad_alpha_alternate)
    deallocate(B_cross_grad_s_dot_grad_alpha)
    deallocate(B_cross_grad_s_dot_grad_alpha_alternate)

    deallocate(normalized_toroidal_flux_full_grid)
    deallocate(normalized_toroidal_flux_half_grid)
    deallocate(theta_vmec)

    if (verbose) print *,"Leaving vmec_to_gs2_geometry_interface."

  contains

    subroutine test_arrays(array1, array2, should_be_0, tolerance, name)
      ! This subroutine is used for verifying the geometry arrays.
      ! When should_be_0 = .true., the subroutine verifies that |array1| = 0 to within
      !     an absolute tolerance specified by 'tolerance'. array2 is ignored in this case.
      ! When should_be_0 = .false., the subroutine verifies that array1 = array2
      !     to within a relative tolerance specified by 'tolerance'.

      implicit none

      real, dimension(nalpha,-nzgrid:nzgrid) :: array1, array2
      real :: tolerance
      character(len=*) :: name
      logical :: should_be_0
      real :: max_value, max_difference

      if (should_be_0) then
         max_value = maxval(abs(array1))
         if (verbose) print *,"  maxval(abs(",trim(name),")):",max_value,"(should be << 1.)"
         if (max_value > tolerance) then
            print *,"Error! ",trim(name)," should be 0, but instead it is:"
            do ialpha = 1,nalpha
               print *,array1(ialpha,:)
            end do
            stop
         end if
      else
         max_difference = maxval(abs(array1 - array2)) / maxval(abs(array1) + abs(array2))
         if (verbose) print *,"  Relative difference between two methods for computing ",trim(name),":",max_difference,"(should be << 1.)"
         if (max_difference > tolerance) then
            print *,"Error! Two methods for computing ",trim(name)," disagree. Here comes method 1:"
            do ialpha = 1,nalpha
               !print *,array1(ialpha,:)
            end do
            print *,"Here comes method 2:"
            do ialpha = 1,nalpha
               !print *,array2(ialpha,:)
            end do
            print *,"Here comes the difference:"
            do ialpha = 1,nalpha
               !print *,array1(ialpha,:) - array2(ialpha,:)
            end do
            stop
         end if
      end if

    end subroutine test_arrays

  end subroutine vmec_to_gs2_geometry_interface

  ! --------------------------------------------------------------------------

  function fzero_residual(theta_vmec_try)

    use read_wout_mod, only: xm, xn, mnmax, lmns, lmnc, lasym
    ! Note that lmns and lmnc use the non-Nyquist xm, xn, and mnmax.
    ! Also note that lmns and lmnc are on the HALF grid.

    implicit none

    real :: theta_vmec_try, fzero_residual
    real :: angle, sinangle, cosangle
    integer :: imn, which_surface

    ! residual = (theta_pest based on theta_vmec_try) - theta_pest_target = theta_vmec_try + Lambda - theta_pest_target
    fzero_residual = theta_vmec_try - theta_pest_target

    do imn = 1, mnmax
       angle = xm(imn)*theta_vmec_try - xn(imn)*zeta0
       sinangle = sin(angle)
       cosangle = cos(angle)
       do which_surface = 1,2
          fzero_residual = fzero_residual + vmec_radial_weight_half(which_surface) * lmns(imn,vmec_radial_index_half(which_surface)) * sinangle
          if (lasym) then
             fzero_residual = fzero_residual + vmec_radial_weight_half(which_surface) * lmnc(imn,vmec_radial_index_half(which_surface)) * cosangle
          end if
       end do
    end do

  end function fzero_residual

end module vmec_to_gs2_geometry_interface_mod

max_stars_repo_name: rogeriojorge/NearAxisGyrokinetics

max_issues_repo_name: rogeriojorge/NearAxisGyrokinetics

max_forks_repo_name: rogeriojorge/NearAxisGyrokinetics