diff --git a/src/gray.f b/src/gray.f index e4a2271..55b3a7c 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1,5 +1,4 @@ implicit real*8 (a-h,o-z) - common/istop/istop common/ierr/ierr common/igrad/igrad common/mode/sox @@ -32,7 +31,6 @@ c postprocessing pabstott=pabstot currtott=currtot - powtr=p0mw-pabstot if (ipass.gt.1) then c second pass into plasma @@ -121,16 +119,9 @@ c common/iieq/iequil common/iipr/iprof common/index_rt/index_rt -c - common/p0/p0mw - common/factb/factb common/taumnx/taumn,taumx,pabstot,currtot - common/scal/iscal - common/facttn/factt,factn - common/pardens/dens0,aln1,aln2 common/parqte/te0,dte0,alt1,alt2 - c c print all ray positions in local reference system c @@ -170,13 +161,12 @@ c call pec(pabs,currt) c return - 99 format(20(1x,e16.8e3)) end c c c subroutine after_onestep(i,istop) - use const_and_precisions, only : pi, cvdr=>degree + use const_and_precisions, only : pi use reflections, only : inside implicit real*8 (a-h,o-z) complex*16 ext,eyt,extr,eytr,exin2,eyin2 @@ -191,7 +181,7 @@ c dimension xwcl(3),xvjk(3,jmx,kmx),anvjk(3,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) real*8 rlim(nbb),zlim(nbb) - logical inside_plasma + logical ins_pl,inside_plasma external inside_plasma c common/pcjki/ppabs,ccci @@ -213,16 +203,12 @@ c common/taumnx/taumn,taumx,pabstot,currtot common/xv/xv common/anv/anv - common/cent/btrcen,rcen c - common/p0/p0mw common/pol0/psipol0,chipol0 common/polcof/psipol,chipol common/evt/ext,eyt - common/densbnd/psdbnd common/yyrfl/yyrfl common/powrfl/powrfl - common/dstvac/dstvac common/strfl11/strfl11 common/dsds/dst common/index_rt/index_rt @@ -267,7 +253,6 @@ c exit psjki(j,k,i)=psinv rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 zzm=xv(3)/100.d0 - if(j.eq.1) rrm11=rrm if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then if(psinv.ge.0.and.psinv.le.1.0d0.and. @@ -285,8 +270,8 @@ c exit call print_output(i,j,k) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - if (mod(iop(j,k),2).eq.0 .and. inside_plasma(rrm,zzm)) then + ins_pl=inside_plasma(rrm,zzm) + if (mod(iop(j,k),2).eq.0 .and. ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) @@ -298,7 +283,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ihcd(j,k)=0 end if else if (mod(iop(j,k),2).eq.1.and. - . .not.inside_plasma(rrm,zzm)) then + . .not.ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if @@ -309,7 +294,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc else if (iow(j,k).eq.1 .and. . .not.inside(rlim,zlim,nlim,rrm,zzm)) then iow(j,k)=2 - if (inside_plasma(rrm,zzm)) then + if (ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if @@ -334,7 +319,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ywrk(4:6,j,k)=anv igrad=0 call gwork(j,k) - if (inside_plasma(rrm,zzm)) then + if (ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) if (index_rt.eq.1) ihcd(j,k)=0 @@ -556,14 +541,10 @@ c central ray only begin rhot=1.0d0 dens11=0.0d0 if(psi.ge.0.0d0) then - if(iequil.eq.2) then if (psi.le.1.0d0) rhot=frhotor(psi) bbr=brr bbz=bzz bpol=sqrt(bbr**2+bbz**2) - else - rhot=sqrt(psi) - end if else tekev=0.0d0 akim=0.0d0 @@ -617,7 +598,7 @@ c If(index_rt.eq.1) then - write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// + write(4,*)' #sst R z phi psi rhot ne Te BR Bphi Bz Jphi '// .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(8,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin' @@ -702,6 +683,7 @@ c common/parqq/q0,qa,alq common/parqte/te0,dte0,alt1,alt2 common/zz/Zeff + common/bound/zbmin,zbmax c common/parbres/bres common/densbnd/psdbnd @@ -844,6 +826,13 @@ c if iequil=0,1 rr0m=rr0/1.0d2 zr0m=zr0/1.0d2 rpam=rpa/1.0d2 + rcen=rr0m + btrcen=b0 + zbmin=(zr0-rpa)/100.d0 + zbmax=(zr0+rpa)/100.d0 + b0=b0*factb + call flux_average_an +c call print_prof_an else read(2,*) dummy,dummy,dummy read(2,*) dummy @@ -954,7 +943,7 @@ c versus psi, rhop, rhot call print_prof end if - if (iequil.eq.1) call surf_anal + if (iequil.eq.1) call bres_anal if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 @@ -1025,7 +1014,6 @@ c set simple limiter as two cylindrical walls at rwallm and r00 915 format('# dst nstep : ',1x,es12.5,i5) 916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2)) - 99 format(20(1x,e16.8e3)) end c c @@ -1068,6 +1056,27 @@ c return 111 format(i4,20(1x,e16.8e3)) end + + subroutine bres_anal + implicit real*8(a-h,o-z) + parameter(pi=3.14159265358979d0) + common/parban/b0,rr0m,zr0m,rpam + common/parbres/bres +c +c print resonant B iequil=1 +c + write(70,*)'#i Btot R z' + rres=b0*rr0m/bres + zmx=zr0m+rpam + zmn=zr0m-rpam + do i=1,21 + zzres=zmn+(i-1)*(zmx-zmn)/2.0d1 + write(70,111) i,bres,rres,zzres + end do + + return +111 format(i4,20(1x,e16.8e3)) + end c c subroutine read_beams @@ -1207,7 +1216,6 @@ c parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh) dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) - dimension iwrkd(ldiwrk) parameter(lwrkf=nnw*4+nrest*16) dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw) dimension fpoli(nnw) @@ -1239,7 +1247,6 @@ c common/fpas/fpolas common/rhotsx/rhotsx common/phitedge/phitedge - common/rrtor/rrtor common/cnt/rup,zup,rlw,zlw common/limiter/rlim,zlim,nlim c @@ -1358,7 +1365,7 @@ c enddo do j=1,nz do i=1,nr if(ipsinorm.eq.0) then - psin(i,j)=(psi(i,j)-psiaxis)/(psia0) + psin(i,j)=(psi(i,j)-psiaxis)/psia0 psi(i,j)=psin(i,j)*psia else psi(i,j)=psin(i,j)*psia @@ -1397,32 +1404,32 @@ c write(79,2021) rv(i),zv(j),psin(i,j) c write(79,*) ' ' enddo end if -2021 format(5(1x,e16.9)) +c2021 format(5(1x,e16.9)) c nur=1 nuz=0 - call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, - . ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier) + call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz, + . cc10,lw10,ier) c nur=0 nuz=1 - call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, - . ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier) + call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz, + . cc01,lw01,ier) c nur=2 nuz=0 - call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, - . ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier) + call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz, + . cc20,lw20,ier) c nur=0 nuz=2 - call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, - . ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier) + call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz, + . cc02,lw02,ier) c nur=1 nuz=1 - call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, - . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier) + call coeff_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz, + . cc11,lw11,ier) c c scaling of f_poloidal c @@ -1579,7 +1586,7 @@ c compute B_toroidal on axis c compute normalized rho_tor from eqdsk q profile call rhotor(nr) phitedge=abs(psia)*rhotsx*2*pi - rrtor=sqrt(phitedge/abs(btrcen)/pi) +c rrtor=sqrt(phitedge/abs(btrcen)/pi) call rhopol c print*,rhotsx,phitedge,rrtor,abs(psia) @@ -1644,6 +1651,7 @@ c end if rf=xvec(1) zf=xvec(2) + call equinum_psi(rf,zf) psinvf=psinv return end @@ -1653,18 +1661,25 @@ c subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/pareq1/psia - call equinum(x(1),x(2)) - if (iflag.ne.2) then + + select case(iflag) + case(1) + call equinum_derpsi(x(1),x(2),iflag) fvec(1) = dpsidr/psia - fvec(2) = dpsidz/psia - else + fvec(2) = dpsidz/psia + case(2) + call equinum_derpsi(x(1),x(2),iflag) fjac(1,1) = ddpsidrr/psia fjac(1,2) = ddpsidrz/psia fjac(2,1) = ddpsidrz/psia fjac(2,2) = ddpsidzz/psia - end if + case default + print*,'iflag undefined' + end select + return end c @@ -1693,20 +1708,29 @@ c subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/psival/psinv common/cnpsi/h common/pareq1/psia - call equinum(x(1),x(2)) - if (iflag.ne.2) then + + select case(iflag) + case(1) + call equinum_psi(x(1),x(2)) + call equinum_derpsi(x(1),x(2),iflag) fvec(1) = psinv-h fvec(2) = dpsidr/psia - else + case(2) + ii=iflag+1 + call equinum_derpsi(x(1),x(2),ii) fjac(1,1) = dpsidr/psia fjac(1,2) = dpsidz/psia fjac(2,1) = ddpsidrr/psia fjac(2,2) = ddpsidrz/psia - end if + case default + print*,'iflag undefined' + end select + return end c @@ -1756,6 +1780,25 @@ c return 111 format(12(1x,e12.5)) end + + subroutine print_prof_an + implicit real*8 (a-h,o-z) + parameter(nst=51) + common/dens/dens,ddens + + write(55,*) ' #psi rhot ne Te' + do i=1,nst + psin=dble(i-1)/dble(nst-1) + rhop=sqrt(psin) + rhot=frhotor(psin) + call density(psin) + te=temp(psin) + write(55,111) psin,rhot,dens,te + end do + + return + 111 format(12(1x,e12.5)) + end c c c @@ -2005,7 +2048,9 @@ c if(irt.eq.nr) irt=nr-1 function frhotor(psi) implicit real*8(a-h,o-z) - frhotor=frhotor_eq(psi) + common/iieq/iequil + if(iequil.eq.2) frhotor=frhotor_eq(psi) + if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) return end @@ -2052,7 +2097,7 @@ c spline interpolation of rhopol versus rhotor iopt=0 xb=0.0d0 xe=1.0d0 - ss=0.00001 + ss=0.00001d0 kspl=3 call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp, . trp,crp,rp,wrkp,lwrkp,iwrkp,ier) @@ -2279,7 +2324,6 @@ c dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c - common/pareq1/psia common/pareq1t/psiant,psinop common/coffeqn/nsr,nsz,nsft common/coffeq/cc @@ -2329,7 +2373,6 @@ c end if return 111 format(i6,12(1x,e12.5)) -99 format(12(1x,e12.5)) end c c @@ -2368,16 +2411,13 @@ c common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/eqnn/nr,nz,npp,nintp - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + common/fpol/fpolv,ffpv c common/pstab/rpstab common/flav/vvol,rri,rbav,bmxpsi,bmnpsi common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc common/cratj/cratja,cratjb,cratjpl common/crhotq/crhotq - common/cnt/rup,zup,rlw,zlw - common/bound/zbmin,zbmax - common/rarea/rarea common/phitedge/phitedge common/cdadrhot/cdadrhot common/cdvdrhot/cdvdrhot @@ -2437,6 +2477,7 @@ c computation of flux surface averaged quantities vratjb(1)=ratio_cdbtor ffc(1)=fc rhot2q=0.0d0 + qqv(1)=1.0d0 C rup=rmaxis C rlw=rmaxis @@ -2464,8 +2505,10 @@ c bmmx=-1.0d+30 bmmn=1.0d+30 - call bfield(rcon(1),zcon(1),bphi,brr,bzz) call tor_curr(rcon(1),zcon(1),ajphi0) + call bpol(rcon(1),zcon(1),brr,bzz) + call equinum_fpol(0) + bphi=fpolv/rcon(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 @@ -2493,9 +2536,10 @@ c compute line integral on the contour psi=height^2 rpsim=rcon(inc1) zpsim=zcon(inc1) - call bfield(rpsim,zpsim,bphi,brr,bzz) + call bpol(rpsim,zpsim,brr,bzz) call tor_curr(rpsim,zpsim,ajphi) - + call equinum_fpol(0) + bphi=fpolv/rpsim btot=sqrt(bphi**2+brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2) dlpv(inc)=dlp @@ -2640,7 +2684,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ . ,vratja(jp),vratjb(jp) end do - rarea=sqrt(varea(nintp)/pi) +c rarea=sqrt(varea(nintp)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd @@ -2720,6 +2764,505 @@ c if(ip.eq.nintp) ip=nintp-1 fdvdrhot=spli(cdvdrhot,nintp,ip,dps) return end + + subroutine flux_average_an + implicit real*8 (a-h,o-z) + real*8 lam + + parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) + parameter(zero=0.0d0,one=1.0d0) + parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) + parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) + parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ + . njest*nnintp+nlest+54) + parameter(kwrk=nnintp+nlam+njest+nlest+3) + parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) + + dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) + dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) + dimension dadrhotv(nnintp),cdadrhot(nnintp,4) + dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4) + dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) + dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) + dimension rcon(2*ncnt+1),zcon(2*ncnt+1) + dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) + dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4) + dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4) + dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4) + dimension cfc(nnintp,4),crhotq(nnintp,4) + dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp) + dimension alam(nlam),fhlam(nnintp,nlam) + dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) + dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)) + dimension iwrk(kwrk),wrk(lwrk) + dimension ch01(lw01),weights(nlam) +c + common/cent/btrcen,rcen + common/parban/b0,rr0m,zr0m,rpam + common/parqq/q0,qa,alq + common/eqnn/nr,nz,npp,nintp + common/fpol/fpolv,ffpv +c common/psival/psin + + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz +c + common/pstab/rpstab + common/flav/vvol,rri,rbav,bmxpsi,bmnpsi + common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc + common/cratj/cratja,cratjb,cratjpl + common/crhotq/crhotq + common/phitedge/phitedge + common/cdadrhot/cdadrhot + common/cdvdrhot/cdvdrhot +c + common/coffvl/ch + common/coffdvl/ch01 + common/coffvlt/tjp,tlm + common/coffvln/njpt,nlmt + +c computation of flux surface averaged quantities + + rmaxis=rr0m + zmaxis=zr0m + btaxis=b0 + + call rhopol_an + phitedge=pi*rpam*rpam*btaxis + + write(71,*)' #i psin R z' + + nintp=nnintp + ninpr=(nintp-1)/10 + + dlam=1.0d0/dble(nlam-1) + do l=1,nlam-1 + alam(l)=dble(l-1)*dlam + fhlam(1,l)=sqrt(1.0d0-alam(l)) + ffhlam(l)=fhlam(1,l) + dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l)) + weights(l)=1.0d0 + end do + weights(1)=0.5d0 + weights(nlam)=0.5d0 + alam(nlam)=1.0d0 + fhlam(1,nlam)=0.0d0 + ffhlam(nlam)=0.0d0 + dffhlam(nlam)=-99999.0d0 + + jp=1 + anorm=2.0d0*pi*rmaxis/abs(btaxis) + b2av=btaxis**2 + dvdpsi=2.0d0*pi*anorm + dadpsi=2.0d0*pi/abs(btaxis) + ratio_cdator=abs(btaxis/btrcen) + ratio_cdbtor=1.0d0 + ratio_pltor=1.0d0 + qq=1.0d0 + fc=1.0d0 + + pstab(1)=0.0d0 + rpstab(1)=0.0d0 + rhotqv(1)=0.0d0 + vcurrp(1)=0.0d0 + vajphiav(1)=0.0d0 + bmxpsi(1)=abs(btaxis) + bmnpsi(1)=abs(btaxis) + bav(1)=abs(btaxis) + rbav(1)=1.0d0 + rri(1)=rmaxis + varea(1)=0.0d0 + vvol(1)=0.0d0 + vratjpl(1)=ratio_pltor + vratja(1)=ratio_cdator + vratjb(1)=ratio_cdbtor + ffc(1)=fc + rhot2q=0.0d0 + dadrhotv(1)=0.0d0 + dvdrhotv(1)=0.0d0 + qqv(1)=1.0d0 + +C rup=rmaxis +C rlw=rmaxis +C zup=(zbmax+zmaxis)/2.0d0 +C zlw=(zmaxis+zbmin)/2.0d0 + + do jp=2,nintp + height=dble(jp-1)/dble(nintp-1) + if(jp.eq.nintp) height=0.9999d0 + height2=height*height + ipr=0 + jpr=mod(jp,ninpr) + if(jpr.eq.1) ipr=1 + call contours_psi_an(height2,ncnt,rcon,zcon,ipr) +c + r2iav=0.0d0 + anorm=0.0d0 + dadpsi=0.0d0 + currp=0.0d0 + b2av=0.0d0 + area=0.0d0 + volume=0.0d0 + ajphiav=0.0d0 + bbav=0.0d0 + bmmx=-1.0d+30 + bmmn=1.0d+30 + + call equian(rcon(1),zcon(1)) + dbvcdc13=-ddpsidzz/rcon(1) + dbvcdc31= ddpsidrr/rcon(1)-bzz/rcon(1) + ajphi=ccj*(dbvcdc13-dbvcdc31) + brr=-dpsidz/rcon(1) + bzz= dpsidr/rcon(1) + bphi=fpolv/rcon(1) + btot0=sqrt(bphi**2+brr**2+bzz**2) + bpoloid0=sqrt(brr**2+bzz**2) + bv(1)=btot0 + bpv(1)=bpoloid0 + rpsim0=rcon(1) + + do inc=1,ncntt-1 + inc1=inc+1 + dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2) + dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2) + dlp=sqrt((rcon(inc1)-rcon(inc))**2+ + . (zcon(inc1)-zcon(inc))**2) + drc=(rcon(inc1)-rcon(inc)) + +c compute length, area and volume defined by psi=height^2 + + ph=0.5d0*(dla+dlb+dlp) + area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) + area=area+sqrt(area2) + rzp=rcon(inc1)*zcon(inc1) + rz=rcon(inc)*zcon(inc) + volume=pi*(rzp+rz)*drc+volume + +c compute line integral on the contour psi=height^2 + + rpsim=rcon(inc1) + zpsim=zcon(inc1) + call equian(rpsim,zpsim) + 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.5d0*dlp + anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0) + dadpsi=dadpsi+dlph* + . (1.0d0/(bpoloid*rpsim)+1.0d0/(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.0d0/(bpoloid*rpsim**2)+1.0d0/(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.0d0*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.0d0*pi*pi)) + qqv(jp)=qq + + rn=frhotor_an(sqrt(pstab(jp))) + qqan=q0+(qa-q0)*rn**alq + + dadr=2.0d0*pi*rn*rpam*rpam + dvdr=4.0d0*pi*pi*rn*rmaxis*rpam*rpam + +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(nintp-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.0d0 + shlam=0.0d0 + do l=nlam,1,-1 + lam=alam(l) + srl=0.0d0 + rl2=1.0d0-lam*bv(1)/bmmx + rl0=0.d0 + if(rl2.gt.0) rl0=sqrt(rl2) + do inc=1,ncntt-1 + rl2=1.0d0-lam*bv(inc+1)/bmmx + rl=0.0d0 + if(rl2.gt.0) rl=sqrt(rl2) + srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) + rl0=rl + end do + srl=srl/anorm + dhlam=0.5d0/srl + fc=fc+lam/srl*weights(l) + if(l.eq.nlam) then + fhlam(jp,l)=0.0d0 + ffhlam(nlam*(jp-1)+l)=0.0d0 + dffhlam(nlam*(jp-1)+l)=-dhlam + dhlam0=dhlam + else + shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam + fhlam(jp,l)=shlam + dffhlam(nlam*(jp-1)+l)=-dhlam + dhlam0=dhlam + end if + end do + fc=0.75d0*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,nintp + rhotqv(jp)=rhotqv(jp)/rhotqv(nintp) + if(jp.eq.nintp) then + rhotqv(jp)=1.0d0 + rpstab(jp)=1.0d0 + pstab(jp)=1.0d0 + end if + rhot_eq=frhotor_an(sqrt(pstab(jp))) + write(56,99) pstab(jp),rhot_eq,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(nintp)/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,nintp,iopt,cvol,ier) + iopt=0 + call difcs(rpstab,rbav,nintp,iopt,crbav,ier) + iopt=0 + call difcs(rpstab,rri,nintp,iopt,crri,ier) + iopt=0 + call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier) + iopt=0 + call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier) + iopt=0 + call difcs(rpstab,vratja,nintp,iopt,cratja,ier) + iopt=0 + call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier) + iopt=0 + call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier) + iopt=0 + call difcs(rpstab,varea,nintp,iopt,carea,ier) + iopt=0 + call difcs(rpstab,ffc,nintp,iopt,cfc,ier) + iopt=0 + call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier) + iopt=0 + call difcs(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier) + iopt=0 + call difcs(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier) + +c spline interpolation of H(lambda,rhop) and dH/dlambda + + iopt=0 + s=0.0d0 + call regrid(iopt,nintp,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 + + call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, + . ch01,lw01,ier) + + + return + 99 format(20(1x,e12.5)) + end + + + subroutine rhopol_an + implicit real*8(a-h,o-z) + parameter(nnr=101,nrest=nnr+4) + parameter(lwrk=nnr*4+nrest*16) + dimension rhop(nnr),rhot(nnr),rhopi(nnr),rhoti(nnr) + dimension psin(nnr) + dimension trp(nrest),crp(nrest) + dimension trot(nrest),crot(nrest) + dimension wrk(lwrk),iwrk(nrest),wp(nrest) + common/coffrtp/trp + common/coffrn/nsrp + common/coffrp/crp + common/coffrptt/trot + common/coffrnt/nsrot + common/coffrpt/crot + common/parqq/q0,qa,alq + common/parban/b0,rr0m,zr0m,rpam + common/pareq1/psia + common/sgnib/sgnbphi,sgniphi + + + dr=1.0d0/dble(nnr-1) + rhot(1)=0.0d0 + psin(1)=0.0d0 + res=0.0d0 + fq0=0.0d0 + do i=2,nnr + rhot(i)=(i-1)*dr + rn=rhot(i) + qq=q0+(qa-q0)*rn**alq + fq1=rn/qq + res=res+0.5d0*(fq1+fq0)*dr + fq0=fq1 + psin(i)=res + end do + + wp=1.0d0 + psin=psin/res + rhop=sqrt(psin) + psia=-sgniphi*abs(res*rpam*rpam*b0) + print*,psia,log(8.0d0*rr0m/rpam)-2.0d0 + wp(1)=1.0d3 + wp(nnr)=1.0d3 + +c spline interpolation of rhopol versus rhotor + iopt=0 + xb=0.0d0 + xe=1.0d0 + ss=0.00001d0 + kspl=3 + call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp, + . trp,crp,rp,wrk,lwrk,iwrk,ier) + call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) + +c spline interpolation of rhotor versus rhopol + iopt=0 + xb=0.0d0 + xe=1.0d0 + ss=0.00001d0 + kspl=3 + call curfit(iopt,nnr,rhop,rhot,wp,xb,xe,kspl,ss,nrest,nsrot, + . trot,crot,rp,wrk,lwrk,iwrk,ier) + call splev(trot,nsrot,crot,3,rhop,rhoti,nnr,ier) + + do i=1,nnr + write(54,*) rhop(i),rhot(i),rhopi(i),rhoti(i) + end do + + return + end + + + function frhotor_an(rhop) + implicit real*8(a-h,o-z) + parameter(nnr=101,nrest=nnr+4) + dimension trot(nrest),crot(nrest),rrs(1),ffspl(1) + common/coffrptt/trot + common/coffrnt/nsrot + common/coffrpt/crot + rrs(1)=rhop + call splev(trot,nsrot,crot,3,rrs,ffspl,1,ier) + frhotor_an=ffspl(1) + return + end + + + subroutine contours_psi_an(h,np,rcn,zcn,ipr) + implicit real*8 (a-h,o-z) + parameter(pi=3.14159265358979d0) + parameter(mest=4,kspl=3) + dimension rcn(2*np+1),zcn(2*np+1) + common/parban/b0,rr0m,zr0m,rpam + + th=pi/dble(np) + rn=frhotor_an(sqrt(h)) + do ic=1,2*np+1 + zcn(ic)=zr0m+rpam*rn*sin(th*(ic-1)) + rcn(ic)=rr0m+rpam*rn*cos(th*(ic-1)) + + if (ipr.gt.0) then + write(71,111) ic,h,rcn(ic),zcn(ic) + end if + end do + write(71,*) + +111 format(i6,12(1x,e12.5)) + return + end + + + subroutine vectinit implicit real*8 (a-h,o-z) @@ -3418,10 +3961,12 @@ c common/ygyg/yg common/dxgyg/derxg,deryg common/iieq/iequil - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz + common/fpol/fpolv,ffpv common/psival/psinv common/sgnib/sgnbphi,sgniphi -c + xg=0.0d0 yg=9.9d1 c @@ -3458,17 +4003,19 @@ c c if(iequil.eq.1) then call equian(rrm,zzm) - yg=fpolv/(rrm*bres) end if c if(iequil.eq.2) then - call equinum(rrm,zzm) + call equinum_psi(rrm,zzm) + call equinum_derpsi(rrm,zzm,3) + call equinum_fpol(1) + end if + call sub_xg_derxg yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) if(psinv.lt.0.0d0) return - end if c c B = f(psi)/R e_phi+ grad psi x e_phi/R c @@ -3573,23 +4120,20 @@ c common/parban/b0,rr0m,zr0m,rpam common/psival/psinv common/pareq1/psia - common/densbnd/psdbnd - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz + common/fpol/fpolv,ffpv common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddens common/sgnib/sgnbphi,sgniphi - common/bmxmn/bmxi,bmni - common/fc/fci - common/iipr/iprof -c - if(iprof.eq.0) psdbnd=1.0d0 -c + dpsidrp=0.0d0 d2psidrp=0.0d0 c -c simple model for equilibrium +c simple model for equilibrium: large aspect ratio +c outside plasma: analytical continuation, not solution Maxwell eqs c rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2) rn=rpm/rpam @@ -3601,9 +4145,14 @@ c cst=(rrm-rr0m)/rpm end if c -c valore fittizio di psinv e di psia: - psinv=rn**2 - psia=sgniphi + rhot=rn + if(rn.le.1) then + rhop=frhopol(rhot) + psinv=rhop*rhop + else + psinv=1.0d0+B0*rpam**2/2.0d0/abs(psia)/qa*(rn*rn-1.0d0) + rhop=sqrt(psinv) + end if c sgn=sgniphi*sgnbphi if(rn.le.1.0d0) then @@ -3612,39 +4161,25 @@ c dqq=alq*(qa-q0)*rn**(alq-1.0d0) d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn else - dpsidrp=B0*rpam/qa/rn*sgn - d2psidrp=-B0/qa/rn**2 *sgn + dpsidrp=B0*rpam/qa*rn*sgn + d2psidrp=B0/qa*sgn end if c fpolv=sgnbphi*b0*rr0m dfpolv=0.0d0 - ffpv=sgniphi*fpolv*dfpolv + ffpv=0.0d0 c dpsidr=dpsidrp*cst dpsidz=dpsidrp*snt ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp -c - dens=0.0d0 - xg=0.0d0 - dxgdpsi=0.0d0 - if(psinv.lt.psdbnd) then - call density(psinv) - end if - xg=xgcn*dens - ddensdrp=2.0d0*rn*ddens - if(dpsidrp.ne.0.0d0) dxgdpsi=xgcn*ddensdrp/(dpsidrp*rpam) -c - bmax=sqrt(fpolv**2+dpsidrp**2)/(rr0m-rpam*rn) - bmin=sqrt(fpolv**2+dpsidrp**2)/(rr0m+rpam*rn) -c + return end -c -c -c - subroutine equinum(rpsim,zpsim) + + + subroutine equinum_psi(rpsim,zpsim) implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) @@ -3652,40 +4187,19 @@ c dimension ccspl(nnw*nnh) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) dimension rrs(1),zzs(1),ffspl(1) - parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) - parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) - parameter(lw11=nnw*3+nnh*3+nnw*nnh) parameter(nrs=1,nzs=1) - dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) - dimension tfp(nrest),cfp(nrest),wrkfd(nrest) c - common/eqnn/nr,nz,npp,nintp common/psival/psinv - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv - common/pareq1/psia common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/pareq1t/psiant,psinop c common/coffeqt/tr,tz - common/coffeqtp/tfp common/coffeq/ccspl - common/coffeqd/cc01,cc10,cc20,cc02,cc11 common/coffeqn/nsrt,nszt,nsft - common/cofffp/cfp - common/fpas/fpolas c psinv=-1.0d0 -c - dpsidr=0.0d0 - dpsidz=0.0d0 - ddpsidrr=0.0d0 - ddpsidzz=0.0d0 - ddpsidrz=0.0d0 c c here lengths are measured in meters -c - fpolv=fpolas - ffpv=0.0d0 c if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return @@ -3697,84 +4211,181 @@ c call fpbisp(tr,nsr,tz,nsz,ccspl,3,3, . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) psinv=(ffspl(1)-psinop)/psiant - if(psinv.lt.0.0d0) - . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim c + return + end + + subroutine equinum_derpsi(rpsim,zpsim,iderpsi) + implicit real*8 (a-h,o-z) + parameter(nnw=501,nnh=501) + parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) + parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) + parameter(lw11=nnw*3+nnh*3+nnw*nnh) + dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) + integer*4 iderpsi +c + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz + common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz + common/coffeqd/cc01,cc10,cc20,cc02,cc11 +c + dpsidr=0.0d0 + dpsidz=0.0d0 + ddpsidrr=0.0d0 + ddpsidzz=0.0d0 + ddpsidrz=0.0d0 + +c here lengths are measured in meters + if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return + if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return + + select case(iderpsi) + + case(1) nur=1 nuz=0 - kkr=3-nur - kkz=3-nuz - iwr=1+(nr-nur-4)*(nz-nuz-4) - iwz=iwr+4-nur - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz, - . rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2)) - dpsidr= ffspl(1)*psia -c + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc10,lw10) + dpsidr=derpsi nur=0 nuz=1 - kkr=3-nur - kkz=3-nuz - iwr=1+(nr-nur-4)*(nz-nuz-4) - iwz=iwr+4-nur - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz, - . rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2)) - dpsidz= ffspl(1)*psia -c + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc01,lw01) + dpsidz=derpsi + + case(2) nur=2 nuz=0 - kkr=3-nur - kkz=3-nuz - iwr=1+(nr-nur-4)*(nz-nuz-4) - iwz=iwr+4-nur - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz, - . rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2)) - ddpsidrr= ffspl(1)*psia -c + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc20,lw20) + ddpsidrr=derpsi nur=0 nuz=2 - kkr=3-nur - kkz=3-nuz - iwr=1+(nr-nur-4)*(nz-nuz-4) - iwz=iwr+4-nur - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz, - . rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2)) - ddpsidzz= ffspl(1)*psia -c + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc02,lw02) + ddpsidzz=derpsi nur=1 nuz=1 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc11,lw11) + ddpsidrz=derpsi + + case(3) + nur=1 + nuz=0 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc10,lw10) + dpsidr=derpsi + nur=0 + nuz=1 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc01,lw01) + dpsidz=derpsi + nur=2 + nuz=0 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc20,lw20) + ddpsidrr=derpsi + nur=0 + nuz=2 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc02,lw02) + ddpsidzz=derpsi + nur=1 + nuz=1 + call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc11,lw11) + ddpsidrz=derpsi + + case default + print*,'iderpsi undefined' + + end select + + return + end + + subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw) + implicit real*8 (a-h,o-z) + parameter(nnw=501,nnh=501) + parameter(nrest=nnw+4,nzest=nnh+4) + parameter(liwrk=2) + parameter(nrs=1,nzs=1) + dimension tr(nrest),tz(nzest),iwrk(liwrk) + dimension rrs(1),zzs(1),ffspl(1) + dimension cc(lw) + + common/pareq1/psia + common/coffeqt/tr,tz + common/coffeqn/nsrt,nszt,nsft + common/eqnn/nr,nz,npp,nintp + + rrs(1)=rpsim + zzs(1)=zpsim + nsr=nsrt + nsz=nszt kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz, - . rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2)) - ddpsidrz= ffspl(1)*psia + call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kkr,kkz + . ,rrs,nrs,zzs,nzs,ffspl,cc(iwr),cc(iwz),iwrk(1),iwrk(2)) + derpsi=ffspl(1)*psia + + return + end + + subroutine equinum_fpol(iderfpol) + implicit real*8 (a-h,o-z) + parameter(nnw=501) + parameter(nrest=nnw+4) + dimension rrs(1),ffspl(1) + dimension tfp(nrest),cfp(nrest),wrkfd(nrest) + integer*4 iderfpol c + common/psival/psinv + common/fpol/fpolv,ffpv + common/pareq1/psia + common/coffeqtp/tfp + common/coffeqn/nsrt,nszt,nsft + common/cofffp/cfp + common/fpas/fpolas + + fpolv=fpolas + dfpolv=0.0d0 + ffpv=0.0d0 + if(iderfpol.lt.0.or.iderfpol.gt.1) return + if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then rrs(1)=psinv call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier) fpolv=ffspl(1) - nu=1 - call splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier) - dfpolv=ffspl(1) - ffpv=fpolv*dfpolv/psia + if(iderfpol.eq.1) then + call splder(tfp,nsft,cfp,3,1,rrs,ffspl,1,wrkfd,ier) + dfpolv=ffspl(1) + ffpv=fpolv*dfpolv/psia + end if end if -c + return end -c -c -c + subroutine bfield(rpsim,zpsim,bphi,brr,bzz) implicit real*8 (a-h,o-z) - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv - call equinum(rpsim,zpsim) + call btor(rpsim,zpsim,bphi) + call bpol(rpsim,zpsim,brr,bzz) + return + end + + subroutine btor(rpsim,zpsim,bphi) + implicit real*8 (a-h,o-z) + common/psival/psinv + common/fpol/fpolv,ffpv + call equinum_psi(rpsim,zpsim) + call equinum_fpol(0) bphi=fpolv/rpsim + return + end + + subroutine bpol(rpsim,zpsim,brr,bzz) + implicit real*8 (a-h,o-z) + common/derip1/dpsidr,dpsidz + call equinum_derpsi(rpsim,zpsim,1) brr=-dpsidz/rpsim bzz= dpsidr/rpsim return end -c + subroutine tor_curr_psi(h,ajphi) implicit real*8 (a-h,o-z) common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz @@ -3782,13 +4393,13 @@ c call tor_curr(r2,zmaxis,ajphi) return end - c subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : pi,ccj=>mu0inv implicit real*8 (a-h,o-z) - common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv - call equinum(rpsim,zpsim) + common/derip1/dpsidr,dpsidz + common/derip2/ddpsidrr,ddpsidzz,ddpsidrz + call equinum_derpsi(rpsim,zpsim,3) bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim @@ -4277,8 +4888,10 @@ c use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree, * ui=>im implicit real*8 (a-h,o-z) + complex*16 ui parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36) + parameter(ui=(0.0d0,1.0d0)) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ytmp(ndim),yptmp(ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) @@ -4478,7 +5091,6 @@ c dimension ytmp(ndim),yptmp(ndim) dimension yyrfl(jmx,kmx,ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) - dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) @@ -4718,18 +5330,12 @@ c rbavi=1.0d0 rrii=rr0m rhop0=sqrt(psjki(j,k,i-1)) - if(iequil.eq.2) then intp=1 psinv=psjki(j,k,i) rhop=sqrt(psinv) call valpsispl(rhop,voli,dervoli,area,rrii, . rbavi,bmxi,bmni,fci,intp) dvvl=dervoli*abs(rhop-rhop0) - else - dvvl=4.0d0*pi*rr0m*pi*abs(rhop*rhop-rhop0*rhop0)*rpam**2 - rbavi=b0/bmni - fci=1.0d0-1.469d0*sqrt(rpam/rr0m) - end if if(dvvl.le.0.0d0) dvvl=1.0d0 c adnm=2.0d0*pi*rrii @@ -4863,13 +5469,13 @@ c end if return -999 format(12(1x,e12.5)) end c c c subroutine dispersion(lrm) implicit real*8(a-h,o-z) + parameter(one=1.0d0) complex*16 cc0,cc2,cc4,rr complex*16 anpr2a,anpra complex*16 anpr,anpr2,ex,ey,ez,den @@ -4885,7 +5491,6 @@ c common/warm/iwarm,ilarm common/nprw/anprr,anpri common/epolar/ex,ey,ez - common/amut/amu common/imx/imx errnpr=1.0d0 @@ -5277,7 +5882,7 @@ c subroutine hermitian_2(rr,lrm) implicit real*8(a-h,o-z) - parameter(tmax=5,npts=500) + parameter(tmax=5) dimension rr(-lrm:lrm,0:2,0:lrm) parameter(epsa=0.0d0,epsr=1.0d-4) parameter (lw=5000,liw=lw/4) @@ -5288,7 +5893,6 @@ c common/amut/amu common/nplr/anpl,anpr common/warm/iwarm,ilarm - common/cri/cr,ci common/nmhermit/n,m,ih c do n=-lrm,lrm @@ -5302,10 +5906,6 @@ c llm=min(3,lrm) c bth2=2.0d0/amu - bth=sqrt(bth2) - amu2=amu*amu - amu4=amu2*amu2 - amu6=amu4*amu2 c n1=1 if(iwarm.gt.10) n1=-llm @@ -5438,7 +6038,6 @@ c upl=bth*x gx=1.0d0+t*t/amu exdxdt=cr*exp(-t*t)*gx/rxt - nn=abs(n) gr=anpl*upl+n*yg zm=-amu*(gx-gr) s=amu*(gx+gr) @@ -6401,17 +7000,17 @@ c yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 rti=sqrt(xti**2+yti**2) -c - dirxt= (dirx*csps1-diry*snps1)/dir - diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir - dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir -c - if(k.eq.1) then - xti1=xti - yti1=yti - zti1=zti - rti1=rti - end if +Cc +C dirxt= (dirx*csps1-diry*snps1)/dir +C diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir +C dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir +Cc +C if(k.eq.1) then +C xti1=xti +C yti1=yti +C zti1=zti +C rti1=rti +C end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 @@ -6467,12 +7066,10 @@ c common/pcjki/ppabs,ccci common/dpjjki/pdjki,currj c - common/cent/btrcen,rcen common/angles/alpha0,beta0 common/iieq/iequil common/parban/b0,rr0m,zr0m,rpam common/taumnx/taumn,taumx,pabstot,currtot - common/jmxmn/rhot1,rhot2,aj1,aj2 c common/polcof/psipol,chipol c @@ -6507,7 +7104,6 @@ c psit=rt**2 psit1=rt1**2 end if - if(iequil.eq.2) then rhotv(it)=frhotor(psit) call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii, . rbavi,bmxi,bmni,fci,intp) @@ -6519,15 +7115,6 @@ c ratjav(it)=ratjai ratjbv(it)=ratjbi ratjplv(it)=ratjpli - else - rhotv(it)=sqrt(psit) - area1=pi*psit1*rpam**2 - voli1=2.0d0*pi*rr0m*area1 - dvolt(it)=abs(voli1-voli0) - darea(it)=abs(areai1-areai0) - voli0=voli1 - areai0=areai1 - end if end do kkk=1 @@ -6790,7 +7377,6 @@ c dPdV [MW/m^3], Jcd [MA/m^2] end do return - 49 format(i5,20(1x,e12.5)) 99 format(30(1x,e12.5)) end c @@ -6874,17 +7460,11 @@ c end if end if c - if(iequil.eq.2) then rhotmx=frhotor(rhopmx**2) rhotmn=frhotor(rhopmn**2) rhote2=frhotor(psie2) rhote1=frhotor(psie1) drhot=rhote2-rhote1 - else - rhotmx=rhopmx - rhotmn=rhopmn - drhot=drhop - end if c if(ipk.eq.2) then drhop=-drhop @@ -6966,7 +7546,7 @@ c real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2 real*8 deno,denx,anpl2,dnl,del0 - real*8 beta0,alpha0,gam + real*8 gam real*8 sox complex*16 exom,eyom,exxm,eyxm,ext,eyt c @@ -7017,7 +7597,6 @@ c gam=atan(sngam/csgam)*180.d0/pi return - 111 format(20(1x,e12.5)) end subroutine stokes(ext,eyt,qq,uu,vv) @@ -7056,7 +7635,7 @@ c ell=bb/aa if(iequil.eq.1) then call equian(rrm,zzm) else - call equinum(rrm,zzm) + call equinum_psi(rrm,zzm) end if if (psinv.ge.0.0d0.and.psinv.lt.psdbnd) then @@ -7151,17 +7730,16 @@ c wave vector and electric field after reflection in lab frame eztr=dot_product(vv3,evrfl) return - 111 format(20(1x,e12.5)) end subroutine vacuum_rt(xvstart,anv,xvend,ivac) use reflections, only : inters_linewall,inside implicit none integer*4 ivac - integer nbb,nlim,i,imax + integer nbb,nlim,i parameter(nbb=5000) - real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax - real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6) + real*8 st,rrm,zzm,dst,dstvac,smax + real*8 anv(3),xvstart(3),xvend(3),walln(3) real*8 xv0(3),anv0(3) real*8 rlim(nbb),zlim(nbb) logical plfound,inside_plasma @@ -7169,8 +7747,6 @@ c wave vector and electric field after reflection in lab frame common/limiter/rlim,zlim,nlim common/dsds/dst - common/psival/psinv - common/densbnd/psdbnd common/dstvac/dstvac c ivac=1 plasma hit before wall reflection @@ -7228,6 +7804,5 @@ c ivac=-1 vessel (and thus plasma) never crossed anv=anv0 return - 111 format(20(1x,e12.5)) end