diff --git a/src/gray.f b/src/gray.f index b976a81..93f3ef9 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1065,12 +1065,14 @@ c parameter(nnw=501,nnh=501) parameter(pi=3.14159265358979d0) parameter(nbb=1000) +c parameter(np=100) character*48 stringa dimension fpol(nnw),pres(nnw),qpsi(nnw) dimension ffprim(nnw),pprim(nnw) dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) dimension rbbbs(nbb),zbbbs(nbb) dimension rlim(nbb),zlim(nbb) +c dimension rcon(2*np+1),zcon(2*np+1) c parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54) @@ -1089,7 +1091,7 @@ c dimension fpoli(nnw) c common/pareq1/psia - common/pareq1a/psiaxis0 + common/pareq1t/psiant,psinop common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/psinr/psinr @@ -1130,9 +1132,9 @@ c end if if(ipsinorm.eq.0) then read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid - read (neqdsk,2020) rmaxis,zmaxis,psiax,psiedge,btrcen - read (neqdsk,2020) current,simag,xdum,rmaxis,xdum - read (neqdsk,2020) zmaxis,xdum,sibry,xdum,xdum + read (neqdsk,2020) rmaxis,zmaxis,psiaxis,psiedge,btrcen + read (neqdsk,2020) current,xdum,xdum,xdum,xdum + read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum read (neqdsk,2020) (fpol(i),i=1,nr) read (neqdsk,2020) (pres(i),i=1,nr) read (neqdsk,2020) (ffprim(i),i=1,nr) @@ -1141,9 +1143,9 @@ c read (neqdsk,2020) (qpsi(i),i=1,nr) else read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid - read (neqdsk,*) rmaxis,zmaxis,psiax,psiedge,btrcen - read (neqdsk,*) current,simag,xdum,rmaxis,xdum - read (neqdsk,*) zmaxis,xdum,sibry,xdum,xdum + read (neqdsk,*) rmaxis,zmaxis,psiaxis,psiedge,btrcen + read (neqdsk,*) current,xdum,xdum,xdum,xdum + read (neqdsk,*) xdum,xdum,xdum,xdum,xdum read (neqdsk,*) (fpol(i),i=1,nr) read (neqdsk,*) (pres(i),i=1,nr) read (neqdsk,*) (ffprim(i),i=1,nr) @@ -1172,7 +1174,7 @@ c c icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip c icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip psiedge=-psiedge - psiax=-psiax + psiaxis=-psiaxis if (ipsinorm.eq.0) then do j=1,nz do i=1,nr @@ -1213,7 +1215,7 @@ c c psi function - psia0=psiedge-psiax + psia0=psiedge-psiaxis c icocos=0: adapt psi to force Ip sign, otherwise maintain psi if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0) current=sign(current,sgniphi) @@ -1222,12 +1224,18 @@ c icocos=0: adapt psi to force Ip sign, otherwise maintain psi c icocos>10: input psi is in Wb c icocos<10: input psi is in Wb/rad (gray convention) if (icocos.ge.10) psia=psia/(2.0d0*pi) +c +c do j=1,nz +c do i=1,nr +c write(80,2021) rv(i),zv(j),psi(i,j) +c enddo +c write(80,*) ' ' +c enddo - psiaxis0=0.0d0 do j=1,nz do i=1,nr if(ipsinorm.eq.0) then - psin(i,j)=(psi(i,j)-psiax)/(psia0) + psin(i,j)=(psi(i,j)-psiaxis)/(psia0) psi(i,j)=psin(i,j)*psia else psi(i,j)=psin(i,j)*psia @@ -1237,12 +1245,13 @@ c icocos<10: input psi is in Wb/rad (gray convention) enddo enddo c -c spline interpolation of psi(i,j) and derivatives +c spline fitting/interpolation of psin(i,j) and derivatives c iopt=0 call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, . wrk,lwrk,iwrk,liwrk,ier) +c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then sspl=0.0d0 call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, @@ -1291,6 +1300,32 @@ c nuz=1 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier) +c +c scaling of f_poloidal +c +c icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol + if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr)) + btrcen=sign(btrcen,sgnbphi) + + do i=1,nr + fpol(i)=sgnbphi*abs(fpol(i))*factb + wf(i)=1.0d0 + end do + wf(nr)=1.0d2 +c +c spline interpolation of fpol(i) +c + iopt=0 + xb=0.0d0 + xe=1.0d0 + ssfp=0.01d0 + call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft, + . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier) +c + call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier) + fpolas=fpoli(nr) +c +c read plasma boundaries from eqdsk file c read (neqdsk,*) nbbbs,nlim if(nbbbs.gt.0) then @@ -1298,6 +1333,9 @@ c . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) if(ipsinorm.eq.0) . read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) +c do i=1,nbbbs +c write(51,*) rbbbs(i),zbbbs(i) +c end do end if if(nlim.gt.0) then if(ipsinorm.eq.1) @@ -1332,37 +1370,18 @@ c if(zbmax.ge.zmxm) zbmax=zbmax-dz if(rbmax.ge.rmxm) rbmax=rbmax-dr c -c scaling of f_poloidal +c start with uncorrected normalized psi c -c icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol - if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr)) - btrcen=sign(btrcen,sgnbphi) - - do i=1,nr - fpol(i)=sgnbphi*abs(fpol(i))*factb - wf(i)=1.0d0 - end do - wf(nr)=1.0d2 - -c -c spline interpolation of fpol(i) -c - iopt=0 - xb=0.0d0 - xe=1.0d0 - ssfp=0.01d0 - call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft, - . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier) -c - call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier) - fpolas=fpoli(nr) + psinop=0.0d0 + psinxp=1.0d0 + psiant=1.0d0 c c search for O-point c - call points_ox(rmaxis,zmaxis,rmop,zmop,psinop,info) + call points_ox(rmaxis,zmaxis,rmop,zmop,psinoptmp,info) rmaxis=rmop zmaxis=zmop - print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinop + print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinoptmp c c search for X-point if ixp not = 0 c @@ -1370,15 +1389,14 @@ c if(ixp.lt.0) then r10=rbmin z10=zbmin - call points_ox(r10,z10,rxp,zxp,psinxp,info) + call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then - print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxp + print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxptmp rbmin=rxp zbmin=zxp - delpsinox=(psinxp-psinop) - psia=psia*delpsinox - deltapsi=abs(psia) - psiaxis0=psia*psinop + psinop=psinoptmp + psinxp=psinxptmp + psiant=psinxp-psinop psin1=1.0d0 r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 @@ -1392,14 +1410,13 @@ c print'(a)','no X-point' else r10=rmop z10=zbmax - call points_ox(r10,z10,rxp,zxp,psinxp,info) + call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then - print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxp + print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxptmp zbmax=zxp - delpsinox=(psinxp-psinop) - psia=psia*delpsinox - deltapsi=abs(psia) - psiaxis0=psia*psinop + psinop=psinoptmp + psinxp=psinxptmp + psiant=psinxp-psinop psin1=1.0d0 z10=(zbmin+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) @@ -1413,10 +1430,8 @@ c print'(a)','no X-point' c if (ixp.eq.0) then psin1=1.0d0 - delpsinox=(psin1-psinop) - psia=psia*delpsinox - deltapsi=abs(psia) - psiaxis0=psia*psinop + psinop=psinoptmp + psiant=psin1-psinop r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) @@ -1440,7 +1455,7 @@ c compute B_toroidal on axis c compute normalized rho_tor from eqdsk q profile call rhotor(nr) -c phitedge=deltapsi*rhotsx*2*pi +c phitedge=abs(psia)*rhotsx*2*pi c rrtor=sqrt(phitedge/abs(btrcen)/pi) c compute flux surface averaged quantities @@ -1450,6 +1465,12 @@ c compute flux surface averaged quantities zup=zmaxis+(zbmax-zmaxis)/10.0d0 zlw=zmaxis-(zmaxis-zbmin)/10.0d0 call flux_average +c ipr=1 +c call contours_psi(1.0d0,np,rcon,zcon,ipr) +c do ii=1,2*np+1 +c write(52,*) rcon(ii), zcon(ii) +c end do +c c locate psi surface for q=1.5 and q=2 @@ -2066,7 +2087,7 @@ c dimension czc(nrest),zeroc(mest) c common/pareq1/psia - common/pareq1a/psiaxis0 + common/pareq1t/psiant,psinop common/coffeqn/nsr,nsz,nsft common/coffeq/cc common/coffeqt/tr,tz @@ -2092,7 +2113,7 @@ c iopt=1 call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier - val=h+psiaxis0/psia + val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) @@ -3401,7 +3422,7 @@ c common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/pareq1/psia common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/pareq1a/psiaxis0 + common/pareq1t/psiant,psinop c common/coffeqt/tr,tz common/coffeqtp/tfp @@ -3433,9 +3454,9 @@ c nsz=nszt 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)-psiaxis0/psia -c if(psinv.lt.0.0d0) -c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim + psinv=(ffspl(1)-psinop)/psiant + if(psinv.lt.0.0d0) + . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim c nur=1 nuz=0 @@ -3541,8 +3562,7 @@ c dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c - common/pareq1/psia - common/pareq1a/psiaxis0 + common/pareq1t/psiant,psinop common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/coffeqn/nsr,nsz,nsft common/coffeq/cc @@ -3552,7 +3572,7 @@ c 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+psiaxis0/psia + val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2)