diff --git a/src/gray.f b/src/gray.f index ae70b6c..6983c6f 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1204,15 +1204,12 @@ c local variables real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,current, . dpsinr,psia0,fp,xb,xe,ssfp,psinxp,rmop,zmop,psinoptmp,r10,z10, . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax, - . rhot15,psi15,rhot2,psi2 + . rhot15,psi15,rhot2,psi2,rhotsx real(wp_), dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol, . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim character(len=48) :: stringa c common/external functions/variables - real(wp_) :: rhotsx real(wp_) :: frhotor -c - common/rhotsx/rhotsx c c read from file eqdsk c (see http://fusion.gat.com/efit/g_eqdsk.html) @@ -1585,7 +1582,7 @@ c compute B_toroidal on axis print'(a,f8.4)','BT_axis = ',btaxis c compute normalized rho_tor from eqdsk q profile - call rhotor(nr) + call rhotor(rhotsx) phitedge=abs(psia)*rhotsx*2*pi c rrtor=sqrt(phitedge/abs(btrcen)/pi) call rhopol @@ -2065,21 +2062,17 @@ c c c c - subroutine rhotor(nnr) + subroutine rhotor(rhotsx) use const_and_precisions, only : wp_ - use interp_eqprof, only : nr,psinr,qpsi,crhot,cq - use simplespline, only : difcsn + use interp_eqprof, only : nr,psinr,rhopnr,qpsi,crhot,cq + use simplespline, only : difcs implicit none c arguments - integer :: nnr + real(wp_), intent(out) :: rhotsx c local variables integer :: iopt,ier,k real(wp_) :: dx,drhot real(wp_), dimension(nr) :: rhotnr -c common/external functions/variables - real(wp_) :: rhotsx -c - common/rhotsx/rhotsx c if(allocated(cq)) deallocate(cq) if(allocated(crhot)) deallocate(crhot) @@ -2088,25 +2081,25 @@ c c normalized toroidal rho : ~ Integral q(psi) dpsi c iopt=0 - call difcsn(psinr,qpsi,nr,nnr,iopt,cq,ier) + call difcs(psinr,qpsi,nr,iopt,cq,ier) c rhotnr(1)=0.0_wp_ - do k=1,nnr-1 + do k=1,nr-1 dx=psinr(k+1)-psinr(k) drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0_wp_+dx*(cq(k,3)/3.0_wp_ . +dx*cq(k,4)/4.0_wp_))) rhotnr(k+1)=rhotnr(k)+drhot end do - rhotsx=rhotnr(nnr) + rhotsx=rhotnr(nr) do k=1,nr - rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nnr)) + rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr)) end do c c spline interpolation of rhotor c iopt=0 - call difcsn(psinr,rhotnr,nr,nnr,iopt,crhot,ier) -c + rhopnr=sqrt(psinr) + call difcs(rhopnr,rhotnr,nr,iopt,crhot,ier) return end c @@ -2134,23 +2127,22 @@ c c c c - function frhotor_eq(psi) + function frhotor_eq(rhop) use const_and_precisions, only : wp_ - use interp_eqprof, only : psinr,nr,crhot - use simplespline, only :spli + use interp_eqprof, only : rhopnr,nr,crhot + use simplespline, only : spli + use utils, only : locate implicit none c arguments - real(wp_) :: psi,frhotor_eq + real(wp_) :: rhop,frhotor_eq c local variables integer :: irt - real(wp_) :: dps + real(wp_) :: drh c - irt=int((nr-1)*psi+1) -c if(irt.eq.0) irt=1 -c if(irt.eq.nr) irt=nr-1 + call locate(rhopnr,nr,rhop,irt) irt=min(max(1,irt),nr-1) - dps=psi-psinr(irt) - frhotor_eq=spli(crhot,nr,irt,dps) + drh=rhop-rhopnr(irt) + frhotor_eq=spli(crhot,nr,irt,drh) c return end @@ -2166,7 +2158,7 @@ c arguments c common/external functions/variables real(wp_) :: frhotor_eq,frhotor_an c - if(iequil.eq.2) frhotor=frhotor_eq(psi) + if(iequil.eq.2) frhotor=frhotor_eq(sqrt(psi)) if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) c return @@ -2779,9 +2771,9 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qqv(jp)=qq - dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) + dadrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2) . *dadpsi/pi - dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) + dvdrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2) . *dvdpsi/pi c computation of rhot from calculated q profile @@ -2846,7 +2838,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ rpstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_ end if - rhot_eq(jp)=frhotor_eq(pstab(jp)) + rhot_eq(jp)=frhotor_eq(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), @@ -3410,6 +3402,7 @@ c function frhotor_an(rhop) use const_and_precisions, only : wp_ use dierckx, only : splev + use utils, only : locate implicit none c local contants integer, parameter :: nnr=101,nrest=nnr+4 diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90 index 549142f..babfa50 100644 --- a/src/interp_eqprof.f90 +++ b/src/interp_eqprof.f90 @@ -10,7 +10,7 @@ module interp_eqprof REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz,tfp REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,cfp REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq01,cceq10,cceq20,cceq02,cceq11 - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,qpsi,rv,zv + REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopnr,qpsi,rv,zv REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin,psi,btotal,crhot,cq REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim,rbbbs,zbbbs @@ -43,8 +43,8 @@ contains call dealloc_equilvec allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), & btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), & - cceq20(lw20),cceq11(lw11),psi(nr,nz),psin(nr,nz),psinr(nr),qpsi(nr), & - stat=ier) + cceq20(lw20),cceq11(lw11),psi(nr,nz),psin(nr,nz),psinr(nr),rhopnr(nr), & + qpsi(nr), stat=ier) if (ier/=0) call dealloc_equilvec end subroutine alloc_equilvec @@ -66,6 +66,7 @@ contains if(allocated(psi)) deallocate(psi) if(allocated(psin)) deallocate(psin) if(allocated(psinr)) deallocate(psinr) + if(allocated(rhopnr)) deallocate(rhopnr) if(allocated(qpsi)) deallocate(qpsi) end subroutine dealloc_equilvec