istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

renamed clashing common block names; fixed wrong calling of splining routines; fixed spline evaluation functions to work with uneven x grid spacing; removed debugging prints

Browse Source
This commit is contained in:
Lorenzo Figini 2014-06-19 08:13:52 +00:00
parent 074f331355
commit 2c90c5f2cf
4 changed files with 102 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -97,7 +97,7 @@ clean:
# Dependencies # Dependencies
# ------------ # ------------
gray_main.o: const_and_precisions.o gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
green_func_p.o: const_and_precisions.o green_func_p.o: const_and_precisions.o
scatterspl.o: const_and_precisions.o scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o beamdata.o: const_and_precisions.o

192
src/gray-externals.f
View File

@ -114,7 +114,7 @@ c
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/istgr/istpr,istpl common/istgr/istpr,istpl
c c
common/psival/psinv common/psinv/psinv
common/psinv11/psinv11 common/psinv11/psinv11
common/ierr/ierr common/ierr/ierr
common/taumnx/taumn,taumx,pabstot,currtot common/taumnx/taumn,taumx,pabstot,currtot
@ -134,7 +134,7 @@ c
common/index_rt/index_rt common/index_rt/index_rt
common/ipass/ipass common/ipass/ipass
common/rwallm/rwallm common/rwallm/rwallm
common/bound/zbmin,zbmax common/zbound/zbmin,zbmax
pabstot=0.0d0 pabstot=0.0d0
currtot=0.0d0 currtot=0.0d0
@ -306,7 +306,7 @@ c
common/btot/btot common/btot/btot
common/xgxg/xg common/xgxg/xg
common/ygyg/yg common/ygyg/yg
common/dens/dens,ddens common/dddens/dens,ddens
common/tete/tekev common/tete/tekev
common/absor/alpha,effjcd,akim,tau0 common/absor/alpha,effjcd,akim,tau0
common/densbnd/psdbnd common/densbnd/psdbnd
@ -732,9 +732,6 @@ c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
anx0c=(anr0c*x00-anphi0c*y00)/r00 anx0c=(anr0c*x00-anphi0c*y00)/r00
any0c=(anr0c*y00+anphi0c*x00)/r00 any0c=(anr0c*y00+anphi0c*x00)/r00
c
print*,' input file read'
! call myflush
c c
c read data for Te , ne , Zeff from file if iprof>0 c read data for Te , ne , Zeff from file if iprof>0
c c
@ -748,8 +745,6 @@ c read profiles from input arguments
call profdata(nrho, psijet, te, dne, zeff) call profdata(nrho, psijet, te, dne, zeff)
c close(nprof) c close(nprof)
end if end if
print*,' profiles fitted'
! call myflush
c c
c read equilibrium data from file if iequil=2 c read equilibrium data from file if iequil=2
c c
@ -762,8 +757,6 @@ c . status= 'unknown', unit=neqdsk)
call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd, call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd,
. rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf) . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf)
c close(neqdsk) c close(neqdsk)
print*,' equilibrium fitted'
! call myflush
c print density, temperature, safecty factor, toroidal current dens c print density, temperature, safecty factor, toroidal current dens
c versus psi, rhop, rhot c versus psi, rhop, rhot
@ -964,16 +957,16 @@ c
if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
call vlocate(alphastv,nisteer,alpha0,k) call vlocate(alphastv,nisteer,alpha0,k)
dal=alpha0-alphastv(k) dal=alpha0-alphastv(k)
betst=fspli(cbeta,nisteer,k,dal) betst=fspli(cbeta,nstrmx,k,dal)
x00=fspli(cx0,nisteer,k,dal) x00=fspli(cx0,nstrmx,k,dal)
y00=fspli(cy0,nisteer,k,dal) y00=fspli(cy0,nstrmx,k,dal)
z00=fspli(cz0,nisteer,k,dal) z00=fspli(cz0,nstrmx,k,dal)
wcsi=fspli(cwaist1,nisteer,k,dal) wcsi=fspli(cwaist1,nstrmx,k,dal)
weta=fspli(cwaist2,nisteer,k,dal) weta=fspli(cwaist2,nstrmx,k,dal)
rcicsi=fspli(crci1,nisteer,k,dal) rcicsi=fspli(crci1,nstrmx,k,dal)
rcieta=fspli(crci2,nisteer,k,dal) rcieta=fspli(crci2,nstrmx,k,dal)
phiw=fspli(cphi1,nisteer,k,dal) phiw=fspli(cphi1,nstrmx,k,dal)
phir=fspli(cphi2,nisteer,k,dal) phir=fspli(cphi2,nstrmx,k,dal)
else else
write(*,*) ' alpha0 outside table range !!!' write(*,*) ' alpha0 outside table range !!!'
if(alpha0.ge.alphastv(nisteer)) ii=nisteer if(alpha0.ge.alphastv(nisteer)) ii=nisteer
@ -1043,7 +1036,7 @@ c
common/ipsn/ipsinorm common/ipsn/ipsinorm
common/sspl/sspl common/sspl/sspl
common/nfile/neqdsk,nprof common/nfile/neqdsk,nprof
common/bound/zbmin,zbmax common/zbound/zbmin,zbmax
common/sgnib/sgnbphi,sgniphi common/sgnib/sgnbphi,sgniphi
common/factb/factb common/factb/factb
common/ixp/ixp common/ixp/ixp
@ -1106,7 +1099,7 @@ c psi function
psia0=psiedge-psiaxis psia0=psiedge-psiaxis
psia=psia0*factb psia=psia0*factb
sgniphi=sign(1.0d0,-psia0) sgniphi=sign(1.0d0,-psia0)
cc c
c do j=1,nz c do j=1,nz
c do i=1,nr c do i=1,nr
c write(620,2021) rv(i),zv(j),psin(i,j) c write(620,2021) rv(i),zv(j),psin(i,j)
@ -1172,21 +1165,21 @@ c
nsr=nr/4+4 nsr=nr/4+4
nsz=nz/4+4 nsz=nz/4+4
call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
. rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier)
c if ier=-1 data are fitted using sspl=0 c if ier=-1 data are fitted using sspl=0
if(ier.eq.-1) then if(ier.eq.-1) then
sspl=0.0d0 sspl=0.0d0
nsr=nr/4+4 nsr=nr/4+4
nsz=nz/4+4 nsz=nz/4+4
call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
. rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier)
end if end if
nsrt=nsr nsrt=nsr
nszt=nsz nszt=nsz
end if end if
cc c
cc re-evaluate psi on the grid using the spline (only for debug and cniteq) c re-evaluate psi on the grid using the spline (only for debug and cniteq)
cc c
c call dierckx_bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, c call dierckx_bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
c c
@ -1354,15 +1347,15 @@ c
c compute B_toroidal on axis c compute B_toroidal on axis
btaxis=fpol(1)/rmaxis btaxis=fpol(1)/rmaxis
btrcen=abs(btrcen)*factb btrcen=btrcen*factb
write(*,'(a,f8.4)') 'factb = ',factb write(*,'(a,f8.4)') 'factb = ',factb
write(*,'(a,f8.4)') '|BT_centr|= ',btrcen write(*,'(a,f8.4)') '|BT_centr|= ',abs(btrcen)
write(*,'(a,f8.4)') '|BT_axis| = ',abs(btaxis) write(*,'(a,f8.4)') '|BT_axis| = ',abs(btaxis)
c compute normalized rho_tor from eqdsk q profile c compute normalized rho_tor from eqdsk q profile
call rhotor(nrho) call rhotor(nrho)
phitedge=abs(psia)*rhotsx*2*pi phitedge=abs(psia)*rhotsx*2*pi
rrtor=sqrt(phitedge/abs(btrcen)/pi) rrtor=sqrt(abs(phitedge/btrcen)/pi)
call rhopol call rhopol
c write(*,*) rhotsx,phitedge,rrtor,abs(psia) c write(*,*) rhotsx,phitedge,rrtor,abs(psia)
@ -1417,7 +1410,7 @@ c
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
external fcnox external fcnox
common/psival/psinv common/psinv/psinv
xvec(1)=rz xvec(1)=rz
xvec(2)=zz xvec(2)=zz
tol = sqrt(comp_eps) tol = sqrt(comp_eps)
@ -1479,7 +1472,7 @@ c
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
dimension x(n),fvec(n),fjac(ldfjac,n) dimension x(n),fvec(n),fjac(ldfjac,n)
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/psival/psinv common/psinv/psinv
common/cnpsi/h common/cnpsi/h
common/pareq1/psia common/pareq1/psia
call equinum(x(1),x(2)) call equinum(x(1),x(2))
@ -1505,7 +1498,7 @@ c
common/psinr/psinr common/psinr/psinr
common/qpsi/qpsi common/qpsi/qpsi
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/dens/dens,ddens common/dddens/dens,ddens
c c
write(645,*) ' #psi rhot ne Te q Jphi' write(645,*) ' #psi rhot ne Te q Jphi'
psin=0.0d0 psin=0.0d0
@ -1739,14 +1732,14 @@ c
rhotnr(k+1)=rhotnr(k)+drhot rhotnr(k+1)=rhotnr(k)+drhot
end do end do
rhotsx=rhotnr(nr) rhotsx=rhotnr(nr)
do k=1,nr do k=2,nr
rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr)) rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr))
end do end do
c c
c spline interpolation of rhotor c spline interpolation of rhotor
c c
iopt=0 iopt=0
call difcsg(psinr,rhotnr,nr,iopt,crhot,ier) call difcsn(psinr,rhotnr,nnw,nr,iopt,crhot,ier)
return return
end end
@ -1757,11 +1750,10 @@ c
common/psinr/psinr common/psinr/psinr
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/cq/cq common/cq/cq
irt=int((nrho-1)*psi+1) call vlocate(psinr,nrho,psi,irt)
if(irt.eq.0) irt=1 irt=max(1,min(irt,nrho-1))
if(irt.eq.nrho) irt=nrho-1
dps=psi-psinr(irt) dps=psi-psinr(irt)
fq_eq=fspli(cq,nrho,irt,dps) fq_eq=fspli(cq,nnw,irt,dps)
return return
end end
@ -1773,11 +1765,10 @@ c
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/crhot/crhot common/crhot/crhot
c c
irt=int((nrho-1)*psi+1) call vlocate(psinr,nrho,psi,irt)
if(irt.eq.0) irt=1 irt=max(1,min(irt,nrho-1))
if(irt.eq.nrho) irt=nrho-1
dps=psi-psinr(irt) dps=psi-psinr(irt)
frhotor_eq=fspli(crhot,nrho,irt,dps) frhotor_eq=fspli(crhot,nnw,irt,dps)
return return
end end
@ -1796,11 +1787,13 @@ c
common/crhotq/crhotq common/crhotq/crhotq
rpsi=sqrt(psi) rpsi=sqrt(psi)
ip=int((nintp-1)*rpsi+1) c ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1 c if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1 c if(ip.eq.nintp) ip=nintp-1
call vlocate(rpstab,nintp,rpsi,ip)
ip=max(1,min(ip,nintp-1))
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
frhotor_av=fspli(crhotq,nintp,ip,dps) frhotor_av=fspli(crhotq,nnintp,ip,dps)
return return
end end
@ -1816,13 +1809,17 @@ c
common/coffrp/crp common/coffrp/crp
dr=1.0d0/dble(nnr-1) dr=1.0d0/dble(nnr-1)
do i=1,nnr do i=2,nnr-1
rhop(i)=(i-1)*dr rhop(i)=(i-1)*dr
psin=rhop(i)*rhop(i) psin=rhop(i)*rhop(i)
rhot(i)=frhotor(psin) rhot(i)=frhotor_eq(psin)
wp(i)=1.0d0 wp(i)=1.0d0
end do end do
rhop(1)=0.0d0
rhot(1)=0.0d0
wp(1)=1.0d3 wp(1)=1.0d3
rhop(nnr)=1.0d0
rhot(nnr)=1.0d0
wp(nnr)=1.0d3 wp(nnr)=1.0d3
c spline interpolation of rhopol versus rhotor c spline interpolation of rhopol versus rhotor
@ -2156,7 +2153,7 @@ c
common/cratj/cratja,cratjb,cratjpl common/cratj/cratja,cratjb,cratjpl
common/crhotq/crhotq common/crhotq/crhotq
common/cnt/rup,zup,rlw,zlw common/cnt/rup,zup,rlw,zlw
common/bound/zbmin,zbmax common/zbound/zbmin,zbmax
common/rarea/rarea common/rarea/rarea
common/phitedge/phitedge common/phitedge/phitedge
common/cdadrhot/cdadrhot common/cdadrhot/cdadrhot
@ -2478,11 +2475,13 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
common/pstab/rpstab common/pstab/rpstab
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/cdadrhot/cdadrhot common/cdadrhot/cdadrhot
ip=int((nintp-1)*rpsi+1) c ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1 c if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1 c if(ip.eq.nintp) ip=nintp-1
call vlocate(rpstab,nintp,rpsi,ip)
ip=max(1,min(ip,nintp-1))
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
fdadrhot=fspli(cdadrhot,nintp,ip,dps) fdadrhot=fspli(cdadrhot,nnintp,ip,dps)
return return
end end
@ -2493,11 +2492,13 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
common/pstab/rpstab common/pstab/rpstab
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/cdvdrhot/cdvdrhot common/cdvdrhot/cdvdrhot
ip=int((nintp-1)*rpsi+1) c ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1 c if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1 c if(ip.eq.nintp) ip=nintp-1
call vlocate(rpstab,nintp,rpsi,ip)
ip=max(1,min(ip,nintp-1))
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
fdvdrhot=fspli(cdvdrhot,nintp,ip,dps) fdvdrhot=fspli(cdvdrhot,nnintp,ip,dps)
return return
end end
@ -3140,7 +3141,7 @@ c
common/dxgyg/derxg,deryg common/dxgyg/derxg,deryg
common/iieq/iequil common/iieq/iequil
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/psival/psinv common/psinv/psinv
common/sgnib/sgnbphi,sgniphi common/sgnib/sgnbphi,sgniphi
c c
xg=0.0d0 xg=0.0d0
@ -3292,14 +3293,14 @@ c
c c
common/parqq/q0,qa,alq common/parqq/q0,qa,alq
common/parban/b0,rr0m,zr0m,rpam common/parban/b0,rr0m,zr0m,rpam
common/psival/psinv common/psinv/psinv
common/pareq1/psia common/pareq1/psia
common/densbnd/psdbnd common/densbnd/psdbnd
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/xgxg/xg common/xgxg/xg
common/dxgdps/dxgdpsi common/dxgdps/dxgdpsi
common/xgcn/xgcn common/xgcn/xgcn
common/dens/dens,ddens common/dddens/dens,ddens
common/sgnib/sgnbphi,sgniphi common/sgnib/sgnbphi,sgniphi
common/bmxmn/bmxi,bmni common/bmxmn/bmxi,bmni
common/fc/fci common/fc/fci
@ -3381,7 +3382,7 @@ c
dimension tfp(nrest),cfp(nrest),wrkfd(nrest) dimension tfp(nrest),cfp(nrest),wrkfd(nrest)
c c
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
common/psival/psinv common/psinv/psinv
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
common/pareq1/psia common/pareq1/psia
common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
@ -3415,7 +3416,7 @@ c
zzs(1)=zpsim zzs(1)=zpsim
nsr=nsrt nsr=nsrt
nsz=nszt nsz=nszt
call fpbisp(tr,nsr,tz,nsz,ccspl,3,3, call dierckx_fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
. rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
psinv=(ffspl(1)-psinop)/psiant psinv=(ffspl(1)-psinop)/psiant
if(psinv.lt.0.0d0) if(psinv.lt.0.0d0)
@ -3427,7 +3428,8 @@ c
kkz=3-nuz kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4) iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz, call dierckx_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)) . rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2))
dpsidr= ffspl(1)*psia dpsidr= ffspl(1)*psia
c c
@ -3437,7 +3439,8 @@ c
kkz=3-nuz kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4) iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz, call dierckx_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)) . rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2))
dpsidz= ffspl(1)*psia dpsidz= ffspl(1)*psia
c c
@ -3447,7 +3450,8 @@ c
kkz=3-nuz kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4) iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz, call dierckx_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)) . rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2))
ddpsidrr= ffspl(1)*psia ddpsidrr= ffspl(1)*psia
c c
@ -3457,7 +3461,8 @@ c
kkz=3-nuz kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4) iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz, call dierckx_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)) . rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2))
ddpsidzz= ffspl(1)*psia ddpsidzz= ffspl(1)*psia
c c
@ -3467,7 +3472,8 @@ c
kkz=3-nuz kkz=3-nuz
iwr=1+(nr-nur-4)*(nz-nuz-4) iwr=1+(nr-nur-4)*(nz-nuz-4)
iwz=iwr+4-nur iwz=iwr+4-nur
call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz, call dierckx_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)) . rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2))
ddpsidrz= ffspl(1)*psia ddpsidrz= ffspl(1)*psia
c c
@ -3548,12 +3554,12 @@ c
c c
subroutine sub_xg_derxg subroutine sub_xg_derxg
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
common/psival/psinv common/psinv/psinv
common/pareq1/psia common/pareq1/psia
common/xgxg/xg common/xgxg/xg
common/dxgdps/dxgdpsi common/dxgdps/dxgdpsi
common/xgcn/xgcn common/xgcn/xgcn
common/dens/dens,ddenspsin common/dddens/dens,ddenspsin
xg=0.0d0 xg=0.0d0
dxgdpsi=0.0d0 dxgdpsi=0.0d0
c if(psinv.le.psdbnd.and.psinv.ge.0) then c if(psinv.le.psdbnd.and.psinv.ge.0) then
@ -3576,7 +3582,7 @@ c
common/denspp/psnpp,aad,bbd,ccd common/denspp/psnpp,aad,bbd,ccd
common/iipr/iprof common/iipr/iprof
common/pardens/dens0,aln1,aln2 common/pardens/dens0,aln1,aln2
common/dens/dens,ddens common/dddens/dens,ddens
common/coffdt/tfn common/coffdt/tfn
common/coffdnst/nsfd common/coffdnst/nsfd
common/cofffn/cfn common/cofffn/cfn
@ -3644,8 +3650,7 @@ c
temperature=(te0-dte0)*proft+dte0 temperature=(te0-dte0)*proft+dte0
else else
call vlocate(psrad,npp,arg,k) call vlocate(psrad,npp,arg,k)
if(k.eq.0) k=1 k=max(1,min(k,npp-1))
if(k.eq.npp) k=npp-1
dps=arg-psrad(k) dps=arg-psrad(k)
temperature=fspli(ct,npmx,k,dps) temperature=fspli(ct,npmx,k,dps)
endif endif
@ -3672,8 +3677,7 @@ c
fzeff=zeff fzeff=zeff
else else
call vlocate(psrad,npp,ps,k) call vlocate(psrad,npp,ps,k)
if(k.eq.0) k=1 k=max(1,min(k,npp-1))
if(k.eq.npp) k=npp-1
dps=ps-psrad(k) dps=ps-psrad(k)
fzeff=fspli(cz,npmx,k,dps) fzeff=fspli(cz,npmx,k,dps)
endif endif
@ -4238,17 +4242,17 @@ c
c c
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
c c
areai=fspli(carea,nintp,ip,dps) areai=fspli(carea,nnintp,ip,dps)
voli=fspli(cvol,nintp,ip,dps) voli=fspli(cvol,nnintp,ip,dps)
dervoli=fsplid(cvol,nintp,ip,dps) dervoli=fsplid(cvol,nnintp,ip,dps)
rrii=fspli(crri,nintp,ip,dps) rrii=fspli(crri,nnintp,ip,dps)
c c
if(intp.eq.0) return if(intp.eq.0) return
c c
rbavi=fspli(crbav,nintp,ip,dps) rbavi=fspli(crbav,nnintp,ip,dps)
bmxi=fspli(cbmx,nintp,ip,dps) bmxi=fspli(cbmx,nnintp,ip,dps)
bmni=fspli(cbmn,nintp,ip,dps) bmni=fspli(cbmn,nnintp,ip,dps)
fci=fspli(cfc,nintp,ip,dps) fci=fspli(cfc,nnintp,ip,dps)
c c
return return
end end
@ -4267,9 +4271,9 @@ c
if(ip.eq.0) ip=1 if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1 if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip) dps=rpsi-rpstab(ip)
ratjai=fspli(cratja,nintp,ip,dps) ratjai=fspli(cratja,nnintp,ip,dps)
ratjbi=fspli(cratjb,nintp,ip,dps) ratjbi=fspli(cratjb,nnintp,ip,dps)
ratjpli=fspli(cratjpl,nintp,ip,dps) ratjpli=fspli(cratjpl,nnintp,ip,dps)
return return
end end
c c
@ -4290,7 +4294,7 @@ c
common/parban/b0,rr0m,zr0m,rpam common/parban/b0,rr0m,zr0m,rpam
common/absor/alpha,effjcd,akim,tau0 common/absor/alpha,effjcd,akim,tau0
c c
common/psival/psinv common/psinv/psinv
common/sgnib/sgnbphi,sgniphi common/sgnib/sgnbphi,sgniphi
common/bmxmn/bmxi,bmni common/bmxmn/bmxi,bmni
common/fc/fci common/fc/fci
@ -4364,7 +4368,7 @@ c
c c
common/absor/alpha,effjcd,akim,tau common/absor/alpha,effjcd,akim,tau
c c
common/psival/psinv common/psinv/psinv
common/tete/tekev common/tete/tekev
common/amut/amu common/amut/amu
common/zz/Zeff common/zz/Zeff
@ -5381,7 +5385,7 @@ c
common/ieccd/ieccd common/ieccd/ieccd
common/tete/tekev common/tete/tekev
common/dens/dens,ddens common/dddens/dens,ddens
common/zz/Zeff common/zz/Zeff
common/btot/btot common/btot/btot
common/bmxmn/bmax,bmin common/bmxmn/bmax,bmin
@ -5713,7 +5717,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/fc/fc common/fc/fc
common/ncl/hb common/ncl/hb
common/psival/psinv common/psinv/psinv
common/amut/amu common/amut/amu
common/tete/tekev common/tete/tekev
common/zz/Zeff common/zz/Zeff
@ -5747,7 +5751,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(kwrk=nnintp+nlam+njest+nlest+3)
parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)
external fpbisp external dierckx_fpbisp
dimension xxs(1),yys(1),ffs(1) dimension xxs(1),yys(1),ffs(1)
dimension ch01(lw01),ch((njest-4)*(nlest-4)) dimension ch01(lw01),ch((njest-4)*(nlest-4))
@ -5764,13 +5768,13 @@ c gg=F(u)/u with F(u) as in Cohen paper
xxs(1)=sqrt(psi) xxs(1)=sqrt(psi)
yys(1)=alam yys(1)=alam
call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs, call dierckx_fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs,
. wrk(1),wrk(5),iwrk(1),iwrk(2)) . wrk(1),wrk(5),iwrk(1),iwrk(2))
fv=ffs(1) fv=ffs(1)
iwp=1+(njp-4)*(nlm-5) iwp=1+(njp-4)*(nlm-5)
iwl=iwp+4 iwl=iwp+4
call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2, call dierckx_fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2,
. xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2)) . xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2))
dfv=ffs(1) dfv=ffs(1)
@ -6580,7 +6584,7 @@ c wave vector and electric field after reflection in lab frame
common/limiter/rlim,zlim,nlim common/limiter/rlim,zlim,nlim
common/anv/anv common/anv/anv
common/dsds/dst common/dsds/dst
common/psival/psinv common/psinv/psinv
common/densbnd/psdbnd common/densbnd/psdbnd
common/dstvac/dstvac common/dstvac/dstvac
c ivac=1 first interface plasma-vacuum c ivac=1 first interface plasma-vacuum

12
src/gray_main.f90
View File

@ -33,21 +33,11 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
! read data plus initialization ! read data plus initialization
index_rt=1 index_rt=1
print*,'GRAY started'
! call myflush
call prfile call prfile
print*,' file headers written'
! call myflush
call paraminit call paraminit
print*,' variables initialized'
! call myflush
call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, & call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, &
nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin) nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin)
print*,' spline computed'
! call myflush
call vectinit call vectinit
print*,' beam arrays allocated'
! call myflush
if(iercom.eq.0) then if(iercom.eq.0) then
if(igrad.eq.0) call ic_rt if(igrad.eq.0) call ic_rt
if(igrad.gt.0) call ic_gb if(igrad.gt.0) call ic_gb
@ -57,8 +47,6 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
write(*,*) ' IERR = ', ierr write(*,*) ' IERR = ', ierr
return return
end if end if
print*,' initial conditions set'
! call myflush
! beam/ray propagation ! beam/ray propagation
call gray_integration call gray_integration

2
src/grayl.f
View File

@ -7733,7 +7733,7 @@ c we partition the working space and evaluate the partial derivative
end end
subroutine dierckx_coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy, subroutine coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,
* wrk,lwrk,ier) * wrk,lwrk,ier)
c subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,... c subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
c ,my the partial derivative ( order nux,nuy) of a bivariate spline c ,my the partial derivative ( order nux,nuy) of a bivariate spline