From 198feb7a1ff470132d4f31334beee668866209ad Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Thu, 5 Nov 2015 18:24:21 +0000 Subject: [PATCH] flux_average now in magsurf_data works anal/num , new subr fluxval, some routines moved into equilibrium, minor bugs fixed --- Makefile | 8 +- src/equilibrium.f90 | 310 ++++--- src/gray-externals.f | 1876 +++++++++++++++++++++--------------------- src/graycore.f90 | 16 +- src/magsurf_data.f90 | 403 ++++++++- 5 files changed, 1556 insertions(+), 1057 deletions(-) diff --git a/Makefile b/Makefile index 3d5f572..a375ee5 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ MAINOBJ=main.o OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ dierckx.o dispersion.o eccd.o eierf.o graycore.o gray-externals.o \ gray_params.o equilibrium.o magsurf_data.o math.o minpack.o numint.o \ - quadpack.o reflections.o simplespline.o utils.o + quadpack.o reflections.o simplespline.o utils.o # Alternative search paths vpath %.f90 src @@ -44,8 +44,10 @@ dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o eierf.o: const_and_precisions.o gray_params.o: const_and_precisions.o utils.o -equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o utils.o -magsurf_data.o: const_and_precisions.o +equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o \ + utils.o gray_params.o +magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \ + simplespline.o utils.o math.o: const_and_precisions.o minpack.o: const_and_precisions.o numint.o: const_and_precisions.o diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index c47a2db..3c66e0b 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -2,32 +2,35 @@ 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 + 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, & + 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 + +! === 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 @@ -146,15 +149,14 @@ contains 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 -! integer, intent(in) :: isgnbphi -! real(wp_), intent(in) :: factb -! real(wp_), intent(out) :: rvac,rax,zax real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q ! local variables integer :: u @@ -199,7 +201,6 @@ contains real(wp_), intent(inout) :: psia real(wp_), dimension(:), intent(inout) :: fpol,q integer, intent(in) :: cocosin, cocosout -! real(wp_), intent(out) :: isign,bsign integer, intent(out), optional :: ierr ! local variables real(wp_) :: isign,bsign @@ -209,7 +210,7 @@ contains call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos) call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout) -! check sign consistency +! check sign consistency isign=sign(one,psia) if (.not.psiincr) isign=-isign bsign=sign(one,fpol(size(fpol))) @@ -222,17 +223,14 @@ contains end if ! convert cocosin to cocosout - if (phiccw.neqv.phiccwout) then -! opposite direction of toroidal angle phi in cocosin and cocosout -! bsign=-bsign -! isign=-isign - fpol=-fpol - end if -! q has opposite sign for given sign of Bphi*Ip + +! 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 +! 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 +! convert Wb to Wb/rad or viceversa psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi) end subroutine change_cocos @@ -301,8 +299,8 @@ contains real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk integer :: ier,ixploc,info -! -! compute array sizes and prepare working space arrays + +! compute array sizes and prepare working space arrays nr=size(rv) nz=size(zv) nrz=nr*nz @@ -310,34 +308,34 @@ contains 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 !!! -! + +! 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 +! 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 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 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, & @@ -345,7 +343,7 @@ contains wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) end if deallocate(fvpsi) -! compute spline coefficients for psi partial derivatives +! 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 @@ -362,10 +360,10 @@ contains 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 + +! 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)) @@ -375,22 +373,22 @@ contains 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 +! set vacuum value used outside 0<=psin<=1 range fpolas=fpoli(1) sgnbphi=sign(one,fpolas) -! free temporary arrays +! free temporary arrays deallocate(wrk,iwrk,wf) -! -! re-normalize psi after spline computation -! -! start with un-corrected psi -! + +! 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 -! + +! 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)) @@ -408,9 +406,9 @@ contains rbmin=rv(2) rbmax=rv(nr-1) end if -! -! search for exact location of the magnetic axis -! + +! search for exact location of the magnetic axis + if(present(rax)) then rax0=rax else @@ -423,9 +421,9 @@ contains 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 -! + +! search for X-point if ixp not = 0 + if(present(ixp)) then ixploc=ixp else @@ -462,21 +460,21 @@ contains if (ixploc==0) then psinop=psinoptmp psiant=one-psinop -! find upper horizontal tangent point +! 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 +! 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) -! + +! 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 @@ -522,9 +520,9 @@ contains integer, dimension(liwrk) :: iwrk real(wp_), dimension(1) :: rrs,zzs,ffspl real(wp_), dimension(lwrk) :: wrk -! -! here lengths are measured in meters -! + +! here lengths are measured in meters + if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. & zpsim.le.zmxm .and. zpsim.ge.zmnm) then @@ -592,7 +590,7 @@ contains 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) @@ -608,6 +606,7 @@ contains 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 @@ -615,11 +614,16 @@ contains ! local variables real(wp_) :: psin,fpol - 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 + 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 @@ -629,11 +633,11 @@ contains use const_and_precisions, only : pi use simplespline, only : difcs implicit none -! arguments +! arguments real(wp_), dimension(:), intent(in) :: psinq,q real(wp_), intent(in) :: psia real(wp_), dimension(:), intent(out), optional :: rhotn -! local variables +! local variables real(wp_), dimension(size(q)) :: phit real(wp_) :: dx integer, parameter :: iopt=0 @@ -646,9 +650,8 @@ contains psinr=psinq call difcs(psinr,q,nq,iopt,cq,ier) -! -! Toroidal flux phi = 2*pi*Integral q dpsi -! + +! Toroidal flux phi = 2*pi*Integral q dpsi phit(1)=0.0_wp_ do k=1,nq-1 dx=psinr(k+1)-psinr(k) @@ -729,10 +732,10 @@ contains use utils, only : locate use simplespline, only : spli implicit none -! arguments +! arguments real(wp_), intent(in) :: rhot real(wp_) :: frhopol -! local variables +! local variables integer :: i real(wp_) :: dr @@ -746,10 +749,10 @@ contains use utils, only : locate use simplespline, only : spli implicit none -! arguments +! arguments real(wp_), dimension(:), intent(in) :: rhot real(wp_), dimension(size(rhot)) :: frhopolv -! local variables +! local variables integer :: i,i0,j real(wp_) :: dr @@ -767,10 +770,10 @@ contains use utils, only : locate use simplespline, only : spli implicit none -! arguments +! arguments real(wp_), intent(in) :: rhop real(wp_) :: frhotor -! local variables +! local variables integer :: i real(wp_) :: dr @@ -818,7 +821,7 @@ contains 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) @@ -826,7 +829,7 @@ contains fvec(2) = dpsidz/psia case(2) call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & - ddpsidrz=ddpsidrz) + ddpsidrz=ddpsidrz) fjac(1,1) = ddpsidrr/psia fjac(1,2) = ddpsidrz/psia fjac(2,1) = ddpsidrz/psia @@ -851,7 +854,6 @@ contains real(wp_), dimension(n) :: xvec,fvec,f0 real(wp_), dimension(lwa) :: wa real(wp_), dimension(ldfjac,n) :: fjac -! common/external functions/variables xvec(1)=rz xvec(2)=zz @@ -860,8 +862,8 @@ contains 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,2f8.4)',' info subr points_tgo =',info, & - ' R,z coord.',xvec + 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) @@ -885,7 +887,7 @@ contains fvec(2) = dpsidr/psia-f0(2) case(2) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & - ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) + ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) fjac(1,1) = dpsidr/psia fjac(1,2) = dpsidz/psia fjac(2,1) = ddpsidrr/psia @@ -927,12 +929,13 @@ contains 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,n + do i=2,nq rn=(i-1)*dr qq=q0+(q1-q0)*rn**qexp fq1=rn/qq @@ -941,11 +944,13 @@ contains 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 @@ -958,11 +963,10 @@ contains 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 -! real(wp_) :: frhopol + 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 +! 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 @@ -987,12 +991,14 @@ contains if(rn <= 1.0_wp_) then qq=q0+(qa-q0)*rn**alq - dpsidrp=btaxis*aminor*rn/qq + btaxqq=btaxis/qq + dpsidrp=btaxqq*rpm dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) - d2psidrp=btaxis/qq*(1.0_wp_-rn*dqq/qq) + d2psidrp=btaxqq*(1.0_wp_-rn*dqq/qq) else - dpsidrp=btaxis*rpm/qa - d2psidrp=btaxis/qa + btaxqq=btaxis/qa + dpsidrp=btaxqq*rpm + d2psidrp=btaxqq end if if(present(fpolv)) fpolv=btaxis*rmaxis @@ -1000,10 +1006,80 @@ contains if(present(dpsidr)) dpsidr=dpsidrp*cst if(present(dpsidz)) dpsidz=dpsidrp*snt - if(present(ddpsidrr)) ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp - if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm) - if(present(ddpsidzz)) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp - + 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 diff --git a/src/gray-externals.f b/src/gray-externals.f index f2a304b..e95dc42 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -588,7 +588,7 @@ c c subroutine print_prof use const_and_precisions, only : wp_ - use equilibrium, only : psinr,nq,fq,frhotor + use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi use coreprofiles, only : density, temp implicit none c local constants @@ -987,780 +987,780 @@ c c c c - subroutine flux_average - use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use magsurf_data - use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, - . equinum_fpol,bfield,fq,frhotor,zbsup,zbinf - use simplespline, only : difcs - use dierckx, only : regrid,coeff_parder - implicit none -c local constants - integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, - . njest=nnintp+ksp+1,nlest=nlam+ksp+1, - . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, - . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam -c local variables - integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr - integer, dimension(kwrk) :: iwrk - real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, - . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp, - . area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, - . dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, - . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, - . bphi,brr,bzz,riav,fp,rup,rlw,zup,zlw - real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl - real(wp_), dimension(2*ncnt) :: dlpv - real(wp_), dimension(2*ncnt+1) :: bv,bpv - real(wp_), dimension(nlam) :: alam,weights - real(wp_), dimension(nnintp,nlam) :: fhlam - real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam - real(wp_), dimension(lwrk) :: wrk - real(wp_), dimension(:), allocatable :: rctemp,zctemp -c common/external functions/variables - real(wp_) :: fpolv -c - npsi=nnintp - ninpr=(npsi-1)/10 - npoints = 2*ncnt+1 -c - call alloc_surfvec(ierr) - if(allocated(tjp)) deallocate(tjp) - if(allocated(tlm)) deallocate(tlm) - if(allocated(ch)) deallocate(ch) - allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), - . rctemp(npoints),zctemp(npoints),stat=ierr) - if (ierr.ne.0) return - -c computation of flux surface averaged quantities - - write(71,*)' #i psin R z' - - dlam=1.0_wp_/dble(nlam-1) - do l=1,nlam-1 - alam(l)=dble(l-1)*dlam - fhlam(1,l)=sqrt(1.0_wp_-alam(l)) - ffhlam(l)=fhlam(1,l) - dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) - weights(l)=1.0_wp_ - end do - weights(1)=0.5_wp_ - weights(nlam)=0.5_wp_ - alam(nlam)=1.0_wp_ - fhlam(1,nlam)=0.0_wp_ - ffhlam(nlam)=0.0_wp_ - dffhlam(nlam)=-99999.0_wp_ - - jp=1 - anorm=2.0_wp_*pi*rmaxis/abs(btaxis) - b2av=btaxis**2 - dvdpsi=2.0_wp_*pi*anorm - dadpsi=2.0_wp_*pi/abs(btaxis) - ratio_cdator=abs(btaxis/btrcen) - ratio_cdbtor=1.0_wp_ - ratio_pltor=1.0_wp_ - qq=1.0_wp_ - fc=1.0_wp_ - - psicon(1)=0.0_wp_ - rcon(1,:)=rmaxis - zcon(1,:)=zmaxis - pstab(1)=0.0_wp_ - rhot_eq(1)=0.0_wp_ - rpstab(1)=0.0_wp_ - rhotqv(1)=0.0_wp_ - vcurrp(1)=0.0_wp_ - vajphiav(1)=0.0_wp_ - bmxpsi(1)=abs(btaxis) - bmnpsi(1)=abs(btaxis) - bav(1)=abs(btaxis) - rbav(1)=1.0_wp_ - rri(1)=rmaxis - varea(1)=0.0_wp_ - vvol(1)=0.0_wp_ - vratjpl(1)=ratio_pltor - vratja(1)=ratio_cdator - vratjb(1)=ratio_cdbtor - ffc(1)=fc - rhot2q=0.0_wp_ - qqv(1)=1.0_wp_ - - dadrhotv(1)=0.0_wp_ - dvdrhotv(1)=0.0_wp_ - - rup=rmaxis - rlw=rmaxis - zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ - zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ - - do jp=2,npsi - height=dble(jp-1)/dble(npsi-1) - psicon(jp)=height - if(jp.eq.npsi) height=0.9999_wp_ - ipr=0 - jpr=mod(jp,ninpr) - if(jpr.eq.1) ipr=1 - height2=height*height - call contours_psi(height2,rup,zup,rlw,zlw,rctemp,zctemp,ipr) - rcon(jp,:) = rctemp - zcon(jp,:) = zctemp -c - r2iav=0.0_wp_ - anorm=0.0_wp_ - dadpsi=0.0_wp_ - currp=0.0_wp_ - b2av=0.0_wp_ - area=0.0_wp_ - volume=0.0_wp_ - ajphiav=0.0_wp_ - bbav=0.0_wp_ - bmmx=-1.0e+30_wp_ - bmmn=1.0e+30_wp_ - - call tor_curr(rctemp(1),zctemp(1),ajphi0) - call bfield(rctemp(1),zctemp(1),br=brr,bz=bzz) - call equinum_fpol(height2,fpolv) - bphi=fpolv/rctemp(1) - btot0=sqrt(bphi**2+brr**2+bzz**2) - bpoloid0=sqrt(brr**2+bzz**2) - bv(1)=btot0 - bpv(1)=bpoloid0 - rpsim0=rctemp(1) - - do inc=1,npoints-1 - inc1=inc+1 - dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) - dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) - dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ - . (zctemp(inc1)-zctemp(inc))**2) - drc=(rctemp(inc1)-rctemp(inc)) - -c compute length, area and volume defined by psi=height^2 - - ph=0.5_wp_*(dla+dlb+dlp) - area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) - area=area+sqrt(area2) - rzp=rctemp(inc1)*zctemp(inc1) - rz=rctemp(inc)*zctemp(inc) - volume=pi*(rzp+rz)*drc+volume - -c compute line integral on the contour psi=height^2 - - rpsim=rctemp(inc1) - zpsim=zctemp(inc1) - call bfield(rpsim,zpsim,br=brr,bz=bzz) - call tor_curr(rpsim,zpsim,ajphi) -! call equinum_fpol(height2,fpolv) - bphi=fpolv/rpsim - btot=sqrt(bphi**2+brr**2+bzz**2) - bpoloid=sqrt(brr**2+bzz**2) - dlpv(inc)=dlp - bv(inc1)=btot - bpv(inc1)=bpoloid - - dlph=0.5_wp_*dlp - anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) - dadpsi=dadpsi+dlph* - . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) - currp=currp+dlph*(bpoloid+bpoloid0) - b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) - bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) - r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) - . +1.0_wp_/(bpoloid0*rpsim0**2)) - ajphiav=ajphiav+dlph* - . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) - - ajphi0=ajphi - rpsim0=rpsim - bpoloid0=bpoloid - btot0=btot - -c computation maximum/minimum B values on given flux surface - - if(btot.le.bmmn) bmmn=btot - if(btot.ge.bmmx) bmmx=btot - end do - -c bav= [T] , b2av= [T^2] , rbav=/b_min -c anorm = int d l_p/B_p = dV/dpsi/(2pi) -c r2iav=<1/R^2> [m^-2] , -c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), -c rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] -c currp = plasma current within psi=const - - bbav=bbav/anorm - r2iav=r2iav/anorm - dvdpsi=2.0_wp_*pi*anorm - riav=dadpsi/anorm - b2av=b2av/anorm - vcurrp(jp)=ccj*currp - vajphiav(jp)=ajphiav/dadpsi - -c area == varea, volume == vvol -c flux surface minor radius == (area/pi)^1/2 -c ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 -c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / -c ratio_pltor = Jcd_||/J_phi Jcd_|| = - - pstab(jp)=height2 - rpstab(jp)=height - vvol(jp)=abs(volume) - varea(jp)=area - bav(jp)=bbav - rbav(jp)=bbav/bmmn - bmxpsi(jp)=bmmx - bmnpsi(jp)=bmmn - rri(jp)=bav(jp)/abs(fpolv*r2iav) - ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) - ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) - ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) - vratjpl(jp)=ratio_pltor - vratja(jp)=ratio_cdator - vratjb(jp)=ratio_cdbtor - qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) - qqv(jp)=qq - - dadrhotv(jp)=phitedge*frhotor(height)/fq(height2) - . *dadpsi/pi - dvdrhotv(jp)=phitedge*frhotor(height)/fq(height2) - . *dvdpsi/pi - -c computation of rhot from calculated q profile - - rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) - . /dble(npsi-1) - rhotqv(jp)=sqrt(rhot2q) - -c computation of fraction of circulating/trapped fraction fc, ft -c and of function H(lambda,rhop) -c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ - - fc=0.0_wp_ - shlam=0.0_wp_ - do l=nlam,1,-1 - lam=alam(l) - srl=0.0_wp_ - rl2=1.0_wp_-lam*bv(1)/bmmx - rl0=0.0_wp_ - if(rl2.gt.0) rl0=sqrt(rl2) - do inc=1,npoints-1 - rl2=1.0_wp_-lam*bv(inc+1)/bmmx - rl=0.0_wp_ - if(rl2.gt.0) rl=sqrt(rl2) - srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) - rl0=rl - end do - srl=srl/anorm - dhlam=0.5_wp_/srl - fc=fc+lam/srl*weights(l) - if(l.eq.nlam) then - fhlam(jp,l)=0.0_wp_ - ffhlam(nlam*(jp-1)+l)=0.0_wp_ - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - else - shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam - fhlam(jp,l)=shlam - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - end if - end do - fc=0.75_wp_*b2av/bmmx**2*fc*dlam - ffc(jp)=fc - - ccfh=bmmn/bmmx/fc - do l=1,nlam - ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) - dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) - end do - end do - - write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// - .' Area Vol |I_pl| qq fc ratioJa ratioJb' - - qqv(1)=qqv(2) - vajphiav(1)=vajphiav(2) - do jp=1,npsi - rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) - if(jp.eq.npsi) then - rhotqv(jp)=1.0_wp_ - rpstab(jp)=1.0_wp_ - pstab(jp)=1.0_wp_ - end if - rhot_eq(jp)=frhotor(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), - . vratja(jp),vratjb(jp) - end do - -c rarea=sqrt(varea(npsi)/pi) - -c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi -c used for computations of dP/dV and J_cd -c spline coefficients of rhot - - iopt=0 - call difcs(rpstab,vvol,npsi,iopt,cvol,ier) - iopt=0 - call difcs(rpstab,rbav,npsi,iopt,crbav,ier) - iopt=0 - call difcs(rpstab,rri,npsi,iopt,crri,ier) - iopt=0 - call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) - iopt=0 - call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) - iopt=0 - call difcs(rpstab,vratja,npsi,iopt,cratja,ier) - iopt=0 - call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) - iopt=0 - call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) - iopt=0 - call difcs(rpstab,varea,npsi,iopt,carea,ier) - iopt=0 - call difcs(rpstab,ffc,npsi,iopt,cfc,ier) - iopt=0 - call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) - iopt=0 - call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) - iopt=0 - call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) - -c spline interpolation of H(lambda,rhop) and dH/dlambda - - iopt=0 - s=0.0_wp_ - call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, - . zero,one,zero,one,ksp,ksp,s, - . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) - njpt=njp - nlmt=nlm - - return - 99 format(20(1x,e12.5)) - end -c -c -c - function fdadrhot(rpsi) - use const_and_precisions, only : wp_ - use magsurf_data, only : rpstab,cdadrhot,npsi - use simplespline, only :spli - implicit none -c arguments - real(wp_) :: rpsi,fdadrhot -c local variables - integer :: ip - real(wp_) :: dps -c - ip=int((npsi-1)*rpsi+1) -c if(ip.eq.0) ip=1 -c if(ip.eq.npsi) ip=npsi-1 - ip=min(max(1,ip),npsi-1) - dps=rpsi-rpstab(ip) - fdadrhot=spli(cdadrhot,npsi,ip,dps) - return - end -c -c -c - function fdvdrhot(rpsi) - use const_and_precisions, only : wp_ - use magsurf_data, only : rpstab,cdvdrhot,npsi - use simplespline, only :spli - implicit none -c arguments - real(wp_) :: rpsi,fdvdrhot -c local variables - integer :: ip - real(wp_) :: dps -c - ip=int((npsi-1)*rpsi+1) - ip=min(max(1,ip),npsi-1) - dps=rpsi-rpstab(ip) - fdvdrhot=spli(cdvdrhot,npsi,ip,dps) -c - return - end +C subroutine flux_average +C use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv +C use magsurf_data +C use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, +C . equinum_fpol,bfield,fq,frhotor,zbsup,zbinf +C use simplespline, only : difcs +C use dierckx, only : regrid,coeff_parder +C implicit none +Cc local constants +C integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, +C . njest=nnintp+ksp+1,nlest=nlam+ksp+1, +C . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, +C . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam +Cc local variables +C integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr +C integer, dimension(kwrk) :: iwrk +C real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, +C . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp, +C . area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, +C . dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, +C . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, +C . bphi,brr,bzz,riav,fp,rup,rlw,zup,zlw +C real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl +C real(wp_), dimension(2*ncnt) :: dlpv +C real(wp_), dimension(2*ncnt+1) :: bv,bpv +C real(wp_), dimension(nlam) :: alam,weights +C real(wp_), dimension(nnintp,nlam) :: fhlam +C real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam +C real(wp_), dimension(lwrk) :: wrk +C real(wp_), dimension(:), allocatable :: rctemp,zctemp +Cc common/external functions/variables +C real(wp_) :: fpolv +Cc +C npsi=nnintp +C ninpr=(npsi-1)/10 +C npoints = 2*ncnt+1 +Cc +C call alloc_surfvec(ierr) +C if(allocated(tjp)) deallocate(tjp) +C if(allocated(tlm)) deallocate(tlm) +C if(allocated(ch)) deallocate(ch) +C allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), +C . rctemp(npoints),zctemp(npoints),stat=ierr) +C if (ierr.ne.0) return +C +Cc computation of flux surface averaged quantities +C +C write(71,*)' #i psin R z' +C +C dlam=1.0_wp_/dble(nlam-1) +C do l=1,nlam-1 +C alam(l)=dble(l-1)*dlam +C fhlam(1,l)=sqrt(1.0_wp_-alam(l)) +C ffhlam(l)=fhlam(1,l) +C dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) +C weights(l)=1.0_wp_ +C end do +C weights(1)=0.5_wp_ +C weights(nlam)=0.5_wp_ +C alam(nlam)=1.0_wp_ +C fhlam(1,nlam)=0.0_wp_ +C ffhlam(nlam)=0.0_wp_ +C dffhlam(nlam)=-99999.0_wp_ +C +C jp=1 +C anorm=2.0_wp_*pi*rmaxis/abs(btaxis) +C b2av=btaxis**2 +C dvdpsi=2.0_wp_*pi*anorm +C dadpsi=2.0_wp_*pi/abs(btaxis) +C ratio_cdator=abs(btaxis/btrcen) +C ratio_cdbtor=1.0_wp_ +C ratio_pltor=1.0_wp_ +C qq=1.0_wp_ +C fc=1.0_wp_ +C +C psicon(1)=0.0_wp_ +C rcon(1,:)=rmaxis +C zcon(1,:)=zmaxis +C pstab(1)=0.0_wp_ +C rhot_eq(1)=0.0_wp_ +C rpstab(1)=0.0_wp_ +C rhotqv(1)=0.0_wp_ +C vcurrp(1)=0.0_wp_ +C vajphiav(1)=0.0_wp_ +C bmxpsi(1)=abs(btaxis) +C bmnpsi(1)=abs(btaxis) +C bav(1)=abs(btaxis) +C rbav(1)=1.0_wp_ +C rri(1)=rmaxis +C varea(1)=0.0_wp_ +C vvol(1)=0.0_wp_ +C vratjpl(1)=ratio_pltor +C vratja(1)=ratio_cdator +C vratjb(1)=ratio_cdbtor +C ffc(1)=fc +C rhot2q=0.0_wp_ +C qqv(1)=1.0_wp_ +C +C dadrhotv(1)=0.0_wp_ +C dvdrhotv(1)=0.0_wp_ +C +C rup=rmaxis +C rlw=rmaxis +C zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ +C zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ +C +C do jp=2,npsi +C height=dble(jp-1)/dble(npsi-1) +C psicon(jp)=height +C if(jp.eq.npsi) height=0.9999_wp_ +C ipr=0 +C jpr=mod(jp,ninpr) +C if(jpr.eq.1) ipr=1 +C height2=height*height +C call contours_psi(height2,rup,zup,rlw,zlw,rctemp,zctemp,ipr) +C rcon(jp,:) = rctemp +C zcon(jp,:) = zctemp +Cc +C r2iav=0.0_wp_ +C anorm=0.0_wp_ +C dadpsi=0.0_wp_ +C currp=0.0_wp_ +C b2av=0.0_wp_ +C area=0.0_wp_ +C volume=0.0_wp_ +C ajphiav=0.0_wp_ +C bbav=0.0_wp_ +C bmmx=-1.0e+30_wp_ +C bmmn=1.0e+30_wp_ +C +C call tor_curr(rctemp(1),zctemp(1),ajphi0) +C call bfield(rctemp(1),zctemp(1),br=brr,bz=bzz) +C call equinum_fpol(height2,fpolv) +C bphi=fpolv/rctemp(1) +C btot0=sqrt(bphi**2+brr**2+bzz**2) +C bpoloid0=sqrt(brr**2+bzz**2) +C bv(1)=btot0 +C bpv(1)=bpoloid0 +C rpsim0=rctemp(1) +C +C do inc=1,npoints-1 +C inc1=inc+1 +C dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) +C dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) +C dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ +C . (zctemp(inc1)-zctemp(inc))**2) +C drc=(rctemp(inc1)-rctemp(inc)) +C +Cc compute length, area and volume defined by psi=height^2 +C +C ph=0.5_wp_*(dla+dlb+dlp) +C area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) +C area=area+sqrt(area2) +C rzp=rctemp(inc1)*zctemp(inc1) +C rz=rctemp(inc)*zctemp(inc) +C volume=pi*(rzp+rz)*drc+volume +C +Cc compute line integral on the contour psi=height^2 +C +C rpsim=rctemp(inc1) +C zpsim=zctemp(inc1) +C call bfield(rpsim,zpsim,br=brr,bz=bzz) +C call tor_curr(rpsim,zpsim,ajphi) +C! call equinum_fpol(height2,fpolv) +C bphi=fpolv/rpsim +C btot=sqrt(bphi**2+brr**2+bzz**2) +C bpoloid=sqrt(brr**2+bzz**2) +C dlpv(inc)=dlp +C bv(inc1)=btot +C bpv(inc1)=bpoloid +C +C dlph=0.5_wp_*dlp +C anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) +C dadpsi=dadpsi+dlph* +C . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) +C currp=currp+dlph*(bpoloid+bpoloid0) +C b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) +C bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) +C r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) +C . +1.0_wp_/(bpoloid0*rpsim0**2)) +C ajphiav=ajphiav+dlph* +C . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) +C +C ajphi0=ajphi +C rpsim0=rpsim +C bpoloid0=bpoloid +C btot0=btot +C +Cc computation maximum/minimum B values on given flux surface +C +C if(btot.le.bmmn) bmmn=btot +C if(btot.ge.bmmx) bmmx=btot +C end do +C +Cc bav= [T] , b2av= [T^2] , rbav=/b_min +Cc anorm = int d l_p/B_p = dV/dpsi/(2pi) +Cc r2iav=<1/R^2> [m^-2] , +Cc riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), +Cc rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] +Cc currp = plasma current within psi=const +C +C bbav=bbav/anorm +C r2iav=r2iav/anorm +C dvdpsi=2.0_wp_*pi*anorm +C riav=dadpsi/anorm +C b2av=b2av/anorm +C vcurrp(jp)=ccj*currp +C vajphiav(jp)=ajphiav/dadpsi +C +Cc area == varea, volume == vvol +Cc flux surface minor radius == (area/pi)^1/2 +Cc ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 +Cc ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / +Cc ratio_pltor = Jcd_||/J_phi Jcd_|| = +C +C pstab(jp)=height2 +C rpstab(jp)=height +C vvol(jp)=abs(volume) +C varea(jp)=area +C bav(jp)=bbav +C rbav(jp)=bbav/bmmn +C bmxpsi(jp)=bmmx +C bmnpsi(jp)=bmmn +C rri(jp)=bav(jp)/abs(fpolv*r2iav) +C ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) +C ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) +C ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) +C vratjpl(jp)=ratio_pltor +C vratja(jp)=ratio_cdator +C vratjb(jp)=ratio_cdbtor +C qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) +C qqv(jp)=qq +C +C dadrhotv(jp)=phitedge*frhotor(height)/fq(height2) +C . *dadpsi/pi +C dvdrhotv(jp)=phitedge*frhotor(height)/fq(height2) +C . *dvdpsi/pi +C +Cc computation of rhot from calculated q profile +C +C rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) +C . /dble(npsi-1) +C rhotqv(jp)=sqrt(rhot2q) +C +Cc computation of fraction of circulating/trapped fraction fc, ft +Cc and of function H(lambda,rhop) +Cc ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ +C +C fc=0.0_wp_ +C shlam=0.0_wp_ +C do l=nlam,1,-1 +C lam=alam(l) +C srl=0.0_wp_ +C rl2=1.0_wp_-lam*bv(1)/bmmx +C rl0=0.0_wp_ +C if(rl2.gt.0) rl0=sqrt(rl2) +C do inc=1,npoints-1 +C rl2=1.0_wp_-lam*bv(inc+1)/bmmx +C rl=0.0_wp_ +C if(rl2.gt.0) rl=sqrt(rl2) +C srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) +C rl0=rl +C end do +C srl=srl/anorm +C dhlam=0.5_wp_/srl +C fc=fc+lam/srl*weights(l) +C if(l.eq.nlam) then +C fhlam(jp,l)=0.0_wp_ +C ffhlam(nlam*(jp-1)+l)=0.0_wp_ +C dffhlam(nlam*(jp-1)+l)=-dhlam +C dhlam0=dhlam +C else +C shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam +C fhlam(jp,l)=shlam +C dffhlam(nlam*(jp-1)+l)=-dhlam +C dhlam0=dhlam +C end if +C end do +C fc=0.75_wp_*b2av/bmmx**2*fc*dlam +C ffc(jp)=fc +C +C ccfh=bmmn/bmmx/fc +C do l=1,nlam +C ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) +C dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) +C end do +C end do +C +C write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// +C .' Area Vol |I_pl| qq fc ratioJa ratioJb' +C +C qqv(1)=qqv(2) +C vajphiav(1)=vajphiav(2) +C do jp=1,npsi +C rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) +C if(jp.eq.npsi) then +C rhotqv(jp)=1.0_wp_ +C rpstab(jp)=1.0_wp_ +C pstab(jp)=1.0_wp_ +C end if +C rhot_eq(jp)=frhotor(sqrt(pstab(jp))) +C write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), +C . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), +C . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), +C . vratja(jp),vratjb(jp) +C end do +C +Cc rarea=sqrt(varea(npsi)/pi) +C +Cc spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi +Cc used for computations of dP/dV and J_cd +Cc spline coefficients of rhot +C +C iopt=0 +C call difcs(rpstab,vvol,npsi,iopt,cvol,ier) +C iopt=0 +C call difcs(rpstab,rbav,npsi,iopt,crbav,ier) +C iopt=0 +C call difcs(rpstab,rri,npsi,iopt,crri,ier) +C iopt=0 +C call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) +C iopt=0 +C call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) +C iopt=0 +C call difcs(rpstab,vratja,npsi,iopt,cratja,ier) +C iopt=0 +C call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) +C iopt=0 +C call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) +C iopt=0 +C call difcs(rpstab,varea,npsi,iopt,carea,ier) +C iopt=0 +C call difcs(rpstab,ffc,npsi,iopt,cfc,ier) +C iopt=0 +C call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) +C iopt=0 +C call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) +C iopt=0 +C call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) +C +Cc spline interpolation of H(lambda,rhop) and dH/dlambda +C +C iopt=0 +C s=0.0_wp_ +C call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, +C . zero,one,zero,one,ksp,ksp,s, +C . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) +C njpt=njp +C nlmt=nlm +C +C return +C 99 format(20(1x,e12.5)) +C end +Cc +Cc +Cc +C function fdadrhot(rpsi) +C use const_and_precisions, only : wp_ +C use magsurf_data, only : rpstab,cdadrhot,npsi +C use simplespline, only :spli +C implicit none +Cc arguments +C real(wp_) :: rpsi,fdadrhot +Cc local variables +C integer :: ip +C real(wp_) :: dps +Cc +C ip=int((npsi-1)*rpsi+1) +Cc if(ip.eq.0) ip=1 +Cc if(ip.eq.npsi) ip=npsi-1 +C ip=min(max(1,ip),npsi-1) +C dps=rpsi-rpstab(ip) +C fdadrhot=spli(cdadrhot,npsi,ip,dps) +C return +C end +Cc +Cc +Cc +C function fdvdrhot(rpsi) +C use const_and_precisions, only : wp_ +C use magsurf_data, only : rpstab,cdvdrhot,npsi +C use simplespline, only :spli +C implicit none +Cc arguments +C real(wp_) :: rpsi,fdvdrhot +Cc local variables +C integer :: ip +C real(wp_) :: dps +Cc +C ip=int((npsi-1)*rpsi+1) +C ip=min(max(1,ip),npsi-1) +C dps=rpsi-rpstab(ip) +C fdvdrhot=spli(cdvdrhot,npsi,ip,dps) +Cc +C return +C end c c c - subroutine flux_average_an - use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use magsurf_data - use equilibrium, only : btrcen,frhotor,btaxis,rmaxis,zmaxis, - . fq,aminor,equian - use simplespline, only : difcs - use dierckx, only : regrid,coeff_parder - implicit none -c local constants - integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, - . njest=nnintp+ksp+1,nlest=nlam+ksp+1, - . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, - . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam -c local variables - integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr - integer, dimension(kwrk) :: iwrk - real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, - . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height2,r2iav,currp,area, - . volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb,dlp, - . drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0,shlam, - . srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi,bphi,brr,bzz,riav,fp, - . dadr,dbvcdc13,dbvcdc31,dvdr,height,phitedge,qqan,rn - real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl - real(wp_), dimension(2*ncnt) :: dlpv - real(wp_), dimension(2*ncnt+1) :: bv,bpv - real(wp_), dimension(nlam) :: alam,weights - real(wp_), dimension(nnintp,nlam) :: fhlam - real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam - real(wp_), dimension(lwrk) :: wrk - real(wp_), dimension(:), allocatable :: rctemp,zctemp - real(wp_) :: psinv !unused. can be removed when equian is modified - !to accept optional arguments -c common/external functions/variables - real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv -c - common/fpol/fpolv,dfpv - common/derip1/dpsidr,dpsidz - common/derip2/ddpsidrr,ddpsidzz,ddpsidrz -c - npsi=nnintp - ninpr=(npsi-1)/10 - npoints = 2*ncnt+1 -c - call alloc_surfvec(ierr) - if(allocated(tjp)) deallocate(tjp) - if(allocated(tlm)) deallocate(tlm) - if(allocated(ch)) deallocate(ch) - allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), - . rctemp(npoints),zctemp(npoints),stat=ierr) - if (ierr.ne.0) return -c -c computation of flux surface averaged quantities - - phitedge=pi*aminor**2*btaxis - - write(71,*)' #i psin R z' - - dlam=1.0_wp_/dble(nlam-1) - do l=1,nlam-1 - alam(l)=dble(l-1)*dlam - fhlam(1,l)=sqrt(1.0_wp_-alam(l)) - ffhlam(l)=fhlam(1,l) - dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) - weights(l)=1.0_wp_ - end do - weights(1)=0.5_wp_ - weights(nlam)=0.5_wp_ - alam(nlam)=1.0_wp_ - fhlam(1,nlam)=0.0_wp_ - ffhlam(nlam)=0.0_wp_ - dffhlam(nlam)=-99999.0_wp_ - - jp=1 - anorm=2.0_wp_*pi*rmaxis/abs(btaxis) - b2av=btaxis**2 - dvdpsi=2.0_wp_*pi*anorm - dadpsi=2.0_wp_*pi/abs(btaxis) - ratio_cdator=abs(btaxis/btrcen) - ratio_cdbtor=1.0_wp_ - ratio_pltor=1.0_wp_ - qq=1.0_wp_ - fc=1.0_wp_ - - psicon(1)=0.0_wp_ - rcon(1,:)=rmaxis - zcon(1,:)=zmaxis - pstab(1)=0.0_wp_ - rhot_eq(1)=0.0_wp_ - rpstab(1)=0.0_wp_ - rhotqv(1)=0.0_wp_ - vcurrp(1)=0.0_wp_ - vajphiav(1)=0.0_wp_ - bmxpsi(1)=abs(btaxis) - bmnpsi(1)=abs(btaxis) - bav(1)=abs(btaxis) - rbav(1)=1.0_wp_ - rri(1)=rmaxis - varea(1)=0.0_wp_ - vvol(1)=0.0_wp_ - vratjpl(1)=ratio_pltor - vratja(1)=ratio_cdator - vratjb(1)=ratio_cdbtor - ffc(1)=fc - rhot2q=0.0_wp_ - dadrhotv(1)=0.0_wp_ - dvdrhotv(1)=0.0_wp_ - qqv(1)=1.0_wp_ - -c rup=rmaxis -c rlw=rmaxis -c zup=(zbmax+zmaxis)/2.0_wp_ -c zlw=(zmaxis+zbmin)/2.0_wp_ - - do jp=2,npsi - height=dble(jp-1)/dble(npsi-1) - psicon(jp)=height - if(jp.eq.npsi) height=0.9999_wp_ - height2=height*height - ipr=0 - jpr=mod(jp,ninpr) - if(jpr.eq.1) ipr=1 - call contours_psi_an(height2,rctemp,zctemp,ipr) - rcon(jp,:) = rctemp - zcon(jp,:) = zctemp -c - r2iav=0.0_wp_ - anorm=0.0_wp_ - dadpsi=0.0_wp_ - currp=0.0_wp_ - b2av=0.0_wp_ - area=0.0_wp_ - volume=0.0_wp_ - ajphiav=0.0_wp_ - bbav=0.0_wp_ - bmmx=-1.0e+30_wp_ - bmmn=1.0e+30_wp_ - - call equian(rctemp(1),zctemp(1),psinv) - dbvcdc13=-ddpsidzz/rctemp(1) - dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) - ajphi=ccj*(dbvcdc13-dbvcdc31) - brr=-dpsidz/rctemp(1) - bzz= dpsidr/rctemp(1) - bphi=fpolv/rctemp(1) - btot0=sqrt(bphi**2+brr**2+bzz**2) - bpoloid0=sqrt(brr**2+bzz**2) - bv(1)=btot0 - bpv(1)=bpoloid0 - rpsim0=rctemp(1) - - do inc=1,npoints-1 - inc1=inc+1 - dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) - dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) - dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ - . (zctemp(inc1)-zctemp(inc))**2) - drc=(rctemp(inc1)-rctemp(inc)) - -c compute length, area and volume defined by psi=height^2 - - ph=0.5_wp_*(dla+dlb+dlp) - area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) - area=area+sqrt(area2) - rzp=rctemp(inc1)*zctemp(inc1) - rz=rctemp(inc)*zctemp(inc) - volume=pi*(rzp+rz)*drc+volume - -c compute line integral on the contour psi=height^2 - - rpsim=rctemp(inc1) - zpsim=zctemp(inc1) - call equian(rpsim,zpsim,psinv) - brr=-dpsidz/rpsim - bzz= dpsidr/rpsim - dbvcdc13=-ddpsidzz/rpsim - dbvcdc31= ddpsidrr/rpsim-bzz/rpsim - ajphi=ccj*(dbvcdc13-dbvcdc31) - bphi=fpolv/rpsim - btot=sqrt(bphi**2+brr**2+bzz**2) - bpoloid=sqrt(brr**2+bzz**2) - dlpv(inc)=dlp - bv(inc1)=btot - bpv(inc1)=bpoloid - - dlph=0.5_wp_*dlp - anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) - dadpsi=dadpsi+dlph* - . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) - currp=currp+dlph*(bpoloid+bpoloid0) - b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) - bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) - r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) - . +1.0_wp_/(bpoloid0*rpsim0**2)) - ajphiav=ajphiav+dlph* - . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) - - ajphi0=ajphi - rpsim0=rpsim - bpoloid0=bpoloid - btot0=btot -c computation maximum/minimum B values on given flux surface - - if(btot.le.bmmn) bmmn=btot - if(btot.ge.bmmx) bmmx=btot - end do - -c bav= [T] , b2av= [T^2] , rbav=/b_min -c anorm = int d l_p/B_p = dV/dpsi/(2pi) -c r2iav=<1/R^2> [m^-2] , -c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), -c rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] -c currp = plasma current within psi=const - - bbav=bbav/anorm - r2iav=r2iav/anorm - dvdpsi=2.0_wp_*pi*anorm - riav=dadpsi/anorm - b2av=b2av/anorm - vcurrp(jp)=ccj*currp - vajphiav(jp)=ajphiav/dadpsi - -c area == varea, volume == vvol -c flux surface minor radius == (area/pi)^1/2 -c ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 -c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / -c ratio_pltor = Jcd_||/J_phi Jcd_|| = - - pstab(jp)=height2 - rpstab(jp)=height - vvol(jp)=abs(volume) - varea(jp)=area - bav(jp)=bbav - rbav(jp)=bbav/bmmn - bmxpsi(jp)=bmmx - bmnpsi(jp)=bmmn - rri(jp)=bav(jp)/abs(fpolv*r2iav) - ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) - ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) - ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) - vratjpl(jp)=ratio_pltor - vratja(jp)=ratio_cdator - vratjb(jp)=ratio_cdbtor - qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) - qqv(jp)=qq - - rn=frhotor(sqrt(pstab(jp))) - qqan=fq(pstab(jp)) - - dadr=2.0_wp_*pi*rn*aminor**2 - dvdr=4.0_wp_*pi**2*rn*rmaxis*aminor**2 - -c dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi -c dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi - dadrhotv(jp)=phitedge*rn*dadpsi/pi/qqan - dvdrhotv(jp)=phitedge*rn*dvdpsi/pi/qqan - -c computation of rhot from calculated q profile - rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) - . /dble(npsi-1) -c print*,jp,rhot2q,qqv(jp),rpstab(jp),qqv(jp-1),rpstab(jp-1) - rhotqv(jp)=sqrt(rhot2q) -c rhotqv(jp)=rn -c - write(57,99) rpstab(jp),rn,rhotqv(jp),qqv(jp),qqan,r2iav, - . dadr,dadrhotv(jp),dadpsi,dvdpsi,fpolv - -c computation of fraction of circulating/trapped fraction fc, ft -c and of function H(lambda,rhop) -c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ - - fc=0.0_wp_ - shlam=0.0_wp_ - do l=nlam,1,-1 - lam=alam(l) - srl=0.0_wp_ - rl2=1.0_wp_-lam*bv(1)/bmmx - rl0=0.0_wp_ - if(rl2.gt.0) rl0=sqrt(rl2) - do inc=1,npoints-1 - rl2=1.0_wp_-lam*bv(inc+1)/bmmx - rl=0.0_wp_ - if(rl2.gt.0) rl=sqrt(rl2) - srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) - rl0=rl - end do - srl=srl/anorm - dhlam=0.5_wp_/srl - fc=fc+lam/srl*weights(l) - if(l.eq.nlam) then - fhlam(jp,l)=0.0_wp_ - ffhlam(nlam*(jp-1)+l)=0.0_wp_ - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - else - shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam - fhlam(jp,l)=shlam - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - end if - end do - fc=0.75_wp_*b2av/bmmx**2*fc*dlam - ffc(jp)=fc - - ccfh=bmmn/bmmx/fc - do l=1,nlam - ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) - dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) - end do - end do - - write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// - .' Area Vol |I_pl| qq fc ratioJa ratioJb dadr dvdr' - - qqv(1)=qqv(2) - vajphiav(1)=vajphiav(2) - do jp=1,npsi - rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) - if(jp.eq.npsi) then - rhotqv(jp)=1.0_wp_ - rpstab(jp)=1.0_wp_ - pstab(jp)=1.0_wp_ - end if - rhot_eq(jp)=frhotor(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), - . vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) - end do - -c rarea=sqrt(varea(npsi)/pi) - -c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi -c used for computations of dP/dV and J_cd -c spline coefficients of rhot - - iopt=0 - call difcs(rpstab,vvol,npsi,iopt,cvol,ier) - iopt=0 - call difcs(rpstab,rbav,npsi,iopt,crbav,ier) - iopt=0 - call difcs(rpstab,rri,npsi,iopt,crri,ier) - iopt=0 - call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) - iopt=0 - call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) - iopt=0 - call difcs(rpstab,vratja,npsi,iopt,cratja,ier) - iopt=0 - call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) - iopt=0 - call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) - iopt=0 - call difcs(rpstab,varea,npsi,iopt,carea,ier) - iopt=0 - call difcs(rpstab,ffc,npsi,iopt,cfc,ier) - iopt=0 - call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) - iopt=0 - call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) - iopt=0 - call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) - -c spline interpolation of H(lambda,rhop) - - iopt=0 - s=0.0_wp_ - call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, - . zero,one,zero,one,ksp,ksp,s, - . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) - njpt=njp - nlmt=nlm - - return - 99 format(20(1x,e12.5)) - end +C subroutine flux_average_an +C use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv +C use magsurf_data +C use equilibrium, only : btrcen,frhotor,btaxis,rmaxis,zmaxis, +C . fq,aminor,equian +C use simplespline, only : difcs +C use dierckx, only : regrid,coeff_parder +C implicit none +Cc local constants +C integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, +C . njest=nnintp+ksp+1,nlest=nlam+ksp+1, +C . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, +C . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam +Cc local variables +C integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr +C integer, dimension(kwrk) :: iwrk +C real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, +C . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height2,r2iav,currp,area, +C . volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb,dlp, +C . drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0,shlam, +C . srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi,bphi,brr,bzz,riav,fp, +C . dadr,dbvcdc13,dbvcdc31,dvdr,height,phitedge,qqan,rn +C real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl +C real(wp_), dimension(2*ncnt) :: dlpv +C real(wp_), dimension(2*ncnt+1) :: bv,bpv +C real(wp_), dimension(nlam) :: alam,weights +C real(wp_), dimension(nnintp,nlam) :: fhlam +C real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam +C real(wp_), dimension(lwrk) :: wrk +C real(wp_), dimension(:), allocatable :: rctemp,zctemp +C real(wp_) :: psinv !unused. can be removed when equian is modified +C !to accept optional arguments +Cc common/external functions/variables +C real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv +Cc +C common/fpol/fpolv,dfpv +C common/derip1/dpsidr,dpsidz +C common/derip2/ddpsidrr,ddpsidzz,ddpsidrz +Cc +C npsi=nnintp +C ninpr=(npsi-1)/10 +C npoints = 2*ncnt+1 +Cc +C call alloc_surfvec(ierr) +C if(allocated(tjp)) deallocate(tjp) +C if(allocated(tlm)) deallocate(tlm) +C if(allocated(ch)) deallocate(ch) +C allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), +C . rctemp(npoints),zctemp(npoints),stat=ierr) +C if (ierr.ne.0) return +Cc +Cc computation of flux surface averaged quantities +C +C phitedge=pi*aminor**2*btaxis +C +C write(71,*)' #i psin R z' +C +C dlam=1.0_wp_/dble(nlam-1) +C do l=1,nlam-1 +C alam(l)=dble(l-1)*dlam +C fhlam(1,l)=sqrt(1.0_wp_-alam(l)) +C ffhlam(l)=fhlam(1,l) +C dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) +C weights(l)=1.0_wp_ +C end do +C weights(1)=0.5_wp_ +C weights(nlam)=0.5_wp_ +C alam(nlam)=1.0_wp_ +C fhlam(1,nlam)=0.0_wp_ +C ffhlam(nlam)=0.0_wp_ +C dffhlam(nlam)=-99999.0_wp_ +C +C jp=1 +C anorm=2.0_wp_*pi*rmaxis/abs(btaxis) +C b2av=btaxis**2 +C dvdpsi=2.0_wp_*pi*anorm +C dadpsi=2.0_wp_*pi/abs(btaxis) +C ratio_cdator=abs(btaxis/btrcen) +C ratio_cdbtor=1.0_wp_ +C ratio_pltor=1.0_wp_ +C qq=1.0_wp_ +C fc=1.0_wp_ +C +C psicon(1)=0.0_wp_ +C rcon(1,:)=rmaxis +C zcon(1,:)=zmaxis +C pstab(1)=0.0_wp_ +C rhot_eq(1)=0.0_wp_ +C rpstab(1)=0.0_wp_ +C rhotqv(1)=0.0_wp_ +C vcurrp(1)=0.0_wp_ +C vajphiav(1)=0.0_wp_ +C bmxpsi(1)=abs(btaxis) +C bmnpsi(1)=abs(btaxis) +C bav(1)=abs(btaxis) +C rbav(1)=1.0_wp_ +C rri(1)=rmaxis +C varea(1)=0.0_wp_ +C vvol(1)=0.0_wp_ +C vratjpl(1)=ratio_pltor +C vratja(1)=ratio_cdator +C vratjb(1)=ratio_cdbtor +C ffc(1)=fc +C rhot2q=0.0_wp_ +C dadrhotv(1)=0.0_wp_ +C dvdrhotv(1)=0.0_wp_ +C qqv(1)=1.0_wp_ +C +Cc rup=rmaxis +Cc rlw=rmaxis +Cc zup=(zbmax+zmaxis)/2.0_wp_ +Cc zlw=(zmaxis+zbmin)/2.0_wp_ +C +C do jp=2,npsi +C height=dble(jp-1)/dble(npsi-1) +C psicon(jp)=height +C if(jp.eq.npsi) height=0.9999_wp_ +C height2=height*height +C ipr=0 +C jpr=mod(jp,ninpr) +C if(jpr.eq.1) ipr=1 +C call contours_psi_an(height2,rctemp,zctemp,ipr) +C rcon(jp,:) = rctemp +C zcon(jp,:) = zctemp +Cc +C r2iav=0.0_wp_ +C anorm=0.0_wp_ +C dadpsi=0.0_wp_ +C currp=0.0_wp_ +C b2av=0.0_wp_ +C area=0.0_wp_ +C volume=0.0_wp_ +C ajphiav=0.0_wp_ +C bbav=0.0_wp_ +C bmmx=-1.0e+30_wp_ +C bmmn=1.0e+30_wp_ +C +C call equian(rctemp(1),zctemp(1),psinv) +C dbvcdc13=-ddpsidzz/rctemp(1) +C dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) +C ajphi=ccj*(dbvcdc13-dbvcdc31) +C brr=-dpsidz/rctemp(1) +C bzz= dpsidr/rctemp(1) +C bphi=fpolv/rctemp(1) +C btot0=sqrt(bphi**2+brr**2+bzz**2) +C bpoloid0=sqrt(brr**2+bzz**2) +C bv(1)=btot0 +C bpv(1)=bpoloid0 +C rpsim0=rctemp(1) +C +C do inc=1,npoints-1 +C inc1=inc+1 +C dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) +C dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) +C dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ +C . (zctemp(inc1)-zctemp(inc))**2) +C drc=(rctemp(inc1)-rctemp(inc)) +C +Cc compute length, area and volume defined by psi=height^2 +C +C ph=0.5_wp_*(dla+dlb+dlp) +C area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) +C area=area+sqrt(area2) +C rzp=rctemp(inc1)*zctemp(inc1) +C rz=rctemp(inc)*zctemp(inc) +C volume=pi*(rzp+rz)*drc+volume +C +Cc compute line integral on the contour psi=height^2 +C +C rpsim=rctemp(inc1) +C zpsim=zctemp(inc1) +C call equian(rpsim,zpsim,psinv) +C brr=-dpsidz/rpsim +C bzz= dpsidr/rpsim +C dbvcdc13=-ddpsidzz/rpsim +C dbvcdc31= ddpsidrr/rpsim-bzz/rpsim +C ajphi=ccj*(dbvcdc13-dbvcdc31) +C bphi=fpolv/rpsim +C btot=sqrt(bphi**2+brr**2+bzz**2) +C bpoloid=sqrt(brr**2+bzz**2) +C dlpv(inc)=dlp +C bv(inc1)=btot +C bpv(inc1)=bpoloid +C +C dlph=0.5_wp_*dlp +C anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) +C dadpsi=dadpsi+dlph* +C . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) +C currp=currp+dlph*(bpoloid+bpoloid0) +C b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) +C bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) +C r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) +C . +1.0_wp_/(bpoloid0*rpsim0**2)) +C ajphiav=ajphiav+dlph* +C . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) +C +C ajphi0=ajphi +C rpsim0=rpsim +C bpoloid0=bpoloid +C btot0=btot +Cc computation maximum/minimum B values on given flux surface +C +C if(btot.le.bmmn) bmmn=btot +C if(btot.ge.bmmx) bmmx=btot +C end do +C +Cc bav= [T] , b2av= [T^2] , rbav=/b_min +Cc anorm = int d l_p/B_p = dV/dpsi/(2pi) +Cc r2iav=<1/R^2> [m^-2] , +Cc riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), +Cc rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] +Cc currp = plasma current within psi=const +C +C bbav=bbav/anorm +C r2iav=r2iav/anorm +C dvdpsi=2.0_wp_*pi*anorm +C riav=dadpsi/anorm +C b2av=b2av/anorm +C vcurrp(jp)=ccj*currp +C vajphiav(jp)=ajphiav/dadpsi +C +Cc area == varea, volume == vvol +Cc flux surface minor radius == (area/pi)^1/2 +Cc ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 +Cc ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / +Cc ratio_pltor = Jcd_||/J_phi Jcd_|| = +C +C pstab(jp)=height2 +C rpstab(jp)=height +C vvol(jp)=abs(volume) +C varea(jp)=area +C bav(jp)=bbav +C rbav(jp)=bbav/bmmn +C bmxpsi(jp)=bmmx +C bmnpsi(jp)=bmmn +C rri(jp)=bav(jp)/abs(fpolv*r2iav) +C ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) +C ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) +C ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) +C vratjpl(jp)=ratio_pltor +C vratja(jp)=ratio_cdator +C vratjb(jp)=ratio_cdbtor +C qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) +C qqv(jp)=qq +C +C rn=frhotor(sqrt(pstab(jp))) +C qqan=fq(pstab(jp)) +C +C dadr=2.0_wp_*pi*rn*aminor**2 +C dvdr=4.0_wp_*pi**2*rn*rmaxis*aminor**2 +C +Cc dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi +Cc dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi +C dadrhotv(jp)=phitedge*rn*dadpsi/pi/qqan +C dvdrhotv(jp)=phitedge*rn*dvdpsi/pi/qqan +C +Cc computation of rhot from calculated q profile +C rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) +C . /dble(npsi-1) +Cc print*,jp,rhot2q,qqv(jp),rpstab(jp),qqv(jp-1),rpstab(jp-1) +C rhotqv(jp)=sqrt(rhot2q) +Cc rhotqv(jp)=rn +Cc +C write(57,99) rpstab(jp),rn,rhotqv(jp),qqv(jp),qqan,r2iav, +C . dadr,dadrhotv(jp),dadpsi,dvdpsi,fpolv +C +Cc computation of fraction of circulating/trapped fraction fc, ft +Cc and of function H(lambda,rhop) +Cc ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ +C +C fc=0.0_wp_ +C shlam=0.0_wp_ +C do l=nlam,1,-1 +C lam=alam(l) +C srl=0.0_wp_ +C rl2=1.0_wp_-lam*bv(1)/bmmx +C rl0=0.0_wp_ +C if(rl2.gt.0) rl0=sqrt(rl2) +C do inc=1,npoints-1 +C rl2=1.0_wp_-lam*bv(inc+1)/bmmx +C rl=0.0_wp_ +C if(rl2.gt.0) rl=sqrt(rl2) +C srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) +C rl0=rl +C end do +C srl=srl/anorm +C dhlam=0.5_wp_/srl +C fc=fc+lam/srl*weights(l) +C if(l.eq.nlam) then +C fhlam(jp,l)=0.0_wp_ +C ffhlam(nlam*(jp-1)+l)=0.0_wp_ +C dffhlam(nlam*(jp-1)+l)=-dhlam +C dhlam0=dhlam +C else +C shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam +C fhlam(jp,l)=shlam +C dffhlam(nlam*(jp-1)+l)=-dhlam +C dhlam0=dhlam +C end if +C end do +C fc=0.75_wp_*b2av/bmmx**2*fc*dlam +C ffc(jp)=fc +C +C ccfh=bmmn/bmmx/fc +C do l=1,nlam +C ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) +C dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) +C end do +C end do +C +C write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// +C .' Area Vol |I_pl| qq fc ratioJa ratioJb dadr dvdr' +C +C qqv(1)=qqv(2) +C vajphiav(1)=vajphiav(2) +C do jp=1,npsi +C rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) +C if(jp.eq.npsi) then +C rhotqv(jp)=1.0_wp_ +C rpstab(jp)=1.0_wp_ +C pstab(jp)=1.0_wp_ +C end if +C rhot_eq(jp)=frhotor(sqrt(pstab(jp))) +C write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), +C . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), +C . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), +C . vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) +C end do +C +Cc rarea=sqrt(varea(npsi)/pi) +C +Cc spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi +Cc used for computations of dP/dV and J_cd +Cc spline coefficients of rhot +C +C iopt=0 +C call difcs(rpstab,vvol,npsi,iopt,cvol,ier) +C iopt=0 +C call difcs(rpstab,rbav,npsi,iopt,crbav,ier) +C iopt=0 +C call difcs(rpstab,rri,npsi,iopt,crri,ier) +C iopt=0 +C call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) +C iopt=0 +C call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) +C iopt=0 +C call difcs(rpstab,vratja,npsi,iopt,cratja,ier) +C iopt=0 +C call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) +C iopt=0 +C call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) +C iopt=0 +C call difcs(rpstab,varea,npsi,iopt,carea,ier) +C iopt=0 +C call difcs(rpstab,ffc,npsi,iopt,cfc,ier) +C iopt=0 +C call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) +C iopt=0 +C call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) +C iopt=0 +C call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) +C +Cc spline interpolation of H(lambda,rhop) +C +C iopt=0 +C s=0.0_wp_ +C call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, +C . zero,one,zero,one,ksp,ksp,s, +C . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) +C njpt=njp +C nlmt=nlm +C +C return +C 99 format(20(1x,e12.5)) +C end c c c @@ -2183,6 +2183,8 @@ c fk1(ieq)=ypwrk(ieq,j,k) yy(ieq)=y(ieq)+fk1(ieq)*hh end do + + if(igrad.eq.1) then gr2=grad2(j,k) do iv=1,3 dgr2(iv)=dgrad2v(iv,j,k) @@ -2191,6 +2193,8 @@ c ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do + end if + call fwork(sox,xgcn,bres,yy,fk2) c do ieq=1,ndim @@ -2530,7 +2534,8 @@ c rrm=1.0e-2_wp_*rr c if(iequil.eq.1) then - call equian(rrm,zzm,psinv) + call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, + . ddpsidrr,ddpsidzz,ddpsidrz) end if c if(iequil.eq.2) then @@ -2638,72 +2643,72 @@ c c c c - subroutine tor_curr_psi(h,ajphi) - use const_and_precisions, only : wp_ - use equilibrium, only : zmaxis - implicit none -c arguments - real(wp_) :: h,ajphi -c local variables - real(wp_) :: r1,r2 -c - call psi_raxis(h,r1,r2) - call tor_curr(r2,zmaxis,ajphi) -c - return - end -c -c -c - subroutine tor_curr(rpsim,zpsim,ajphi) - use const_and_precisions, only : wp_,ccj=>mu0inv - use equilibrium, only : equinum_psi - implicit none -c arguments - real(wp_) :: rpsim,zpsim,ajphi -c local variables - real(wp_) :: bzz,dbvcdc13,dbvcdc31 - real(wp_) :: dpsidr,ddpsidrr,ddpsidzz -c - call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,ddpsidrr=ddpsidrr, - . ddpsidzz=ddpsidzz) - bzz= dpsidr/rpsim - dbvcdc13=-ddpsidzz/rpsim - dbvcdc31= ddpsidrr/rpsim-bzz/rpsim - ajphi=ccj*(dbvcdc13-dbvcdc31) -c - return - end -c -c -c - subroutine psi_raxis(h,r1,r2) - use const_and_precisions, only : wp_ - use equilibrium, only : psiant,psinop,zmaxis,nsr, - . nsz,cc=>cceq,tr,tz,kspl - use dierckx, only : profil,sproota - implicit none -c local constants - integer, parameter :: mest=4 -c arguments - real(wp_) :: h,r1,r2 -c local variables - integer :: iopt,ier,m - real(wp_) :: zc,val - real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(nsr) :: czc -c - iopt=1 - zc=zmaxis - call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) - if(ier.gt.0) print*,' profil =',ier - val=h*psiant+psinop - call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) - r1=zeroc(1) - r2=zeroc(2) -c - return - end +C subroutine tor_curr_psi(h,ajphi) +C use const_and_precisions, only : wp_ +C use equilibrium, only : zmaxis +C implicit none +Cc arguments +C real(wp_) :: h,ajphi +Cc local variables +C real(wp_) :: r1,r2 +Cc +C call psi_raxis(h,r1,r2) +C call tor_curr(r2,zmaxis,ajphi) +Cc +C return +C end +Cc +Cc +Cc +C subroutine tor_curr(rpsim,zpsim,ajphi) +C use const_and_precisions, only : wp_,ccj=>mu0inv +C use equilibrium, only : equinum_psi +C implicit none +Cc arguments +C real(wp_) :: rpsim,zpsim,ajphi +Cc local variables +C real(wp_) :: bzz,dbvcdc13,dbvcdc31 +C real(wp_) :: dpsidr,ddpsidrr,ddpsidzz +Cc +C call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,ddpsidrr=ddpsidrr, +C . ddpsidzz=ddpsidzz) +C bzz= dpsidr/rpsim +C dbvcdc13=-ddpsidzz/rpsim +C dbvcdc31= ddpsidrr/rpsim-bzz/rpsim +C ajphi=ccj*(dbvcdc13-dbvcdc31) +Cc +C return +C end +Cc +Cc +Cc +C subroutine psi_raxis(h,r1,r2) +C use const_and_precisions, only : wp_ +C use equilibrium, only : psiant,psinop,zmaxis,nsr, +C . nsz,cc=>cceq,tr,tz,kspl +C use dierckx, only : profil,sproota +C implicit none +Cc local constants +C integer, parameter :: mest=4 +Cc arguments +C real(wp_) :: h,r1,r2 +Cc local variables +C integer :: iopt,ier,m +C real(wp_) :: zc,val +C real(wp_), dimension(mest) :: zeroc +C real(wp_), dimension(nsr) :: czc +Cc +C iopt=1 +C zc=zmaxis +C call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) +C if(ier.gt.0) print*,' profil =',ier +C val=h*psiant+psinop +C call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) +C r1=zeroc(1) +C r2=zeroc(2) +Cc +C return +C end c c c @@ -3408,53 +3413,76 @@ c end c c +Cc +C subroutine valpsisplpec(rhop,voli,areai) +C use const_and_precisions, only : wp_ +C use utils, only : locate +C use simplespline, only :spli +C use magsurf_data, only : rpstab,cvol,carea,npsi +C implicit none +Cc arguments +C real(wp_), intent(in) :: rhop +C real(wp_), intent(out) :: voli,areai +Cc local variables +C integer :: ip +C real(wp_) :: drh +Cc +C call locate(rpstab,npsi,rhop,ip) +C ip=min(max(1,ip),npsi-1) +C drh=rhop-rpstab(ip) +Cc +C voli=spli(cvol,npsi,ip,drh) +C areai=spli(carea,npsi,ip,drh) +Cc +C return +C end c - subroutine valpsisplpec(rhop,voli,areai) - use const_and_precisions, only : wp_ - use utils, only : locate - use simplespline, only :spli - use magsurf_data, only : rpstab,cvol,carea,npsi - implicit none -c arguments - real(wp_), intent(in) :: rhop - real(wp_), intent(out) :: voli,areai -c local variables - integer :: ip - real(wp_) :: drh -c - call locate(rpstab,npsi,rhop,ip) - ip=min(max(1,ip),npsi-1) - drh=rhop-rpstab(ip) -c - voli=spli(cvol,npsi,ip,drh) - areai=spli(carea,npsi,ip,drh) -c - return - end -c -c -c - subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) - use const_and_precisions, only : wp_ - use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi - use utils, only : locate - use simplespline, only :spli - implicit none -c arguments - real(wp_) :: rpsi,ratjai,ratjbi,ratjpli -c local variables - integer :: ip - real(wp_) :: dps -c - call locate(rpstab,npsi,rpsi,ip) - ip=min(max(1,ip),npsi-1) - dps=rpsi-rpstab(ip) - ratjai=spli(cratja,npsi,ip,dps) - ratjbi=spli(cratjb,npsi,ip,dps) - ratjpli=spli(cratjpl,npsi,ip,dps) -c - return - end +C subroutine valpsisplpec(rhop,voli,areai) +C use const_and_precisions, only : wp_ +C use utils, only : locate +C use simplespline, only :spli +C use magsurf_data, only : rpstab,cvol,carea,npsi +C implicit none +Cc arguments +C real(wp_), intent(in) :: rhop +C real(wp_), intent(out) :: voli,areai +Cc local variables +C integer :: ip +C real(wp_) :: drh +Cc +C call locate(rpstab,npsi,rhop,ip) +C ip=min(max(1,ip),npsi-1) +C drh=rhop-rpstab(ip) +Cc +C voli=spli(cvol,npsi,ip,drh) +C areai=spli(carea,npsi,ip,drh) +Cc +C return +C end +Cc +Cc +Cc +C subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) +C use const_and_precisions, only : wp_ +C use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi +C use utils, only : locate +C use simplespline, only :spli +C implicit none +Cc arguments +C real(wp_) :: rpsi,ratjai,ratjbi,ratjpli +Cc local variables +C integer :: ip +C real(wp_) :: dps +Cc +C call locate(rpstab,npsi,rpsi,ip) +C ip=min(max(1,ip),npsi-1) +C dps=rpsi-rpstab(ip) +C ratjai=spli(cratja,npsi,ip,dps) +C ratjbi=spli(cratjb,npsi,ip,dps) +C ratjpli=spli(cratjpl,npsi,ip,dps) +Cc +C return +C end c c c @@ -3467,6 +3495,8 @@ c use dispersion, only : harmnumber, warmdisp use equilibrium, only : rmaxis,sgnbphi use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl + use magsurf_data, only : fluxval + implicit none c arguments integer, intent(in) :: i,j,k @@ -3529,7 +3559,8 @@ c effjcd=0.0_wp_ c tekev=temp(psinv) - call valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) + call fluxval(rhop,dervol=dervoli,rri=rrii,rbav=rbavi, + . bmn=bmni,bmx=bmxi,fc=fci) c if(tekev.gt.0.0_wp_.and.tau0.le.taucr) then c @@ -3605,31 +3636,31 @@ c c c c - subroutine valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) - use const_and_precisions, only : wp_ - use utils, only : locate - use simplespline, only :spli,splid - use magsurf_data, only : rpstab,crri,crbav,cvol,cbmn,cbmx,cfc,npsi - implicit none -c arguments - real(wp_), intent(in) :: rhop - real(wp_), intent(out) :: rrii,rbavi,dervoli,bmni,bmxi,fci -c local variables - integer :: ip - real(wp_) :: drh -c - call locate(rpstab,npsi,rhop,ip) - ip=min(max(1,ip),npsi-1) - drh=rhop-rpstab(ip) - rrii=spli(crri,npsi,ip,drh) - rbavi=spli(crbav,npsi,ip,drh) - dervoli=splid(cvol,npsi,ip,drh) - bmni=spli(cbmn,npsi,ip,drh) - bmxi=spli(cbmx,npsi,ip,drh) - fci=spli(cfc,npsi,ip,drh) -c - return - end +C subroutine valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) +C use const_and_precisions, only : wp_ +C use utils, only : locate +C use simplespline, only :spli,splid +C use magsurf_data, only : rpstab,crri,crbav,cvol,cbmn,cbmx,cfc,npsi +C implicit none +Cc arguments +C real(wp_), intent(in) :: rhop +C real(wp_), intent(out) :: rrii,rbavi,dervoli,bmni,bmxi,fci +Cc local variables +C integer :: ip +C real(wp_) :: drh +Cc +C call locate(rpstab,npsi,rhop,ip) +C ip=min(max(1,ip),npsi-1) +C drh=rhop-rpstab(ip) +C rrii=spli(crri,npsi,ip,drh) +C rbavi=spli(crbav,npsi,ip,drh) +C dervoli=splid(cvol,npsi,ip,drh) +C bmni=spli(cbmn,npsi,ip,drh) +C bmxi=spli(cbmx,npsi,ip,drh) +C fci=spli(cfc,npsi,ip,drh) +Cc +C return +C end c c c @@ -3724,6 +3755,7 @@ c use const_and_precisions, only : wp_,zero,one,izero use gray_params, only : nnd use equilibrium, only : frhotor,frhopol + use magsurf_data, only : fluxval implicit none integer :: it,ipec real(wp_) :: drt,rt,rt1,rhop1 @@ -3740,7 +3772,7 @@ c ipec negative read input grid c c ipec=+/-1 rho_pol grid c ipec=+/-2 rho_tor grid - + voli0=zero areai0=zero rtabpsi1(0)=zero @@ -3773,13 +3805,14 @@ c radial coordinate of i-(i+1) interval mid point c psi grid at mid points, dimension nnd+1, for use in pec_tab rtabpsi1(it)=rhop1**2 - call valpsisplpec(rhop1,voli1,areai1) + call fluxval(rhop1,area=areai1,vol=voli1) dvol(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 - call ratioj(rhop_tab(it),ratjai,ratjbi,ratjpli) + call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi, + . ratjpl=ratjpli) ratjav(it)=ratjai ratjbv(it)=ratjbi ratjplv(it)=ratjpli @@ -3965,6 +3998,7 @@ c radial average values over power and current density profile use const_and_precisions, only : wp_,one,zero,izero,pi use gray_params, only : nnd,ieccd use equilibrium, only : frhopol + use magsurf_data, only : fluxval implicit none real(wp_), intent(in) :: pabs,currt real(wp_), intent(out) :: rhotpav,rhotjava @@ -3975,10 +4009,9 @@ c radial average values over power and current density profile real(wp_) :: ratjamx,ratjbmx,ratjplmx real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv - real(wp_), external :: fdadrhot,fdvdrhot real(wp_) :: sccsa - real(wp_) :: rhotjav,rhot2pav,rhot2java + real(wp_) :: rhotjav,rhot2pav,rhot2java,dvdrhotav,dadrhotava rhotpav=zero rhot2pav=zero @@ -4010,10 +4043,15 @@ c factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile drhotjava=sqrt(8._wp_*(rhot2java-rhotjava**2)) rhoppav=frhopol(rhotpav) + call fluxval(rhoppav,dvdrhot=dvdrhotav,ratja=ratjamx, + . ratjb=ratjbmx,ratjpl=ratjplmx) + rhopjava=frhopol(rhotjava) - dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhoppav)) - ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhopjava)) - call ratioj(rhopjava,ratjamx,ratjbmx,ratjplmx) + call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx, + . ratjb=ratjbmx,ratjpl=ratjplmx) + + dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav) + ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava) call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp) if(ieccd.ne.0) then diff --git a/src/graycore.f90 b/src/graycore.f90 index 73d2220..6e3949c 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -13,6 +13,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl + use magsurf_data, only : flux_average use beamdata, only : init_rtr, dealloc_beam, nstep, dst use reflections, only : set_lim implicit none @@ -49,7 +50,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & if(iequil<2) then call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) - call flux_average_an + call flux_average else call set_eqspl(rv,zv,psin, psia, psinr,fpol, eqp%ssplps,eqp%ssplf, rvac, & rax,zax, rbnd,zbnd, eqp%ixp) @@ -97,6 +98,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & call prfile call paraminit call vectinit + if(igrad==0) then call ic_rt(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) @@ -104,20 +106,22 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & call ic_gb(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) end if -! if(ierr==0) return -! beam/ray propagation + +! beam/ray propagation + ! st0=0.0_wp_ ! if(index_rt.gt.1) st0=strfl11 do i=1,nstep istep=i st=i*dst !+st0 -! advance one step +! advance one step call rkint4(sox,xgcn,bres) -! calculations after one step: +! calculations after one step: call after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ierr) if(istop.eq.1) exit end do -! postprocessing + +! postprocessing call after_gray_integration(pec,icd,dpdv,jcd) ! ======= main loop END ====== diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index b34554f..fa69d23 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -2,23 +2,23 @@ module magsurf_data use const_and_precisions, only : wp_ implicit none - INTEGER, SAVE :: npsi, npoints !# sup mag, # punti per sup - INTEGER, SAVE :: njpt, nlmt + integer, save :: npsi, npoints !# sup mag, # punti per sup + integer, save :: njpt, nlmt - REAL(wp_), SAVE :: rarea + real(wp_), save :: rarea - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psicon,pstab,rhot_eq, & + real(wp_), dimension(:), allocatable, save :: psicon,pstab,rhot_eq, & rhotqv,bav,varea,vcurrp,vajphiav,qqv,ffc,vratja,vratjb - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rpstab - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: vvol,rri,rbav,bmxpsi,bmnpsi - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tjp,tlm,ch,ch01 + real(wp_), dimension(:), allocatable, save :: rpstab + real(wp_), dimension(:), allocatable, save :: vvol,rri,rbav,bmxpsi,bmnpsi + real(wp_), dimension(:), allocatable, save :: tjp,tlm,ch,ch01 - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: rcon,zcon - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cdadrhot,cdvdrhot + real(wp_), dimension(:,:), allocatable, save :: rcon,zcon + real(wp_), dimension(:,:), allocatable, save :: cdadrhot,cdvdrhot - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cvol,crri,crbav,cbmx,cbmn,carea,cfc - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhotq - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cratja,cratjb,cratjpl + real(wp_), dimension(:,:), allocatable, save :: cvol,crri,crbav,cbmx,cbmn,carea,cfc + real(wp_), dimension(:,:), allocatable, save :: crhotq + real(wp_), dimension(:,:), allocatable, save :: cratja,cratjb,cratjpl contains @@ -103,4 +103,383 @@ contains if(allocated(cratjb)) deallocate(cratjb) end subroutine dealloc_surfvec + + + subroutine flux_average + use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv + use gray_params, only : iequil + use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & + equian,equinum_psi,bfield,frhotor,fq,tor_curr + use simplespline, only : difcs + use dierckx, only : regrid,coeff_parder + use utils, only : get_free_unit + implicit none +! local constants + integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, & + njest=nnintp+ksp+1,nlest=nlam+ksp+1, & + lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, & + kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam +! local variables + integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr,u56,unit + integer, dimension(kwrk) :: iwrk + real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, & + ratio_cdbtor,ratio_pltor,fc,height,height2,r2iav,currp, & + area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, & + dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, & + shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, & + bphi,brr,bzz,riav,fp,psinjp,rhopjp,rhotjp,qq,rup,rlw,zup,zlw + real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl + real(wp_), dimension(2*ncnt) :: dlpv + real(wp_), dimension(2*ncnt+1) :: bv,bpv + real(wp_), dimension(nlam) :: alam,weights + real(wp_), dimension(nnintp,nlam) :: fhlam + real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam + real(wp_), dimension(lwrk) :: wrk + real(wp_), dimension(:), allocatable :: rctemp,zctemp +! common/external functions/variables + real(wp_) :: fpolv,ddpsidrr,ddpsidzz + + u56=get_free_unit() + open(file='f56.txt',unit=u56) + + npsi=nnintp + ninpr=(npsi-1)/10 + npoints = 2*ncnt+1 + + call alloc_surfvec(ierr) + if(allocated(tjp)) deallocate(tjp) + if(allocated(tlm)) deallocate(tlm) + if(allocated(ch)) deallocate(ch) + allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), & + rctemp(npoints),zctemp(npoints),stat=ierr) + if (ierr.ne.0) return + +! computation of flux surface averaged quantities + + write(71,*)' #i psin R z' + + dlam=1.0_wp_/dble(nlam-1) + do l=1,nlam-1 + alam(l)=dble(l-1)*dlam + fhlam(1,l)=sqrt(1.0_wp_-alam(l)) + ffhlam(l)=fhlam(1,l) + dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) + weights(l)=1.0_wp_ + end do + weights(1)=0.5_wp_ + weights(nlam)=0.5_wp_ + alam(nlam)=1.0_wp_ + fhlam(1,nlam)=0.0_wp_ + ffhlam(nlam)=0.0_wp_ + dffhlam(nlam)=-99999.0_wp_ + + jp=1 + anorm=2.0_wp_*pi*rmaxis/abs(btaxis) + dvdpsi=2.0_wp_*pi*anorm + dadpsi=2.0_wp_*pi/abs(btaxis) + b2av=btaxis**2 + ratio_cdator=abs(btaxis/btrcen) + ratio_cdbtor=1.0_wp_ + ratio_pltor=1.0_wp_ + fc=1.0_wp_ + if(iequil < 2) then + call equian(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) + else + call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) + end if + qq=btaxis/sqrt(ddpsidrr*ddpsidzz) + ajphiav=-ccj*(ddpsidrr+ddpsidzz)/rmaxis + + psicon(1)=0.0_wp_ + rcon(1,:)=rmaxis + zcon(1,:)=zmaxis + pstab(1)=0.0_wp_ + rpstab(1)=0.0_wp_ + vcurrp(1)=0.0_wp_ + vajphiav(1)=ajphiav + bmxpsi(1)=abs(btaxis) + bmnpsi(1)=abs(btaxis) + bav(1)=abs(btaxis) + rbav(1)=1.0_wp_ + rri(1)=rmaxis + varea(1)=0.0_wp_ + vvol(1)=0.0_wp_ + vratjpl(1)=ratio_pltor + vratja(1)=ratio_cdator + vratjb(1)=ratio_cdbtor + ffc(1)=fc + qqv(1)=qq + dadrhotv(1)=0.0_wp_ + dvdrhotv(1)=0.0_wp_ + + rup=rmaxis + rlw=rmaxis + zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ + zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ + + do jp=2,npsi + height=dble(jp-1)/dble(npsi-1) + if(jp.eq.npsi) height=0.9999_wp_ + ipr=0 + jpr=mod(jp,ninpr) + if(jpr.eq.1) ipr=1 + rhopjp=height + psinjp=height*height + rhotjp=frhotor(rhopjp) + psicon(jp)=height + + if(iequil<2) then + call contours_psi_an(psinjp,rctemp,zctemp,ipr) + else + call contours_psi(psinjp,rup,zup,rlw,zlw,rctemp,zctemp,ipr) + end if + rcon(jp,:) = rctemp + zcon(jp,:) = zctemp + + r2iav=0.0_wp_ + anorm=0.0_wp_ + dadpsi=0.0_wp_ + currp=0.0_wp_ + b2av=0.0_wp_ + area=0.0_wp_ + volume=0.0_wp_ + ajphiav=0.0_wp_ + bbav=0.0_wp_ + bmmx=-1.0e+30_wp_ + bmmn=1.0e+30_wp_ + + call tor_curr(rctemp(1),zctemp(1),ajphi0) + call bfield(rctemp(1),zctemp(1),bphi,br=brr,bz=bzz) + fpolv=bphi*rctemp(1) + btot0=sqrt(bphi**2+brr**2+bzz**2) + bpoloid0=sqrt(brr**2+bzz**2) + bv(1)=btot0 + bpv(1)=bpoloid0 + rpsim0=rctemp(1) + + do inc=1,npoints-1 + inc1=inc+1 + dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) + dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) + dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2) + drc=(rctemp(inc1)-rctemp(inc)) + +! compute length, area and volume defined by psi=psinjp=height^2 + ph=0.5_wp_*(dla+dlb+dlp) + area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) + area=area+sqrt(area2) + rzp=rctemp(inc1)*zctemp(inc1) + rz=rctemp(inc)*zctemp(inc) + volume=pi*(rzp+rz)*drc+volume + +! compute line integrals on the contour psi=psinjp=height^2 + rpsim=rctemp(inc1) + zpsim=zctemp(inc1) + call bfield(rpsim,zpsim,br=brr,bz=bzz) + call tor_curr(rpsim,zpsim,ajphi) + bphi=fpolv/rpsim + btot=sqrt(bphi**2+brr**2+bzz**2) + bpoloid=sqrt(brr**2+bzz**2) + dlpv(inc)=dlp + bv(inc1)=btot + bpv(inc1)=bpoloid + + dlph=0.5_wp_*dlp + anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) + dadpsi=dadpsi+dlph*(1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) + currp=currp+dlph*(bpoloid+bpoloid0) + b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) + bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) + r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2)+1.0_wp_/(bpoloid0*rpsim0**2)) + ajphiav=ajphiav+dlph*(ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) + + ajphi0=ajphi + rpsim0=rpsim + bpoloid0=bpoloid + btot0=btot + +! computation maximum/minimum B values on given flux surface + if(btot.le.bmmn) bmmn=btot + if(btot.ge.bmmx) bmmx=btot + end do + +! bav= [T] , b2av= [T^2] , rbav=/b_min +! anorm = int d l_p/B_p = dV/dpsi/(2pi) +! r2iav=<1/R^2> [m^-2] , +! riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), +! rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] +! currp = plasma current within psi=const + + bbav=bbav/anorm + r2iav=r2iav/anorm + dvdpsi=2.0_wp_*pi*anorm + riav=dadpsi/anorm + b2av=b2av/anorm + vcurrp(jp)=ccj*currp + vajphiav(jp)=ajphiav/dadpsi + +! area == varea, volume == vvol +! flux surface minor radius == (area/pi)^1/2 +! ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 +! ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / +! ratio_pltor = Jcd_||/J_phi Jcd_|| = + pstab(jp)=psinjp + rpstab(jp)=rhopjp + vvol(jp)=abs(volume) + varea(jp)=area + bav(jp)=bbav + rbav(jp)=bbav/bmmn + bmxpsi(jp)=bmmx + bmnpsi(jp)=bmmn + rri(jp)=bav(jp)/abs(fpolv*r2iav) + ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) + ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) + ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) + vratjpl(jp)=ratio_pltor + vratja(jp)=ratio_cdator + vratjb(jp)=ratio_cdbtor + qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) + qqv(jp)=qq + dadrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dadpsi/pi + dvdrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dvdpsi/pi + +! computation of fraction of circulating/trapped fraction fc, ft +! and of function H(lambda,rhop) +! ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ + fc=0.0_wp_ + shlam=0.0_wp_ + do l=nlam,1,-1 + lam=alam(l) + srl=0.0_wp_ + rl2=1.0_wp_-lam*bv(1)/bmmx + rl0=0.0_wp_ + if(rl2.gt.0) rl0=sqrt(rl2) + do inc=1,npoints-1 + rl2=1.0_wp_-lam*bv(inc+1)/bmmx + rl=0.0_wp_ + if(rl2.gt.0) rl=sqrt(rl2) + srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) + rl0=rl + end do + srl=srl/anorm + dhlam=0.5_wp_/srl + fc=fc+lam/srl*weights(l) + if(l.eq.nlam) then + fhlam(jp,l)=0.0_wp_ + ffhlam(nlam*(jp-1)+l)=0.0_wp_ + dffhlam(nlam*(jp-1)+l)=-dhlam + dhlam0=dhlam + else + shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam + fhlam(jp,l)=shlam + dffhlam(nlam*(jp-1)+l)=-dhlam + dhlam0=dhlam + end if + end do + fc=0.75_wp_*b2av/bmmx**2*fc*dlam + ffc(jp)=fc + + ccfh=bmmn/bmmx/fc + do l=1,nlam + ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) + dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) + end do + end do + + write(u56,*)' #rhop rhot || |Bmx| |Bmn| Area Vol |I_pl| fc ratJa ratJb' + + do jp=1,npsi + if(jp.eq.npsi) then + rpstab(jp)=1.0_wp_ + pstab(jp)=1.0_wp_ + end if + rhotjp=frhotor(rpstab(jp)) + write(u56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp), & + varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp),ffc(jp), & + vratja(jp),vratjb(jp) + end do + + close(u56) + +! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs +! used for computations of dP/dV and J_cd + iopt=0 + call difcs(rpstab,vvol,npsi,iopt,cvol,ier) + iopt=0 + call difcs(rpstab,rbav,npsi,iopt,crbav,ier) + iopt=0 + call difcs(rpstab,rri,npsi,iopt,crri,ier) + iopt=0 + call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) + iopt=0 + call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) + iopt=0 + call difcs(rpstab,vratja,npsi,iopt,cratja,ier) + iopt=0 + call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) + iopt=0 + call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) + iopt=0 + call difcs(rpstab,varea,npsi,iopt,carea,ier) + iopt=0 + call difcs(rpstab,ffc,npsi,iopt,cfc,ier) + iopt=0 + call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) + iopt=0 + call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) +! iopt=0 +! call difcs(rpstab,qqv,npsi,iopt,cqq,ier) + +! spline interpolation of H(lambda,rhop) and dH/dlambda + iopt=0 + s=0.0_wp_ + call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam,zero,one,zero,one, & + ksp,ksp,s,njest,nlest,njp,tjp,nlm,tlm,ch,fp, & + wrk,lwrk,iwrk,kwrk,ier) + njpt=njp + nlmt=nlm + + 99 format(20(1x,e12.5)) + + end subroutine flux_average + + + + subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, & + rri,rbav,bmn,bmx,fc,ratja,ratjb,ratjpl) + use const_and_precisions, only : wp_ + use utils, only : locate + use simplespline, only :spli,splid + implicit none +! arguments + real(wp_), intent(in) :: rhop + real(wp_), intent(out), optional :: vol,area,rri,rbav,dervol,bmn,bmx,fc, & + ratja,ratjb,ratjpl,dadrhot,dvdrhot +! local variables + integer :: ip + real(wp_) :: drh + + call locate(rpstab,npsi,rhop,ip) + ip=min(max(1,ip),npsi-1) + drh=rhop-rpstab(ip) + + if (present(area)) area=spli(carea,npsi,ip,drh) + if (present(vol)) vol=spli(cvol,npsi,ip,drh) + + if (present(dervol)) dervol=splid(cvol,npsi,ip,drh) + if (present(dadrhot)) dadrhot=spli(cdadrhot,npsi,ip,drh) + if (present(dvdrhot)) dvdrhot=spli(cdvdrhot,npsi,ip,drh) + + if (present(rri)) rri=spli(crri,npsi,ip,drh) + if (present(rbav)) rbav=spli(crbav,npsi,ip,drh) + if (present(bmn)) bmn=spli(cbmn,npsi,ip,drh) + if (present(bmx)) bmx=spli(cbmx,npsi,ip,drh) + if (present(fc)) fc=spli(cfc,npsi,ip,drh) + + if (present(ratja)) ratja=spli(cratja,npsi,ip,drh) + if (present(ratjb)) ratjb=spli(cratjb,npsi,ip,drh) + if (present(ratjpl)) ratjpl=spli(cratjpl,npsi,ip,drh) + + end subroutine fluxval + end module magsurf_data