module equilibrium use const_and_precisions, only : wp_ implicit none real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi real(wp_), save :: btrcen ! used only for jcd_astra def. real(wp_), save :: rcen ! computed as fpol(a)/btrcen real(wp_), save :: rmnm,rmxm,zmnm,zmxm real(wp_), save :: zbinf,zbsup ! real(wp_), save :: rup,zup,rlw,zlw integer, parameter :: kspl=3,ksplp=kspl+1 ! === 2d spline psi(r,z), normalization and derivatives ========== integer, save :: nsr, nsz real(wp_), save :: psia, psiant, psinop real(wp_), dimension(:), allocatable, save :: tr,tz real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, & cceq20, cceq02, cceq11 ! === 1d spline fpol(psi) ======================================== ! integer, save :: npsiest integer, save :: nsf real(wp_), dimension(:), allocatable, save :: tfp, cfp real(wp_), save :: fpolas ! === 1d spline rhot(rhop), rhop(rhot), q(psi) =================== ! computed on psinr,rhopnr [,rhotnr] arrays integer, save :: nq,nrho real(wp_), dimension(:), allocatable, save :: psinr,rhopr,rhotr real(wp_), dimension(:,:), allocatable, save :: cq,crhop,crhot real(wp_), save :: phitedge,aminor real(wp_), save :: q0,qa,alq contains subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, & rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit) use const_and_precisions, only : one use utils, only : get_free_unit implicit none ! arguments character(len=*), intent(in) :: filenm real(wp_), intent(out) :: psia,rvac,rax,zax real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim real(wp_), dimension(:,:), allocatable, intent(out) :: psin integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit ! local variables integer, parameter :: indef=0,iddef=1,iffdef=0 integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim character(len=48) :: string real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis ! set default values if optional arguments are absent in=indef; id=iddef; iffmt=iffdef if(present(ipsinorm)) in=ipsinorm if(present(idesc)) id=idesc if(present(ifreefmt)) iffmt=ifreefmt if (present(unit)) then u=unit else u=get_free_unit() end if ! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html) open(file=trim(filenm),status='old',action='read',unit=u) ! get size of main arrays and allocate them if (id==1) then read (u,'(a48,3i4)') string,idum,nr,nz else read (u,*) nr,nz end if if (allocated(rv)) deallocate(rv) if (allocated(zv)) deallocate(zv) if (allocated(psin)) deallocate(psin) if (allocated(psinr)) deallocate(psinr) if (allocated(fpol)) deallocate(fpol) if (allocated(q)) deallocate(q) allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr)) ! store 0D data and main arrays if (iffmt==1) then read (u,*) dr,dz,rvac,rleft,zmid read (u,*) rax,zax,psiaxis,psiedge,xdum read (u,*) xdum,xdum,xdum,xdum,xdum read (u,*) xdum,xdum,xdum,xdum,xdum read (u,*) (fpol(i),i=1,nr) read (u,*) (xdum,i=1,nr) read (u,*) (xdum,i=1,nr) read (u,*) (xdum,i=1,nr) read (u,*) ((psin(i,j),i=1,nr),j=1,nz) read (u,*) (q(i),i=1,nr) else read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum read (u,'(5e16.9)') (fpol(i),i=1,nr) read (u,'(5e16.9)') (xdum,i=1,nr) read (u,'(5e16.9)') (xdum,i=1,nr) read (u,'(5e16.9)') (xdum,i=1,nr) read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz) read (u,'(5e16.9)') (q(i),i=1,nr) end if ! get size of boundary and limiter arrays and allocate them read (u,*) nbnd,nlim if (allocated(rbnd)) deallocate(rbnd) if (allocated(zbnd)) deallocate(zbnd) if (allocated(rlim)) deallocate(rlim) if (allocated(zlim)) deallocate(zlim) ! store boundary and limiter data if(nbnd>0) then allocate(rbnd(nbnd),zbnd(nbnd)) if (iffmt==1) then read(u,*) (rbnd(i),zbnd(i),i=1,nbnd) else read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd) end if end if if(nlim>0) then allocate(rlim(nlim),zlim(nlim)) if (iffmt==1) then read(u,*) (rlim(i),zlim(i),i=1,nlim) else read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim) end if end if ! reading of G EQDSK file completed close(u) ! build rv,zv,psinr arrays and normalize psin zleft=zmid-0.5_wp_*dz dr=dr/(nr-1) dz=dz/(nz-1) dps=one/(nr-1) do i=1,nr psinr(i)=(i-1)*dps rv(i)=rleft+(i-1)*dr end do do i=1,nz zv(i)=zleft+(i-1)*dz end do psia=psiedge-psiaxis if(in==0) psin=(psin-psiaxis)/psia end subroutine read_eqdsk subroutine read_equil_an(filenm,rv,zv,fpol,q,unit) use utils, only : get_free_unit implicit none ! arguments character(len=*), intent(in) :: filenm integer, optional, intent(in) :: unit real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q ! local variables integer :: u real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen if (present(unit)) then u=unit else u=get_free_unit() end if open(file=trim(filenm),status='old',action='read',unit=u) read(u,*) rr0m,zr0m,rpam read(u,*) b0 read(u,*) q0,qa,alq ! rcen=rr0m ! btrcen=b0 ! b0=isgnbphi*b0*factb ! rvac=rr0m ! rax=rr0m ! zax=zr0m ! zbmin=(zr0-rpa)/1.0e2_wp_ ! zbmax=(zr0+rpa)/1.0e2_wp_ if(allocated(rv)) deallocate(rv) if(allocated(zv)) deallocate(zv) if(allocated(fpol)) deallocate(fpol) if(allocated(q)) deallocate(q) allocate(rv(2),zv(1),fpol(1),q(3)) rv(1)=rr0m rv(2)=rpam zv(1)=zr0m fpol(1)=b0*rr0m q(1)=q0 q(2)=qa q(3)=alq close(u) end subroutine read_equil_an subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr) use const_and_precisions, only : zero,one,pi implicit none ! arguments real(wp_), intent(inout) :: psia real(wp_), dimension(:), intent(inout) :: fpol,q integer, intent(in) :: cocosin, cocosout integer, intent(out), optional :: ierr ! local variables real(wp_) :: isign,bsign integer :: exp2pi,exp2piout logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos) call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout) ! check sign consistency isign=sign(one,psia) if (.not.psiincr) isign=-isign bsign=sign(one,fpol(size(fpol))) if (qpos.neqv.isign*bsign*q(size(q))>zero) then ! warning: sign inconsistency found among q, Ipla and Bref q=-q if(present(ierr)) ierr=1 else if(present(ierr)) ierr=0 end if ! convert cocosin to cocosout ! opposite direction of toroidal angle phi in cocosin and cocosout if (phiccw.neqv.phiccwout) fpol=-fpol ! q has opposite sign for given sign of Bphi*Ip if (qpos .neqv. qposout) q=-q ! psi and Ip signs don't change accordingly if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia ! convert Wb to Wb/rad or viceversa psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi) end subroutine change_cocos subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos) implicit none integer, intent(in) :: cocos integer, intent(out) :: exp2pi logical, intent(out) :: phiccw,psiincr,qpos integer :: cmod10,cmod4 cmod10=mod(cocos,10) cmod4=mod(cmod10,4) ! cocos>10 psi in Wb, cocos<10 psi in Wb/rad exp2pi=(cocos-cmod10)/10 ! cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW phiccw=(mod(cmod10,2)==1) ! cocos mod 10 = 1,2,5,6: psi increasing with positive Ip psiincr=(cmod4==1 .or. cmod4==2) ! cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip qpos=(cmod10<3 .or. cmod10>6) end subroutine decode_cocos subroutine eq_scal(psia,fpol,isign,bsign,factb) use const_and_precisions, only : one implicit none real(wp_), intent(inout) :: psia integer, intent(inout) :: isign,bsign real(wp_), dimension(:), intent(inout) :: fpol real(wp_), intent(in) :: factb ! isign and bsign ignored on input if equal to 0 ! used to assign the direction of Bphi and Ipla BEFORE scaling ! cocos=3 assumed: CCW direction is >0 ! Bphi and Ipla scaled by the same factor factb to keep q unchanged ! factb<0 reverses the directions of Bphi and Ipla if(isign/=0) psia=sign(psia,real(-isign,wp_)) if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol psia=psia*factb fpol=fpol*factb isign=int(sign(one,-psia)) bsign=int(sign(one,fpol(size(fpol)))) end subroutine eq_scal subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,sspl,ssfp, & r0,rax,zax,rbnd,zbnd,ixp) use const_and_precisions, only : zero,one use dierckx, only : regrid,coeff_parder,curfit,splev use utils, only : vmaxmin,vmaxmini implicit none ! arguments real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol real(wp_), dimension(:,:), intent(in) :: psin real(wp_), intent(in) :: psiwbrad real(wp_), intent(in) :: sspl,ssfp real(wp_), intent(in), optional :: r0,rax,zax real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd integer, intent(in), optional :: ixp ! local constants integer, parameter :: iopt=0 ! local variables integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 real(wp_), dimension(1) :: fpoli real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk integer :: ier,ixploc,info ! compute array sizes and prepare working space arrays nr=size(rv) nz=size(zv) nrz=nr*nz nrest=nr+ksplp nzest=nz+ksplp lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest) liwrk=nz+nr+nrest+nzest+kspl npsi=size(psinr) npsest=npsi+ksplp lwrkf=npsi*ksplp+npsest*(7+3*kspl) allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest))) ! spline fitting/interpolation of psin(i,j) and derivatives ! length in m !!! rmnm=rv(1) rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) ! allocate knots and spline coefficients arrays if (allocated(tr)) deallocate(tr) if (allocated(tz)) deallocate(tz) allocate(tr(nrest),tz(nzest),cceq(nrz)) ! allocate work arrays ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs allocate(fvpsi(nrz)) fvpsi=reshape(transpose(psin),(/nrz/)) sspln=sspl call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) ! if ier=-1 data are re-fitted using sspl=0 if(ier==-1) then sspln=0.0_wp_ call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) end if deallocate(fvpsi) ! compute spline coefficients for psi partial derivatives lw10 = nr*(ksplp-1) + nz*ksplp + nrz lw01 = nr*ksplp + nz*(ksplp-1) + nrz lw20 = nr*(ksplp-2) + nz*ksplp + nrz lw02 = nr*ksplp + nz*(ksplp-2) + nrz lw11 = nr*(ksplp-1) + nz*(ksplp-1) + nrz if (allocated(cceq10)) deallocate(cceq10) if (allocated(cceq01)) deallocate(cceq01) if (allocated(cceq20)) deallocate(cceq20) if (allocated(cceq02)) deallocate(cceq02) if (allocated(cceq11)) deallocate(cceq11) allocate(cceq10(lw10),cceq01(lw01),cceq20(lw20),cceq02(lw02),cceq11(lw11)) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,0,cceq10,lw10,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,1,cceq01,lw01,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier) ! spline interpolation of fpol(i) ! allocate knots and spline coefficients arrays if (allocated(tfp)) deallocate(tfp) if (allocated(cfp)) deallocate(cfp) allocate(tfp(npsest),cfp(npsest)) allocate(wf(npsi)) wf(1:npsi-1)=one wf(npsi)=1.0e2_wp_ call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, & tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier) call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier) ! set vacuum value used outside 0<=psin<=1 range fpolas=fpoli(1) sgnbphi=sign(one,fpolas) ! free temporary arrays deallocate(wrk,iwrk,wf) ! re-normalize psi after spline computation ! start with un-corrected psi psia=psiwbrad psinop=0.0_wp_ psiant=1.0_wp_ ! use provided boundary to set an initial guess for the search of O/X points nbnd=0 if(present(rbnd).and.present(zbnd)) then nbnd=min(size(rbnd),size(zbnd)) end if if (nbnd>0) then call vmaxmini(zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) rbinf=rbnd(ibinf) rbsup=rbnd(ibsup) call vmaxmin(rbnd,nbnd,rbmin,rbmax) else zbinf=zv(2) zbsup=zv(nz-1) rbinf=rv((nr+1)/2) rbsup=rbinf rbmin=rv(2) rbmax=rv(nr-1) end if ! search for exact location of the magnetic axis if(present(rax)) then rax0=rax else rax0=0.5_wp_*(rbmin+rbmax) end if if(present(zax)) then zax0=zax else zax0=0.5_wp_*(zbinf+zbsup) end if call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info) print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp ! search for X-point if ixp not = 0 if(present(ixp)) then ixploc=ixp else ixploc=0 end if if(ixploc/=0) then if(ixploc<0) then call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) if(psinxptmp/=-1.0_wp_) then print'(a,2f8.4,es12.5)','X-point',r1,z1,psinxptmp zbinf=z1 psinop=psinoptmp psiant=psinxptmp-psinop call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) zbsup=z1 else ixploc=0 end if else call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info) if(psinxptmp.ne.-1.0_wp_) then print'(a,2f8.4,e16.8)','X-point',r1,z1,psinxptmp zbsup=z1 psinop=psinoptmp psiant=psinxptmp-psinop call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) zbinf=z1 else ixploc=0 end if end if end if if (ixploc==0) then psinop=psinoptmp psiant=one-psinop ! find upper horizontal tangent point call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) zbsup=z1 rbsup=r1 ! find lower horizontal tangent point call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) zbinf=z1 rbinf=r1 print'(a,4f8.4)','no X-point ',rbinf,zbinf,rbsup,zbsup end if print*,' ' ! save Bt value on axis (required in flux_average and used in Jcd def) ! and vacuum value B0 at ref. radius R0 (used in Jcd_astra def) call equinum_fpol(0.0_wp_,btaxis) btaxis=btaxis/rmaxis if(present(r0)) then btrcen=fpolas/r0 rcen=r0 else btrcen=fpolas/rmaxis rcen=rmaxis end if print'(a,f8.4)','BT_centr= ',btrcen print'(a,f8.4)','BT_axis = ',btaxis end subroutine set_eqspl subroutine unset_eqspl implicit none if(allocated(tr)) deallocate(tr) if(allocated(tz)) deallocate(tz) if(allocated(tfp)) deallocate(tfp) if(allocated(cfp)) deallocate(cfp) if(allocated(cceq)) deallocate(cceq) if(allocated(cceq01)) deallocate(cceq01) if(allocated(cceq10)) deallocate(cceq10) if(allocated(cceq02)) deallocate(cceq02) if(allocated(cceq20)) deallocate(cceq20) if(allocated(cceq11)) deallocate(cceq11) nsr=0 nsz=0 nsf=0 end subroutine unset_eqspl subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) use dierckx, only : fpbisp implicit none ! local constants integer, parameter :: lwrk=8,liwrk=2 ! arguments real(wp_), intent(in) :: rpsim,zpsim real(wp_), intent(out), optional :: psinv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz ! local variables integer, dimension(liwrk) :: iwrk real(wp_), dimension(1) :: rrs,zzs,ffspl real(wp_), dimension(lwrk) :: wrk ! here lengths are measured in meters if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. & zpsim.le.zmxm .and. zpsim.ge.zmnm) then if (present(psinv)) then rrs(1)=rpsim zzs(1)=zpsim call fpbisp(tr,nsr,tz,nsz,cceq,3,3,rrs,1,zzs,1,ffspl, & wrk(1),wrk(5),iwrk(1),iwrk(2)) psinv=(ffspl(1)-psinop)/psiant end if if (present(dpsidr)) then call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10) end if if (present(dpsidz)) then call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01) end if if (present(ddpsidrr)) then call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20) end if if (present(ddpsidzz)) then call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02) end if if (present(ddpsidrz)) then call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11) end if else if(present(psinv)) psinv=-1.0_wp_ if(present(dpsidr)) dpsidr=0.0_wp_ if(present(dpsidz)) dpsidz=0.0_wp_ if(present(ddpsidrr)) ddpsidrr=0.0_wp_ if(present(ddpsidzz)) ddpsidzz=0.0_wp_ if(present(ddpsidrz)) ddpsidrz=0.0_wp_ end if end subroutine equinum_psi subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc) use dierckx, only : fpbisp implicit none ! arguments real(wp_), intent(in) :: rpsim,zpsim integer, intent(in) :: nur,nuz real(wp_), intent(out) :: derpsi real(wp_), dimension(:) :: cc ! local variables integer, dimension(1) :: iwrkr,iwrkz real(wp_), dimension(1) :: rrs,zzs,ffspl real(wp_), dimension(1,ksplp) :: wrkr real(wp_), dimension(1,ksplp) :: wrkz rrs(1)=rpsim zzs(1)=zpsim call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kspl-nur,kspl-nuz, & rrs,1,zzs,1,ffspl,wrkr,wrkz,iwrkr,iwrkz) derpsi=ffspl(1)*psia end subroutine sub_derpsi subroutine equinum_fpol(psinv,fpolv,dfpv) use dierckx, only : splev,splder implicit none ! arguments real(wp_), intent(in) :: psinv real(wp_), intent(out) :: fpolv real(wp_), intent(out), optional :: dfpv ! local variables integer :: ier real(wp_), dimension(1) :: rrs,ffspl real(wp_), dimension(nsf) :: wrkfd ! if(psinv.le.1.0_wp_.and.psinv.ge.0.0_wp_) then rrs(1)=psinv call splev(tfp,nsf,cfp,3,rrs,ffspl,1,ier) fpolv=ffspl(1) if(present(dfpv)) then call splder(tfp,nsf,cfp,3,1,rrs,ffspl,1,wrkfd,ier) dfpv=ffspl(1)/psia end if else fpolv=fpolas if (present(dfpv)) dfpv=0._wp_ end if end subroutine equinum_fpol subroutine bfield(rpsim,zpsim,bphi,br,bz) use gray_params, only : iequil implicit none ! arguments real(wp_), intent(in) :: rpsim,zpsim real(wp_), intent(out), optional :: bphi,br,bz ! local variables real(wp_) :: psin,fpol if (iequil < 2) then call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) bphi=bphi/rpsim else call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) then psin=bphi call equinum_fpol(psin,fpol) bphi=fpol/rpsim end if end if if (present(br)) br=-br/rpsim if (present(bz)) bz= bz/rpsim end subroutine bfield subroutine setqphi_num(psinq,q,psia,rhotn) use const_and_precisions, only : pi use simplespline, only : difcs implicit none ! arguments real(wp_), dimension(:), intent(in) :: psinq,q real(wp_), intent(in) :: psia real(wp_), dimension(:), intent(out), optional :: rhotn ! local variables real(wp_), dimension(size(q)) :: phit real(wp_) :: dx integer, parameter :: iopt=0 integer :: k,ier nq=size(q) if(allocated(psinr)) deallocate(psinr) if(allocated(cq)) deallocate(cq) allocate(psinr(nq),cq(nq,4)) psinr=psinq call difcs(psinr,q,nq,iopt,cq,ier) ! Toroidal flux phi = 2*pi*Integral q dpsi phit(1)=0.0_wp_ do k=1,nq-1 dx=psinr(k+1)-psinr(k) phit(k+1)=phit(k) + 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_) ) ) end do phitedge=phit(nq) if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge) phitedge=2*pi*psia*phitedge end subroutine setqphi_num subroutine unset_q implicit none if(allocated(psinr)) deallocate(psinr) if(allocated(cq)) deallocate(cq) nq=0 end subroutine unset_q function fq(psin) use const_and_precisions, only : wp_ use gray_params, only : iequil use simplespline, only :spli use utils, only : locate implicit none ! arguments real(wp_), intent(in) :: psin real(wp_) :: fq ! local variables integer :: i real(wp_) :: dps,rn if (iequil<2) then rn=frhotor(sqrt(psin)) fq=q0+(qa-q0)*rn**alq else call locate(psinr,nq,psin,i) i=min(max(1,i),nq-1) dps=psin-psinr(i) fq=spli(cq,nq,i,dps) end if end function fq subroutine set_rhospl(rhop,rhot) use simplespline, only : difcs implicit none ! arguments real(wp_), dimension(:), intent(in) :: rhop, rhot ! local variables integer, parameter :: iopt=0 integer :: ier nrho=size(rhop) if(allocated(rhopr)) deallocate(rhopr) if(allocated(rhotr)) deallocate(rhotr) if(allocated(crhop)) deallocate(crhop) if(allocated(crhot)) deallocate(crhot) allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4)) rhopr=rhop rhotr=rhot call difcs(rhotr,rhopr,nrho,iopt,crhop,ier) call difcs(rhopr,rhotr,nrho,iopt,crhot,ier) end subroutine set_rhospl subroutine unset_rhospl implicit none if(allocated(rhopr)) deallocate(rhopr) if(allocated(rhotr)) deallocate(rhotr) if(allocated(crhop)) deallocate(crhop) if(allocated(crhot)) deallocate(crhot) nrho=0 end subroutine unset_rhospl function frhopol(rhot) use utils, only : locate use simplespline, only : spli implicit none ! arguments real(wp_), intent(in) :: rhot real(wp_) :: frhopol ! local variables integer :: i real(wp_) :: dr call locate(rhotr,nrho,rhot,i) i=min(max(1,i),nrho-1) dr=rhot-rhotr(i) frhopol=spli(crhop,nrho,i,dr) end function frhopol function frhopolv(rhot) use utils, only : locate use simplespline, only : spli implicit none ! arguments real(wp_), dimension(:), intent(in) :: rhot real(wp_), dimension(size(rhot)) :: frhopolv ! local variables integer :: i,i0,j real(wp_) :: dr i0=1 do j=1,size(rhot) call locate(rhotr(i0:),nrho-i0+1,rhot(j),i) i=min(max(1,i),nrho-i0)+i0-1 dr=rhot(j)-rhotr(i) frhopolv(j)=spli(crhop,nrho,i,dr) i0=i end do end function frhopolv function frhotor(rhop) use utils, only : locate use simplespline, only : spli implicit none ! arguments real(wp_), intent(in) :: rhop real(wp_) :: frhotor ! local variables integer :: i real(wp_) :: dr call locate(rhopr,nrho,rhop,i) i=min(max(1,i),nrho-1) dr=rhop-rhopr(i) frhotor=spli(crhot,nrho,i,dr) end function frhotor subroutine points_ox(rz,zz,rf,zf,psinvf,info) use const_and_precisions, only : comp_eps use minpack, only : hybrj1 implicit none ! local constants integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 ! arguments real(wp_), intent(in) :: rz,zz real(wp_), intent(out) :: rf,zf,psinvf integer, intent(out) :: info ! local variables real(wp_) :: tol real(wp_), dimension(n) :: xvec,fvec real(wp_), dimension(lwa) :: wa real(wp_), dimension(ldfjac,n) :: fjac xvec(1)=rz xvec(2)=zz tol = sqrt(comp_eps) call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then print'(a,i2,a,2f8.4)',' info subr points_ox =',info, & ' O/X coord.',xvec end if rf=xvec(1) zf=xvec(2) call equinum_psi(rf,zf,psinvf) end subroutine points_ox subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) implicit none ! arguments integer, intent(in) :: n,iflag,ldfjac real(wp_), dimension(n), intent(in) :: x real(wp_), dimension(n), intent(inout) :: fvec real(wp_), dimension(ldfjac,n), intent(inout) :: fjac ! local variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz select case(iflag) case(1) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) fvec(1) = dpsidr/psia fvec(2) = dpsidz/psia case(2) call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & ddpsidrz=ddpsidrz) fjac(1,1) = ddpsidrr/psia fjac(1,2) = ddpsidrz/psia fjac(2,1) = ddpsidrz/psia fjac(2,2) = ddpsidzz/psia case default print*,'iflag undefined' end select end subroutine fcnox subroutine points_tgo(rz,zz,rf,zf,psin0,info) use const_and_precisions, only : comp_eps use minpack, only : hybrj1mv implicit none ! local constants integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 ! arguments real(wp_), intent(in) :: rz,zz,psin0 real(wp_), intent(out) :: rf,zf integer, intent(out) :: info ! local variables real(wp_) :: tol real(wp_), dimension(n) :: xvec,fvec,f0 real(wp_), dimension(lwa) :: wa real(wp_), dimension(ldfjac,n) :: fjac xvec(1)=rz xvec(2)=zz f0(1)=psin0 f0(2)=0.0_wp_ tol = sqrt(comp_eps) call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then print'(a,i2,a,5f8.4)',' info subr points_tgo =',info, & ' R,z coord.',xvec,rz,zz,psin0 end if rf=xvec(1) zf=xvec(2) end subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) use const_and_precisions, only : wp_ implicit none ! arguments integer, intent(in) :: n,ldfjac,iflag real(wp_), dimension(n), intent(in) :: x,f0 real(wp_), dimension(n), intent(inout) :: fvec real(wp_), dimension(ldfjac,n), intent(inout) :: fjac ! internal variables real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz select case(iflag) case(1) call equinum_psi(x(1),x(2),psinv,dpsidr) fvec(1) = psinv-f0(1) fvec(2) = dpsidr/psia-f0(2) case(2) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) fjac(1,1) = dpsidr/psia fjac(1,2) = dpsidz/psia fjac(2,1) = ddpsidrr/psia fjac(2,2) = ddpsidrz/psia case default print*,'iflag undefined' end select end subroutine fcntgo subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n) use const_and_precisions, only : pi,zero,one implicit none ! arguments real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp integer, intent(in), optional :: n ! local variables integer, parameter :: nqdef=101 integer :: i real(wp_) :: dr,fq0,fq1,qq,res,rn real(wp_), dimension(:), allocatable :: rhotn,rhopn btaxis=bax rmaxis=rax zmaxis=zax btrcen=bax rcen=rax aminor=a zbinf=zmaxis-a zbsup=zmaxis+a q0=qax qa=q1 alq=qexp sgnbphi=sign(one,bax) if (present(n)) then nq=n else nq=nqdef end if if (allocated(psinr)) deallocate(psinr) allocate(psinr(nq),rhotn(nq),rhopn(nq)) dr=one/(nq-1) rhotn(1)=zero psinr(1)=zero res=zero fq0=zero do i=2,nq rn=(i-1)*dr qq=q0+(q1-q0)*rn**qexp fq1=rn/qq res=res+0.5_wp_*(fq1+fq0)*dr fq0=fq1 rhotn(i)=rn psinr(i)=res end do phitedge=btaxis*aminor**2 ! temporary psia=res*phitedge phitedge=pi*phitedge ! final psinr=psinr/res rhopn=sqrt(psinr) call set_rhospl(rhopn,rhotn) end subroutine set_equian subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) use const_and_precisions, only : wp_ implicit none ! arguments real(wp_), intent(in) :: rrm,zzm real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz ! local variables real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot,btaxqq ! simple model for equilibrium: large aspect ratio ! outside plasma: analytical continuation, not solution Maxwell eqs rpm=sqrt((rrm-rmaxis)**2+(zzm-zmaxis)**2) !!! rpm==rho_tor[m], rn=rho_tor_norm rn=rpm/aminor snt=0.0_wp_ cst=1.0_wp_ if (rpm > 0.0_wp_) then snt=(zzm-zmaxis)/rpm cst=(rrm-rmaxis)/rpm end if if (present(psinv)) then rhot=rn if(rn <= 1.0_wp_) then rhop=frhopol(rhot) psinv=rhop**2 else psinv=1.0_wp_+btaxis/(2.0_wp_*psia*qa)*(rpm**2-aminor**2) rhop=sqrt(psinv) end if end if if(rn <= 1.0_wp_) then qq=q0+(qa-q0)*rn**alq btaxqq=btaxis/qq dpsidrp=btaxqq*rpm dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) d2psidrp=btaxqq*(1.0_wp_-rn*dqq/qq) else btaxqq=btaxis/qa dpsidrp=btaxqq*rpm d2psidrp=btaxqq end if if(present(fpolv)) fpolv=btaxis*rmaxis if(present(dfpv)) dfpv=0.0_wp_ if(present(dpsidr)) dpsidr=dpsidrp*cst if(present(dpsidz)) dpsidz=dpsidrp*snt if(present(ddpsidrr)) ddpsidrr=btaxqq*snt**2+cst**2*d2psidrp if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-btaxqq) if(present(ddpsidzz)) ddpsidzz=btaxqq*cst**2+snt**2*d2psidrp end subroutine equian subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : wp_,ccj=>mu0inv use gray_params, only : iequil implicit none ! arguments real(wp_) :: rpsim,zpsim,ajphi ! local variables real(wp_) :: bzz,dbvcdc13,dbvcdc31 real(wp_) :: dpsidr,ddpsidrr,ddpsidzz if(iequil < 2) then call equian(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) else call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) end if bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) end subroutine tor_curr subroutine psi_raxis(psin,r1,r2) use const_and_precisions, only : wp_ use gray_params, only : iequil use dierckx, only : profil,sproota implicit none ! local constants integer, parameter :: mest=4 ! arguments real(wp_) :: psin,r1,r2 ! local variables integer :: iopt,ier,m real(wp_) :: zc,val real(wp_), dimension(mest) :: zeroc real(wp_), dimension(nsr) :: czc if (iequil < 2) then val=frhotor(sqrt(psin)) r1=rmaxis-val*aminor r2=rmaxis+val*aminor else iopt=1 zc=zmaxis call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier val=psin*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) end if end subroutine psi_raxis subroutine tor_curr_psi(psin,ajphi) use const_and_precisions, only : wp_ implicit none ! arguments real(wp_) :: psin,ajphi ! local variables real(wp_) :: r1,r2 call psi_raxis(psin,r1,r2) call tor_curr(r2,zmaxis,ajphi) end subroutine tor_curr_psi end module equilibrium