VMEC_ops.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
module vmec_ops
#include <PB3D_macros.h>
    use str_utilities
    use output_ops
    use messages
    use num_vars, only: &
        &dp, max_str_ln, pi
    use read_wout_mod, only: read_wout_file, read_wout_deallocate, &            ! from LIBSTELL
        &n_r_vmec => ns, &                                                      ! n_r
        &xn, xm, &                                                              ! xn, xm
        &lrfp, &                                                                ! whether or not the poloidal flux is used as radial variable
        &phi, phipf, &                                                          ! toroidal flux (FM), norm. deriv. of toroidal flux (FM)
        &iotaf, &                                                               ! rot. transf. (tor. flux) or saf. fac. (pol. flux) (FM)
        &presf, &                                                               ! pressure (FM)
        &bsubumns, bsubumnc, bsubvmns, bsubvmnc, bsubsmns, bsubsmnc, &          ! B_theta (HM), B_zeta (HM), B_r (FM)
        &jcuru, jcurv, &                                                        ! poloidal (FM) and toroidal (FM) current
        &bmns, bmnc, &                                                          ! magnitude of B (HM)
        &lmns, lmnc, rmns, rmnc, zmns, zmnc, &                                  ! lambda (HM), R (FM), Z(FM)
        &gmnc, gmns                                                             ! Jacobian in VMEC coordinates
    use vmec_vars
    
    implicit none
    private
    public read_vmec, normalize_vmec
 
contains
    integer function read_vmec(n_r_in,use_pol_flux_V) result(ierr)
        use num_utilities, only: calc_int, spline
        use num_vars, only: eq_name, max_deriv, norm_disc_prec_eq
        
        character(*), parameter :: rout_name = 'read_VMEC'
        
        ! input / output
        integer, intent(inout) :: n_r_in
        logical, intent(inout) :: use_pol_flux_v
        
        ! local variables
        integer :: id, kd                                                       ! counters
        real(dp), allocatable :: r_v(:)                                         ! normal coordinate
        real(dp), allocatable :: l_c_h(:,:,:)                                   ! temporary HM variable
        real(dp), allocatable :: l_s_h(:,:,:)                                   ! temporary HM variable
        real(dp), allocatable :: jac_c_h(:,:,:)                                 ! temporary HM variable
        real(dp), allocatable :: jac_s_h(:,:,:)                                 ! temporary HM variable
        character(len=max_str_ln) :: err_msg                                    ! error message
        character(len=8) :: flux_name                                           ! either poloidal or toroidal
        real(dp), allocatable :: b_v_sub_c_m(:,:,:), b_v_sub_s_m(:,:,:)         ! Coeff. of B_i in (co)sine series (r,theta,phi) (FM, HM, HM)
#if ldebug
        real(dp), allocatable :: b_v_c_h(:,:), b_v_s_h(:,:)                     ! Coeff. of magnitude of B (HM)
#endif
        
        ! initialize ierr
        ierr = 0
        
        call writo('Reading data from VMEC output "' &
            &// trim(eq_name) // '"')
        call lvl_ud(1)
        
        ! read VMEC output using LIBSTELL
        if (scan(eq_name, ".") .ne. scan(eq_name, ".", .true.)) &
            &call writo('There seems to be a dot "."  in "'//trim(eq_name)//&
            &'". This can create problems.', alert=.true.)
        call read_wout_file(eq_name,ierr)                                       ! read the VMEC file
        chckerr('Failed to read the VMEC file')
        
        ! set some variables
        n_r_in = n_r_vmec
        allocate(flux_t_v(n_r_in,0:max_deriv+1))
        allocate(flux_p_v(n_r_in,0:max_deriv+1))
        allocate(rot_t_v(n_r_in,0:max_deriv+1))
        allocate(q_saf_v(n_r_in,0:max_deriv+1))
        allocate(pres_v(n_r_in,0:max_deriv+1))
        use_pol_flux_v = lrfp
        flux_t_v(:,0) = phi
        flux_t_v(:,1) = phipf
        flux_p_v(:,1) = iotaf*phipf
        ierr = calc_int(flux_p_v(:,1),1.0_dp/(n_r_in-1.0_dp),flux_p_v(:,0))
        chckerr('')
        pres_v(:,0) = presf
        rot_t_v(:,0) = iotaf
        q_saf_v(:,0) = 1._dp/iotaf
        if (use_pol_flux_v) then
            flux_name = 'poloidal'
        else
            flux_name = 'toroidal'
        end if
        b_0_v = sum(1.5*bmnc(:,2)-0.5_dp*bmnc(:,3))                             ! convert to full mesh
        
        call writo('VMEC version is ' // trim(r2str(vmec_version)))
        if (is_freeb_v) then
            call writo('Free boundary VMEC')
        else
            call writo('Fixed boundary VMEC')
        end if
        if (is_asym_v) then
            call writo('No stellarator symmetry')
        else
            call writo('Stellarator symmetry')
        end if
        call writo('VMEC has '//trim(i2str(mpol_v))//' poloidal and '&
            &//trim(i2str(ntor_v))//' toroidal modes,')
        if (ntor_v.gt.0) call writo('(basic toroidal field period NFP = '//&
            &trim(i2str(nfp_v))//')')
        call writo('defined on '//trim(i2str(n_r_in))//' '//flux_name//&
            &' flux surfaces')
        
        call writo('Running tests')
        call lvl_ud(1)
        if (mnmax_v.ne.((ntor_v+1)+(2*ntor_v+1)*(mpol_v-1))) then               ! for mpol_V = 0, only ntor_V+1 values needed
            err_msg = 'Inconsistency in ntor_V, mpol_V and mnmax_V'
            ierr = 1
            chckerr(err_msg)
        end if
        call lvl_ud(-1)
        
        call writo('Updating variables')
        allocate(mn_v(mnmax_v,2))
        mn_v(:,1) = nint(xm(1:mnmax_v))
        mn_v(:,2) = nint(xn(1:mnmax_v))
        
        call lvl_ud(-1)
        call writo('Data from VMEC output successfully read')
        
        call writo('Saving VMEC output to PB3D')
        call lvl_ud(1)
        
        ! Allocate the Fourier coefficients
        allocate(r_v_c(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(r_v_s(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(z_v_c(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(z_v_s(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(l_v_c(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(l_v_s(mnmax_v,n_r_in,0:max_deriv+1))
        allocate(jac_v_c(mnmax_v,n_r_in,0:max_deriv))
        allocate(jac_v_s(mnmax_v,n_r_in,0:max_deriv))
        allocate(l_c_h(mnmax_v,n_r_in,0:0))
        allocate(l_s_h(mnmax_v,n_r_in,0:0))
        allocate(jac_c_h(mnmax_v,n_r_in,0:0))
        allocate(jac_s_h(mnmax_v,n_r_in,0:0))
        
        ! factors R_V_c,s; Z_V_c,s and L_C,s and HM varieties
        r_v_c(:,:,0) = rmnc(1:mnmax_v,:)
        z_v_s(:,:,0) = zmns(1:mnmax_v,:)
        l_s_h(:,:,0) = lmns(1:mnmax_v,:)
        jac_c_h(:,:,0) = gmnc(1:mnmax_v,:)
        if (is_asym_v) then                                                     ! following only needed in asymmetric situations
            r_v_s(:,:,0) = rmns(1:mnmax_v,:)
            z_v_c(:,:,0) = zmnc(1:mnmax_v,:)
            l_c_h(:,:,0) = lmnc(1:mnmax_v,:)
            jac_s_h(:,:,0) = gmns(1:mnmax_v,:)
        else
            r_v_s = 0._dp
            z_v_c = 0._dp
            l_c_h = 0._dp
            jac_s_h = 0._dp
        end if
        
        ! calculate data for normal derivatives  with the toroidal (or poloidal)
        ! flux, normalized wrt. to the  maximum flux, equidistantly, so the step
        ! size is 1/(n_r_in-1).
        allocate(r_v(n_r_in))
        r_v = [((kd-1._dp)/(n_r_in-1),kd=1,n_r_in)]
        
        ! up to order 2, full mesh
        do kd = 1,2
            ierr = spline(r_v,q_saf_v(:,0),r_v,q_saf_v(:,kd),&
                &ord=norm_disc_prec_eq,deriv=kd)
            chckerr('')
            ierr = spline(r_v,rot_t_v(:,0),r_v,rot_t_v(:,kd),&
                &ord=norm_disc_prec_eq,deriv=kd)
            chckerr('')
            ierr = spline(r_v,pres_v(:,0),r_v,pres_v(:,kd),&
                &ord=norm_disc_prec_eq,deriv=kd)
            chckerr('')
            ierr = spline(r_v,flux_t_v(:,0),r_v,flux_t_v(:,kd),&
                &ord=norm_disc_prec_eq,deriv=kd)
            chckerr('')
            ierr = spline(r_v,flux_p_v(:,0),r_v,flux_p_v(:,kd),&
                &ord=norm_disc_prec_eq,deriv=kd)
            chckerr('')
        end do
        
        ! up to order 2, half mesh: also extrapolate
        do kd = 0,2
            do id = 1,mnmax_v
                ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                    &jac_c_h(id,2:n_r_in,0),r_v,jac_v_c(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=kd,extrap=.true.)
                chckerr('')
                if (is_asym_v) then
                    ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                        &jac_s_h(id,2:n_r_in,0),r_v,jac_v_s(id,:,kd),&
                        &ord=norm_disc_prec_eq,deriv=kd,extrap=.true.)
                    chckerr('')
                end if
            end do
        end do
        
        ! up to order 2, full mesh
        do kd = 1,3
            do id = 1,mnmax_v
                ! full mesh
                ierr = spline(r_v,r_v_c(id,:,0),r_v,r_v_c(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=kd)
                chckerr('')
                ierr = spline(r_v,z_v_s(id,:,0),r_v,z_v_s(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=kd)
                chckerr('')
                if (is_asym_v) then
                    ierr = spline(r_v,r_v_s(id,:,0),r_v,r_v_s(id,:,kd),&
                        &ord=norm_disc_prec_eq,deriv=kd)
                    chckerr('')
                    ierr = spline(r_v,z_v_c(id,:,0),r_v,z_v_c(id,:,kd),&
                        &ord=norm_disc_prec_eq,deriv=kd)
                    chckerr('')
                end if
            end do
        end do
        
        ! up to order 2, half mesh: also extrapolate
        do kd = 0,3
            do id = 1,mnmax_v
                ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                    &l_s_h(id,2:n_r_in,0),r_v,l_v_s(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=kd,extrap=.true.)
                chckerr('')
                if (is_asym_v) then
                    ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                        &l_c_h(id,2:n_r_in,0),r_v,l_v_c(id,:,kd),&
                        &ord=norm_disc_prec_eq,deriv=kd,extrap=.true.)
                    chckerr('')
                end if
            end do
        end do
        
        ! These are only necessary if the compiler searches for uninitialized
        ! memory
        flux_t_v(:,3:) = 0._dp
        flux_p_v(:,3:) = 0._dp
        q_saf_v(:,3:) = 0._dp
        rot_t_v(:,3:) = 0._dp
        pres_v(:,3:) = 0._dp
        
        ! allocate helper variables
        allocate(b_v_sub_c_m(mnmax_v,n_r_in,3)); b_v_sub_c_m = 0._dp
        allocate(b_v_sub_s_m(mnmax_v,n_r_in,3)); b_v_sub_s_m = 0._dp
#if ldebug
        allocate(b_v_c_h(mnmax_v,n_r_in)); b_v_c_h = 0._dp
        allocate(b_v_s_h(mnmax_v,n_r_in)); b_v_s_h = 0._dp
#endif
        
        ! store in helper variables
        b_v_sub_s_m(:,:,1) = bsubsmns(1:mnmax_v,:)
        b_v_sub_c_m(:,:,2) = bsubumnc(1:mnmax_v,:)
        b_v_sub_c_m(:,:,3) = bsubvmnc(1:mnmax_v,:)
        if (is_asym_v) then                                                     ! following only needed in asymmetric situations
            b_v_sub_c_m(:,:,1) = bsubsmnc(1:mnmax_v,:)
            b_v_sub_s_m(:,:,2) = bsubumns(1:mnmax_v,:)
            b_v_sub_s_m(:,:,3) = bsubvmns(1:mnmax_v,:)
        else
            b_v_sub_c_m(:,:,1) = 0._dp
            b_v_sub_s_m(:,:,2) = 0._dp
            b_v_sub_s_m(:,:,3) = 0._dp
        end if
#if ldebug
        b_v_c_h(:,:) = bmnc(1:mnmax_v,:)
        if (is_asym_v) then                                                     ! following only needed in asymmetric situations
            b_v_s_h(:,:) = bmns(1:mnmax_v,:)
        else
            b_v_s_h(:,:) = 0._dp
        end if
#endif
        
        ! allocate FM variables
        allocate(b_v_sub_c(mnmax_v,n_r_in,3))
        allocate(b_v_sub_s(mnmax_v,n_r_in,3))
#if ldebug
        allocate(b_v_c(mnmax_v,n_r_in))
        allocate(b_v_s(mnmax_v,n_r_in))
        allocate(j_v_sup_int(n_r_in,2))
#endif
        
        ! half mesh: extrapolate
        do id = 1,mnmax_v
            do kd = 2,3
                ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                    &b_v_sub_c_m(id,2:n_r_in,kd),r_v,b_v_sub_c(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=0,extrap=.true.)
                chckerr('')
                ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                    &b_v_sub_s_m(id,2:n_r_in,kd),r_v,b_v_sub_s(id,:,kd),&
                    &ord=norm_disc_prec_eq,deriv=0,extrap=.true.)
                chckerr('')
            end do
        end do
#if ldebug
        do id = 1,mnmax_v
            ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                &b_v_c_h(id,2:n_r_in),r_v,b_v_c(id,:),&
                &ord=norm_disc_prec_eq,deriv=0,extrap=.true.)
            chckerr('')
            ierr = spline(-0.5_dp/n_r_in+r_v(2:n_r_in),&
                &b_v_s_h(id,2:n_r_in),r_v,b_v_s(id,:),&
                &ord=norm_disc_prec_eq,deriv=0,extrap=.true.)
            chckerr('')
        end do
        j_v_sup_int(:,1) = jcuru
        j_v_sup_int(:,2) = jcurv
#endif
        
        ! deallocate repacked variables
        call read_wout_deallocate()
        
        call lvl_ud(-1)
        call writo('Conversion complete')
    end function read_vmec
    
    subroutine normalize_vmec
        use eq_vars, only: pres_0, psi_0, r_0, b_0
        
        ! scale the VMEC quantities
        pres_v = pres_v/pres_0
        flux_t_v = flux_t_v/psi_0
        flux_p_v = flux_p_v/psi_0
        r_v_c = r_v_c/r_0
        r_v_s = r_v_s/r_0
        z_v_c = z_v_c/r_0
        z_v_s = z_v_s/r_0
        l_v_c = l_v_c
        l_v_s = l_v_s
        jac_v_c = jac_v_c/(r_0**3)
        jac_v_s = jac_v_s/(r_0**3)
        b_v_sub_s = b_v_sub_s/(r_0*b_0)
        b_v_sub_c = b_v_sub_c/(r_0*b_0)
#if ldebug
        b_v_c = b_v_c/b_0
        b_v_s = b_v_s/b_0
        j_v_sup_int = j_v_sup_int*psi_0/pres_0
#endif
    end subroutine normalize_vmec
end module vmec_ops