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

changes in rhotor interpolation merged from trunk into nocommon

Browse Source
This commit is contained in:
Lorenzo Figini 2015-06-19 10:22:49 +00:00
commit 771cdb3822
2 changed files with 30 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
src/gray.f
View File

@ -1204,15 +1204,12 @@ c local variables
real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,current, real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,current,
. dpsinr,psia0,fp,xb,xe,ssfp,psinxp,rmop,zmop,psinoptmp,r10,z10, . dpsinr,psia0,fp,xb,xe,ssfp,psinxp,rmop,zmop,psinoptmp,r10,z10,
. rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax, . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax,
. rhot15,psi15,rhot2,psi2 . rhot15,psi15,rhot2,psi2,rhotsx
real(wp_), dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol, real(wp_), dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol,
. wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim
character(len=48) :: stringa character(len=48) :: stringa
c common/external functions/variables c common/external functions/variables
real(wp_) :: rhotsx
real(wp_) :: frhotor real(wp_) :: frhotor
c
common/rhotsx/rhotsx
c c
c read from file eqdsk c read from file eqdsk
c (see http://fusion.gat.com/efit/g_eqdsk.html) c (see http://fusion.gat.com/efit/g_eqdsk.html)
@ -1585,7 +1582,7 @@ c compute B_toroidal on axis
print'(a,f8.4)','BT_axis = ',btaxis print'(a,f8.4)','BT_axis = ',btaxis
c compute normalized rho_tor from eqdsk q profile c compute normalized rho_tor from eqdsk q profile
call rhotor(nr) call rhotor(rhotsx)
phitedge=abs(psia)*rhotsx*2*pi phitedge=abs(psia)*rhotsx*2*pi
c rrtor=sqrt(phitedge/abs(btrcen)/pi) c rrtor=sqrt(phitedge/abs(btrcen)/pi)
call rhopol call rhopol
@ -2065,21 +2062,17 @@ c
c c
c c
c c
subroutine rhotor(nnr) subroutine rhotor(rhotsx)
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use interp_eqprof, only : nr,psinr,qpsi,crhot,cq use interp_eqprof, only : nr,psinr,rhopnr,qpsi,crhot,cq
use simplespline, only : difcsn use simplespline, only : difcs
implicit none implicit none
c arguments c arguments
integer :: nnr real(wp_), intent(out) :: rhotsx
c local variables c local variables
integer :: iopt,ier,k integer :: iopt,ier,k
real(wp_) :: dx,drhot real(wp_) :: dx,drhot
real(wp_), dimension(nr) :: rhotnr real(wp_), dimension(nr) :: rhotnr
c common/external functions/variables
real(wp_) :: rhotsx
c
common/rhotsx/rhotsx
c c
if(allocated(cq)) deallocate(cq) if(allocated(cq)) deallocate(cq)
if(allocated(crhot)) deallocate(crhot) if(allocated(crhot)) deallocate(crhot)
@ -2088,25 +2081,25 @@ c
c normalized toroidal rho : ~ Integral q(psi) dpsi c normalized toroidal rho : ~ Integral q(psi) dpsi
c c
iopt=0 iopt=0
call difcsn(psinr,qpsi,nr,nnr,iopt,cq,ier) call difcs(psinr,qpsi,nr,iopt,cq,ier)
c c
rhotnr(1)=0.0_wp_ rhotnr(1)=0.0_wp_
do k=1,nnr-1 do k=1,nr-1
dx=psinr(k+1)-psinr(k) dx=psinr(k+1)-psinr(k)
drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0_wp_+dx*(cq(k,3)/3.0_wp_ drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0_wp_+dx*(cq(k,3)/3.0_wp_
. +dx*cq(k,4)/4.0_wp_))) . +dx*cq(k,4)/4.0_wp_)))
rhotnr(k+1)=rhotnr(k)+drhot rhotnr(k+1)=rhotnr(k)+drhot
end do end do
rhotsx=rhotnr(nnr) rhotsx=rhotnr(nr)
do k=1,nr do k=1,nr
rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nnr)) 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 difcsn(psinr,rhotnr,nr,nnr,iopt,crhot,ier) rhopnr=sqrt(psinr)
c call difcs(rhopnr,rhotnr,nr,iopt,crhot,ier)
return return
end end
c c
@ -2134,23 +2127,22 @@ c
c c
c c
c c
function frhotor_eq(psi) function frhotor_eq(rhop)
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use interp_eqprof, only : psinr,nr,crhot use interp_eqprof, only : rhopnr,nr,crhot
use simplespline, only :spli use simplespline, only : spli
use utils, only : locate
implicit none implicit none
c arguments c arguments
real(wp_) :: psi,frhotor_eq real(wp_) :: rhop,frhotor_eq
c local variables c local variables
integer :: irt integer :: irt
real(wp_) :: dps real(wp_) :: drh
c c
irt=int((nr-1)*psi+1) call locate(rhopnr,nr,rhop,irt)
c if(irt.eq.0) irt=1
c if(irt.eq.nr) irt=nr-1
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)
c c
return return
end end
@ -2166,7 +2158,7 @@ c arguments
c common/external functions/variables c common/external functions/variables
real(wp_) :: frhotor_eq,frhotor_an real(wp_) :: frhotor_eq,frhotor_an
c c
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))
c c
return return
@ -2779,9 +2771,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 computation of rhot from calculated q profile c computation of rhot from calculated q profile
@ -2846,7 +2838,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
rpstab(jp)=1.0_wp_ rpstab(jp)=1.0_wp_
pstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_
end if end if
rhot_eq(jp)=frhotor_eq(pstab(jp)) rhot_eq(jp)=frhotor_eq(sqrt(pstab(jp)))
write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), write(56,99) pstab(jp),rhot_eq(jp),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),
@ -3410,6 +3402,7 @@ c
function frhotor_an(rhop) function frhotor_an(rhop)
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use dierckx, only : splev use dierckx, only : splev
use utils, only : locate
implicit none implicit none
c local contants c local contants
integer, parameter :: nnr=101,nrest=nnr+4 integer, parameter :: nnr=101,nrest=nnr+4

7
src/interp_eqprof.f90
View File

@ -10,7 +10,7 @@ module interp_eqprof
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz,tfp REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz,tfp
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,cfp REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,cfp
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq01,cceq10,cceq20,cceq02,cceq11 REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq01,cceq10,cceq20,cceq02,cceq11
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,qpsi,rv,zv REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopnr,qpsi,rv,zv
REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin,psi,btotal,crhot,cq REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin,psi,btotal,crhot,cq
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim,rbbbs,zbbbs REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim,rbbbs,zbbbs
@ -43,8 +43,8 @@ contains
call dealloc_equilvec call dealloc_equilvec
allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), & allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), &
btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), & btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), &
cceq20(lw20),cceq11(lw11),psi(nr,nz),psin(nr,nz),psinr(nr),qpsi(nr), & cceq20(lw20),cceq11(lw11),psi(nr,nz),psin(nr,nz),psinr(nr),rhopnr(nr), &
stat=ier) qpsi(nr), stat=ier)
if (ier/=0) call dealloc_equilvec if (ier/=0) call dealloc_equilvec
end subroutine alloc_equilvec end subroutine alloc_equilvec
@ -66,6 +66,7 @@ contains
if(allocated(psi)) deallocate(psi) if(allocated(psi)) deallocate(psi)
if(allocated(psin)) deallocate(psin) if(allocated(psin)) deallocate(psin)
if(allocated(psinr)) deallocate(psinr) if(allocated(psinr)) deallocate(psinr)
if(allocated(rhopnr)) deallocate(rhopnr)
if(allocated(qpsi)) deallocate(qpsi) if(allocated(qpsi)) deallocate(qpsi)
end subroutine dealloc_equilvec end subroutine dealloc_equilvec