split profdata routine and moved in coreprof module
This commit is contained in:
parent
667f6fd111
commit
b9b6d3e8f4
@ -1,50 +1,15 @@
|
||||
module coreprofiles
|
||||
use const_and_precisions, only : wp_
|
||||
use const_and_precisions, only : wp_,zero,one
|
||||
implicit none
|
||||
|
||||
INTEGER, SAVE :: npp,nsfd
|
||||
REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
|
||||
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psrad,derad,terad,zfc
|
||||
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn
|
||||
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad
|
||||
REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
|
||||
|
||||
contains
|
||||
|
||||
subroutine alloc_profvec(ier)
|
||||
implicit none
|
||||
integer, intent(out) :: ier
|
||||
integer :: npest
|
||||
|
||||
ier=0
|
||||
if(npp.le.0) then
|
||||
ier = -1
|
||||
return
|
||||
end if
|
||||
|
||||
npest=npp+4
|
||||
|
||||
call dealloc_profvec
|
||||
allocate(psrad(npp),terad(npp),derad(npp),zfc(npp),ct(npp,4), &
|
||||
cz(npp,4),tfn(npest),cfn(npest), &
|
||||
stat=ier)
|
||||
if (ier/=0) call dealloc_profvec
|
||||
end subroutine alloc_profvec
|
||||
|
||||
subroutine dealloc_profvec
|
||||
implicit none
|
||||
if(allocated(psrad)) deallocate(psrad)
|
||||
if(allocated(terad)) deallocate(terad)
|
||||
if(allocated(derad)) deallocate(derad)
|
||||
if(allocated(zfc)) deallocate(zfc)
|
||||
if(allocated(ct)) deallocate(ct)
|
||||
if(allocated(cz)) deallocate(cz)
|
||||
if(allocated(tfn)) deallocate(tfn)
|
||||
if(allocated(cfn)) deallocate(cfn)
|
||||
|
||||
end subroutine dealloc_profvec
|
||||
|
||||
subroutine density(psin,dens,ddens)
|
||||
use const_and_precisions, only : wp_
|
||||
use graydata_flags, only : iprof
|
||||
use graydata_anequil, only : dens0,aln1,aln2
|
||||
use dierckx, only : splev,splder
|
||||
@ -62,16 +27,16 @@ contains
|
||||
!
|
||||
! computation of density [10^19 m^-3] and derivative wrt psi
|
||||
!
|
||||
dens=0.0_wp_
|
||||
ddens=0.0_wp_
|
||||
if(psin.ge.psdbnd.or.psin.lt.0.0_wp_) return
|
||||
dens=zero
|
||||
ddens=zero
|
||||
if(psin.ge.psdbnd.or.psin.lt.zero) return
|
||||
!
|
||||
if(iprof.eq.0) then
|
||||
if(psin.gt.1.0_wp_) return
|
||||
profd=(1.0_wp_-psin**aln1)**aln2
|
||||
if(psin.gt.one) return
|
||||
profd=(one-psin**aln1)**aln2
|
||||
dens=dens0*profd
|
||||
dprofd=-aln1*aln2*psin**(aln1-1.0_wp_) &
|
||||
*(1.0_wp_-psin**aln1)**(aln2-1.0_wp_)
|
||||
dprofd=-aln1*aln2*psin**(aln1-one) &
|
||||
*(one-psin**aln1)**(aln2-one)
|
||||
ddens=dens0*dprofd
|
||||
else
|
||||
if(psin.gt.psnpp) then
|
||||
@ -86,8 +51,8 @@ contains
|
||||
fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
|
||||
dfp=ddenpp+dpsib*d2denpp
|
||||
tt=dpsib/(psdbnd-psnpp)
|
||||
fh=(1.0_wp_-tt)**3*(1.0_wp_+3.0_wp_*tt+6.0_wp_*tt**2)
|
||||
dfh=-30.0_wp_*(1.0_wp_-tt)**2*tt**2/(psdbnd-psnpp)
|
||||
fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2)
|
||||
dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp)
|
||||
dens=fp*fh
|
||||
ddens=dfp*fh+fp*dfh
|
||||
else
|
||||
@ -100,9 +65,9 @@ contains
|
||||
call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
|
||||
ddens=ffs(1)
|
||||
if(ier.gt.0) print*,ier
|
||||
if(abs(dens).lt.1.0e-10_wp_) dens=0.0_wp_
|
||||
if(abs(dens).lt.1.0e-10_wp_) dens=zero
|
||||
end if
|
||||
if(dens.lt.0.0_wp_) print*,' DENSITY NEGATIVE',dens
|
||||
if(dens.lt.zero) print*,' DENSITY NEGATIVE',dens
|
||||
end if
|
||||
end subroutine density
|
||||
|
||||
@ -159,4 +124,147 @@ contains
|
||||
endif
|
||||
end function fzeff
|
||||
|
||||
subroutine read_profiles(filenm,psin,te,ne,zeff,unit)
|
||||
use utils, only : get_free_unit
|
||||
implicit none
|
||||
! arguments
|
||||
character(len=*), intent(in) :: filenm
|
||||
real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff
|
||||
integer, optional, intent(in) :: unit
|
||||
! local variables
|
||||
integer :: u, i, n
|
||||
|
||||
if (present(unit)) then
|
||||
u=unit
|
||||
else
|
||||
u=get_free_unit()
|
||||
end if
|
||||
open(file=trim(filenm),status='old',unit=u)
|
||||
read(u,*) n
|
||||
if(allocated(psin)) deallocate(psin)
|
||||
if(allocated(te)) deallocate(te)
|
||||
if(allocated(ne)) deallocate(ne)
|
||||
if(allocated(zeff)) deallocate(zeff)
|
||||
allocate(psin(n),te(n),ne(n),zeff(n))
|
||||
do i=1,n
|
||||
read(u,*) psin(i),te(i),ne(i),zeff(i)
|
||||
end do
|
||||
psin(1)=max(psin(1),zero)
|
||||
close(u)
|
||||
end subroutine read_profiles
|
||||
|
||||
subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal)
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), dimension(:), intent(inout) :: te,ne
|
||||
real(wp_), intent(in) :: tfact,nfact,bfact
|
||||
integer, intent(in) :: iscal
|
||||
! local variables
|
||||
real(wp_) :: aat,aan,ffact
|
||||
|
||||
if (iscal==0) then
|
||||
aat=2.0_wp_/3.0_wp_
|
||||
aan=4.0_wp_/3.0_wp_
|
||||
else
|
||||
aat=1.0_wp_
|
||||
aan=1.0_wp_
|
||||
end if
|
||||
if(iscal==2) then
|
||||
ffact=1.0_wp_
|
||||
else
|
||||
ffact=bfact
|
||||
end if
|
||||
te(:)=te(:)*ffact**aat*tfact
|
||||
ne(:)=ne(:)*ffact**aan*nfact
|
||||
end subroutine tene_scal
|
||||
|
||||
subroutine set_prfspl(psin,te,ne,zeff,ssplne)
|
||||
use simplespline, only : difcs
|
||||
use dierckx, only : curfit, splev, splder
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff
|
||||
real(wp_), intent(in) :: ssplne
|
||||
! local variables
|
||||
integer, parameter :: iopt=0, kspl=3
|
||||
integer :: n, npest, lwrkf, ier
|
||||
real(wp_) :: xb, xe, fp, xnv, ynv
|
||||
real(wp_), dimension(:), allocatable :: wf, wrkf
|
||||
integer, dimension(:), allocatable :: iwrkf
|
||||
real(wp_), dimension(1) :: dedge,ddedge,d2dedge
|
||||
|
||||
n=size(psin)
|
||||
npest=n+4
|
||||
lwrkf=n*4+npest*16
|
||||
allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
|
||||
|
||||
! if necessary, reallocate spline arrays
|
||||
if(.not.allocated(psrad)) then
|
||||
allocate(psrad(n),ct(n,4),cz(n,4))
|
||||
else
|
||||
if(size(psrad)<n) then
|
||||
deallocate(psrad,ct,cz)
|
||||
allocate(psrad(n),ct(n,4),cz(n,4))
|
||||
end if
|
||||
end if
|
||||
if(.not.allocated(cfn)) then
|
||||
allocate(tfn(npest),cfn(npest))
|
||||
else
|
||||
if(size(cfn)<npest) then
|
||||
deallocate(tfn,cfn)
|
||||
allocate(tfn(npest),cfn(npest))
|
||||
end if
|
||||
end if
|
||||
|
||||
! spline approximation of temperature and Zeff
|
||||
call difcs(psin,te, n,iopt,ct,ier)
|
||||
call difcs(psin,zeff,n,iopt,cz,ier)
|
||||
psrad=psin
|
||||
npp=n
|
||||
|
||||
! spline approximation of density
|
||||
xb=zero
|
||||
xe=psin(n)
|
||||
wf(:)=one
|
||||
call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne,npest, &
|
||||
nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
|
||||
|
||||
! compute polinomial extrapolation matching the spline boundary up to the
|
||||
! 2nd order derivative, extending the profile up to psi=psdbnd where
|
||||
! ne=ne'=ne''=0
|
||||
! spline value and derivatives at the edge
|
||||
call splev(tfn,nsfd,cfn,3,psin(n:n),dedge(1:1),1,ier)
|
||||
call splder(tfn,nsfd,cfn,3,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
|
||||
call splder(tfn,nsfd,cfn,3,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
|
||||
! determination of the boundary
|
||||
psnpp=psin(n)
|
||||
denpp=dedge(1)
|
||||
ddenpp=ddedge(1)
|
||||
d2denpp=d2dedge(1)
|
||||
psdbnd=psnpp
|
||||
if(ddenpp.lt.zero) then
|
||||
xnv=psnpp-ddenpp/d2denpp
|
||||
ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
|
||||
if(d2denpp.gt.zero.and.ynv.ge.zero) then
|
||||
psdbnd=xnv
|
||||
else
|
||||
psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
|
||||
end if
|
||||
print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
|
||||
print*,psdbnd
|
||||
end if
|
||||
|
||||
deallocate(iwrkf,wrkf,wf)
|
||||
end subroutine set_prfspl
|
||||
|
||||
subroutine unset_prfspl
|
||||
implicit none
|
||||
|
||||
if(allocated(psrad)) deallocate(psrad)
|
||||
if(allocated(ct)) deallocate(ct)
|
||||
if(allocated(cz)) deallocate(cz)
|
||||
if(allocated(tfn)) deallocate(tfn)
|
||||
if(allocated(cfn)) deallocate(cfn)
|
||||
end subroutine unset_prfspl
|
||||
|
||||
end module coreprofiles
|
||||
126
src/gray.f
126
src/gray.f
@ -657,15 +657,16 @@ c
|
||||
use graydata_flags
|
||||
use graydata_par
|
||||
use graydata_anequil
|
||||
use coreprofiles, only : psdbnd
|
||||
use coreprofiles, only : psdbnd,read_profiles,tene_scal,set_prfspl
|
||||
use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen
|
||||
use reflections, only : rlim,zlim,nlim,alloc_lim
|
||||
use beamdata, only : nrayr,nrayth,nstep
|
||||
implicit none
|
||||
c local variables
|
||||
integer :: nfil,iox,ierr,leqq,lprf
|
||||
integer :: nfil,iox,ierr,leqq
|
||||
real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
|
||||
. w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta
|
||||
real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc
|
||||
character(len=8) :: wdat
|
||||
character(len=10) :: wtim
|
||||
c common/external functions/variables
|
||||
@ -908,11 +909,11 @@ c
|
||||
|
||||
if (iprof.eq.1) then
|
||||
nprof=98
|
||||
lprf=len_trim(filenmprf)
|
||||
open(file=filenmprf(1:lprf)//'.prf',
|
||||
. status= 'unknown',unit=nprof)
|
||||
call profdata
|
||||
close(nprof)
|
||||
call read_profiles(trim(filenmprf)//'.prf',psrad,terad,derad,
|
||||
. zfc,nprof)
|
||||
call tene_scal(terad,derad,factt,factn,factb,iscal)
|
||||
call set_prfspl(psrad,terad,derad,zfc,0.001_wp_)
|
||||
deallocate(psrad,terad,derad,zfc)
|
||||
end if
|
||||
c
|
||||
c read equilibrium data from file if iequil=2
|
||||
@ -1905,117 +1906,6 @@ c
|
||||
113 format(i6,12(1x,e12.5))
|
||||
end
|
||||
c
|
||||
c
|
||||
subroutine profdata
|
||||
use const_and_precisions, only : wp_,zero
|
||||
use graydata_flags, only : iprof,iscal,nprof
|
||||
use graydata_par, only : factb,factt,factn
|
||||
use coreprofiles, only : nsfd,npp,psnpp,denpp,ddenpp,d2denpp,
|
||||
. psdbnd,psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec
|
||||
use simplespline, only : difcs
|
||||
use dierckx, only : curfit,splev,splder
|
||||
implicit none
|
||||
c local variables
|
||||
integer :: ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2
|
||||
integer, dimension(:), allocatable :: iwrkf
|
||||
real(wp_) :: aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi,
|
||||
. teradi,deradi,zfci,xb,xe,sspl,dpsb,fp
|
||||
real(wp_) :: xnv,ynv
|
||||
real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd
|
||||
real(wp_), dimension(1) :: densi,ddensi,d2densi
|
||||
c
|
||||
c read plasma profiles from file if iprof>0
|
||||
c
|
||||
if (iscal.eq.0) then
|
||||
aat=2.0_wp_/3.0_wp_
|
||||
aan=4.0_wp_/3.0_wp_
|
||||
else
|
||||
aat=1.0_wp_
|
||||
aan=1.0_wp_
|
||||
end if
|
||||
ffact=factb
|
||||
if(iscal.eq.2) ffact=1.0_wp_
|
||||
c
|
||||
if (iprof.gt.0) then
|
||||
read(nprof,*) npp
|
||||
c
|
||||
call alloc_profvec(ierr)
|
||||
if (ierr.ne.0) stop
|
||||
c
|
||||
npest=npp+4
|
||||
lwrkf=npp*4+npest*16
|
||||
if(allocated(wrkf)) deallocate(wrkf)
|
||||
if(allocated(iwrkf)) deallocate(iwrkf)
|
||||
if(allocated(wf)) deallocate(wf)
|
||||
if(allocated(wrkfd)) deallocate(wrkfd)
|
||||
allocate(wrkf(lwrkf),iwrkf(npest),wf(npp),wrkfd(npest))
|
||||
c
|
||||
read(nprof,*) psrad0,terad0,derad0,zfc0
|
||||
if(psrad0.ne.0.0_wp_) psrad0=0.0_wp_
|
||||
psrad(1)=psrad0
|
||||
terad(1)=terad0*ffact**aat*factt
|
||||
derad(1)=derad0*ffact**aan*factn
|
||||
zfc(1)=zfc0
|
||||
wf(1)=1.0_wp_
|
||||
do i=2,npp
|
||||
read(nprof,*) psradi,teradi,deradi,zfci
|
||||
psrad(i)=psradi
|
||||
terad(i)=teradi*ffact**aat*factt
|
||||
derad(i)=deradi*ffact**aan*factn
|
||||
zfc(i)=zfci
|
||||
wf(i)=1.0_wp_
|
||||
end do
|
||||
|
||||
c spline approximation of temperature and Zeff
|
||||
|
||||
iopt=0
|
||||
call difcs(psrad,terad,npp,iopt,ct,ier)
|
||||
|
||||
iopt=0
|
||||
call difcs(psrad,zfc,npp,iopt,cz,ier)
|
||||
|
||||
c spline approximation of density
|
||||
|
||||
iopt=0
|
||||
xb=0.0_wp_
|
||||
xe=psrad(npp)
|
||||
kspl=3
|
||||
sspl=.001_wp_
|
||||
|
||||
call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
|
||||
. tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
|
||||
|
||||
call splev(tfn,nsfd,cfn,3,psrad(npp),densi(1),1,ier)
|
||||
nu=1
|
||||
call splder(tfn,nsfd,cfn,3,nu,psrad(npp),ddensi(1),1,wrkfd,ier)
|
||||
nu=2
|
||||
call splder(tfn,nsfd,cfn,3,nu,psrad(npp),d2densi(1),1,wrkfd,ier)
|
||||
|
||||
psnpp=psrad(npp)
|
||||
denpp=densi(1)
|
||||
ddenpp=ddensi(1)
|
||||
d2denpp=d2densi(1)
|
||||
psdbnd=psnpp
|
||||
if(ddenpp.lt.0.0_wp_) then
|
||||
xnv=psnpp-ddenpp/d2denpp
|
||||
ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
|
||||
if(d2denpp.gt.zero.and.ynv.ge.zero) then
|
||||
psdbnd=xnv
|
||||
else
|
||||
psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
|
||||
end if
|
||||
print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
|
||||
print*,psdbnd
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
deallocate(iwrkf,wrkf,wf,wrkfd)
|
||||
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c
|
||||
subroutine rhotor(rhotsx)
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
@ -246,4 +246,33 @@ contains
|
||||
end do
|
||||
end subroutine bubble
|
||||
|
||||
function get_free_unit(umin,umax) result(i)
|
||||
implicit none
|
||||
integer :: i
|
||||
integer, intent(in), optional :: umin, umax
|
||||
integer, parameter :: max_allowed = 999
|
||||
integer :: ierr, iend
|
||||
logical :: ex, op
|
||||
|
||||
if (present(umin)) then
|
||||
i = max(0,umin) ! start searching from unit min
|
||||
else
|
||||
i = 0
|
||||
end if
|
||||
if (present(umax)) then
|
||||
iend = min(max(0,umax),max_allowed)
|
||||
else
|
||||
iend = max_allowed
|
||||
end if
|
||||
do
|
||||
if (i>iend) then
|
||||
i=-1 ! no free units found
|
||||
exit
|
||||
end if
|
||||
inquire(unit=i,exist=ex,opened=op,iostat=ierr)
|
||||
if (ierr==0.and.ex.and..not.op) exit ! unit i exists and is not open
|
||||
i = i + 1
|
||||
end do
|
||||
end function get_free_unit
|
||||
|
||||
end module utils
|
||||
Loading…
Reference in New Issue
Block a user