diff --git a/gpec/gpdiag.f b/gpec/gpdiag.f index e029e89f..9f1c7e3c 100644 --- a/gpec/gpdiag.f +++ b/gpec/gpdiag.f @@ -936,7 +936,7 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, $ llagbpar_mn,divx_mn,curv_mn,divx2_mn,curv2_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, - $ bvt_fun,bvz_fun,xmt_fun,xmz_fun,xvt_fun,xvz_fun,xmp1_fun + $ bvt_fun,bvz_fun,xmt_fun,xmp1_fun REAL(r8), DIMENSION(:,:), POINTER :: rs,zs,eqfunx,eqfuny,eqfuns COMPLEX(r8), DIMENSION(:,:), POINTER :: eulbparmns,lagbparmns, @@ -989,9 +989,6 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) CALL iscdftb(mfac,mpert,bvt_fun,mthsurf,bvt_mn) CALL iscdftb(mfac,mpert,bvz_fun,mthsurf,bvz_mn) CALL iscdftb(mfac,mpert,xmt_fun,mthsurf,xmt_mn) - CALL iscdftb(mfac,mpert,xmz_fun,mthsurf,xmz_mn) - CALL iscdftb(mfac,mpert,xvt_fun,mthsurf,xvt_mn) - CALL iscdftb(mfac,mpert,xvz_fun,mthsurf,xvz_mn) CALL iscdftb(mfac,mpert,xmp1_fun,mthsurf,xmp1_mn) xspmns(istep,:)=xsp_mn @@ -1010,12 +1007,12 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) CALL bicube_eval(rzphi,psifac(istep),theta(itheta),0) jac=rzphi%f(4) - cspl%fs(itheta,1)=xmt_fun(itheta)- - $ (chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) - cspl%fs(itheta,2)=xmz_fun(itheta)- - $ sq%f(4)*(chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) +c Form J*xi_perp directly from xi^psi and xi^alpha=xi^s/chi1. + cspl%fs(itheta,1)=chi1**2* + $ (-eqfun%f(2)*xsp_fun(itheta)+ + $ eqfun%f(3)*xms_fun(itheta)/chi1) + cspl%fs(itheta,2)=sq%f(4)*cspl%fs(itheta,1)- + $ jac*xms_fun(itheta)/chi1 cspl%fs(itheta,3)=xmt_fun(itheta) ENDDO @@ -1060,10 +1057,9 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) $ xsp_fun(itheta)*eqfun%fx(1)+ $ xms_fun(itheta)/(chi1*sq%f(4))*eqfun%fy(1) llagbparfun(istep,itheta)= - $ lagbparfun(istep,itheta)+eqfun%fy(1)* - $ (xmz_fun(itheta)/(jac*sq%f(4))- - $ (chi1/(jac*eqfun%f(1)))**2* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta))) + $ eulbparfun(istep,itheta)+ + $ xsp_fun(itheta)*eqfun%fx(1)+ + $ eqfun%fy(1)*cspl%f(1)/jac divxfun(istep,itheta)=xmp1_fun(itheta)+ $ (jac1/jac)*xsp_fun(itheta)+ @@ -1074,10 +1070,7 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) curvfun(istep,itheta)=-xsp_fun(itheta)/eqfun%f(1)*sq%f1(2)- $ xsp_fun(itheta)*eqfun%fx(1)-eqfun%fy(1)* - $ (xms_fun(itheta)/(chi1*sq%f(4))+ - $ xmz_fun(itheta)/(jac*sq%f(4))- - $ (chi1/(jac*eqfun%f(1)))**2* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta))) + $ cspl%f(1)/jac curv2fun(istep,itheta)=-jac*eqfun%f(1)*curvfun(istep,itheta) ENDDO diff --git a/gpec/gpout.f b/gpec/gpout.f index ce376b52..069ae193 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -2955,7 +2955,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, $ divxprp_mn,curv_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, - $ bvt_fun,bvz_fun,xmz_fun,xvt_fun,xvz_fun,xmt_fun,xsp1_fun + $ bvt_fun,bvz_fun,xsp1_fun COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmns, $ lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns @@ -3012,10 +3012,6 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL iscdftb(mfac,mpert,xms_fun,mthsurf,xms_mn) CALL iscdftb(mfac,mpert,bvt_fun,mthsurf,bvt_mn) CALL iscdftb(mfac,mpert,bvz_fun,mthsurf,bvz_mn) - CALL iscdftb(mfac,mpert,xmt_fun,mthsurf,xmt_mn) - CALL iscdftb(mfac,mpert,xmz_fun,mthsurf,xmz_mn) - CALL iscdftb(mfac,mpert,xvt_fun,mthsurf,xvt_mn) - CALL iscdftb(mfac,mpert,xvz_fun,mthsurf,xvz_mn) CALL spline_eval(sq,psifac(istep),1) @@ -3027,12 +3023,12 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL bicube_eval(rzphi,psifac(istep),theta(itheta),0) jac=rzphi%f(4) - cspl%fs(itheta,1)=xmt_fun(itheta)- - $ (chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) - cspl%fs(itheta,2)=xmz_fun(itheta)- - $ sq%f(4)*(chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) +c Form J*xi_perp directly from xi^psi and xi^alpha=xi^s/chi1. + cspl%fs(itheta,1)=chi1**2* + $ (-eqfun%f(2)*xsp_fun(itheta)+ + $ eqfun%f(3)*xms_fun(itheta)/chi1) + cspl%fs(itheta,2)=sq%f(4)*cspl%fs(itheta,1)- + $ jac*xms_fun(itheta)/chi1 ENDDO CALL cspline_fit(cspl,"periodic") @@ -3062,10 +3058,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ divxprpfun(istep,itheta) curvfun(istep,itheta)=-xsp_fun(itheta)/eqfun%f(1)*sq%f1(2)- $ (xsp_fun(itheta)*eqfun%fx(1)+eqfun%fy(1)* - $ (xms_fun(itheta)/(chi1*sq%f(4))+ - $ xmz_fun(itheta)/(jac*sq%f(4))- - $ (chi1/(jac*eqfun%f(1)))**2* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)))) + $ cspl%f(1)/jac) equilbfun(istep,itheta) = eqfun%f(1) dphi(itheta) = rzphi%f(3) ENDDO diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index e98d64a9..d764fffd 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -434,9 +434,8 @@ subroutine read_pmodb(file,jac_in,jsurf_in,tmag_in,debug,op_powin) +xj*reshape(table(:,10),(/npsi,nm/),order=(/1,2/)) endif - ! For consistency with GPEC-0.3.0 and smaller values near rationals - if(verbose) print *," -> Using div(xi_perp) = -(dB/B+kappa.xi_perp)" - divxmni = -(lagbmni+kapxmni) + ! Keep the supplied divergence and use its perpendicular gauge. + lagbmni = -(divxmni+kapxmni) ! convert to chebyshev coordinates if(jac_in=="" .or. jac_in=="default")then @@ -540,7 +539,7 @@ subroutine read_pmodb(file,jac_in,jsurf_in,tmag_in,debug,op_powin) write(out_unit,'(a8,es16.8e3,2/)') " chi1 = ",chi1 write(out_unit,'(1x,a16,a5,8(1x,a16))')"psi_n","m",& 'real(deltaB/B)','imag(deltaB/B)','real(divxprp)','imag(divxprp)',& - 'real(deltaB)','imag(deltaB)','real(Bdivxprp)','imag(Bdivxprp)' + 'real(deltaBprp)','imag(deltaBprp)','real(Bdivxprp)','imag(Bdivxprp)' do i=0,npsi-1 do j=1,mpert write(out_unit,'(1x,es16.8e3,i5,14(1x,es16.8e3))') &