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

frhotor_eq now is a function of rhop

Browse Source
This commit is contained in:
Daniela Farina 2015-06-09 17:09:53 +00:00
parent a0275d8691
commit 75b77b01c1
1 changed files with 16 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/gray.f
View File

@ -1983,8 +1983,10 @@ c
subroutine rhotor(nr) subroutine rhotor(nr)
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nnw=501) parameter(nnw=501)
dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4) dimension psinr(nnw),rhopnr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),
* crhot(nnw,4)
common/psinr/psinr common/psinr/psinr
common/rhopnr/rhopnr
common/qpsi/qpsi common/qpsi/qpsi
common/rhotsx/rhotsx common/rhotsx/rhotsx
common/crhot/crhot common/crhot/crhot
@ -2010,7 +2012,8 @@ c
c spline interpolation of rhotor c spline interpolation of rhotor
c c
iopt=0 iopt=0
call difcs(psinr,rhotnr,nr,iopt,crhot,ier) rhopnr=sqrt(psinr)
call difcs(rhopnr,rhotnr,nr,iopt,crhot,ier)
return return
end end
@ -2029,27 +2032,28 @@ c
return return
end end
function frhotor_eq(psi) function frhotor_eq(rhop)
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nnw=501) parameter(nnw=501)
dimension psinr(nnw),crhot(nnw,4) dimension rhopnr(nnw),crhot(nnw,4)
common/psinr/psinr common/rhopnr/rhopnr
common/eqnn/nr,nz,npp,nintp common/eqnn/nr,nz,npp,nintp
common/crhot/crhot common/crhot/crhot
c c
irt=int((nr-1)*psi+1) c irt=int((nr-1)*psi+1)
c if(irt.eq.0) irt=1 c if(irt.eq.0) irt=1
c if(irt.eq.nr) irt=nr-1 c if(irt.eq.nr) irt=nr-1
call locate(rhopnr,nr,rhop,irt)
irt=min(max(1,irt),nr-1) irt=min(max(1,irt),nr-1)
dps=psi-psinr(irt) drh=rhop-rhopnr(irt)
frhotor_eq=spli(crhot,nr,irt,dps) frhotor_eq=spli(crhot,nr,irt,drh)
return return
end end
function frhotor(psi) function frhotor(psi)
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
common/iieq/iequil common/iieq/iequil
if(iequil.eq.2) frhotor=frhotor_eq(psi) if(iequil.eq.2) frhotor=frhotor_eq(sqrt(psi))
if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi))
return return
end end
@ -2609,9 +2613,9 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = <Jcd.b>
qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
qqv(jp)=qq qqv(jp)=qq
dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) dadrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2)
. *dadpsi/pi . *dadpsi/pi
dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) dvdrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2)
. *dvdpsi/pi . *dvdpsi/pi
c c
c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)
@ -2678,7 +2682,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
rpstab(jp)=1.0d0 rpstab(jp)=1.0d0
pstab(jp)=1.0d0 pstab(jp)=1.0d0
end if end if
rhot_eq=frhotor_eq(pstab(jp)) rhot_eq=frhotor_eq(sqrt(pstab(jp)))
write(56,99) pstab(jp),rhot_eq,rhotqv(jp), write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
. bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp),
. vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)