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

removed dependency on itm-modules. fixed uncorrect use of ipec in pec (ipec=-1). changed dubious assignments (isev, iind) in pec. cleaned dependencies in Makefile. added missing beamdata.f90 file to the repository. added main program to test with standard input files (eqdsk, prf). renamed ei function and dierckx subroutine to avoid conflicts with flush library.

Browse Source
This commit is contained in:
Lorenzo Figini 2014-06-16 16:41:43 +00:00
parent 9ca1ccd817
commit 074f331355
11 changed files with 643 additions and 330 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Makefile
View File

@ -35,7 +35,7 @@ FC=$(F90)
ifeq ("$(DBG)","")
FFLAGS= $(AUTO) -O3
else
FFLAGS= $(AUTO) -O0 -g
FFLAGS= $(AUTO) -O0 -g -Minform=inform -Mbounds -Mchkptr
endif
# Set include directories
@ -96,13 +96,11 @@ clean:
# Dependencies
# ------------
gray_main.o: gray-externals.o itm_types.o
gray-externals.o: green_func_p.o reflections.o scatterspl.o beamdata.o
gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o
green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o
itm_constants.o: itm_types.o
scatterspl.o: itm_types.o
beamdata.o: itm_types.o
scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o
# Special rule to preprocess source file and include svn revision
# ---------------------------------------------------------------

16
Makefile.standalone
View File

@ -1,10 +1,10 @@
# Executable name
EXE=gray
EXE=gray-jintrac
LIBFILE=lib$(EXE).a
# Objects list
OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \
beamdata.o green_func_p.o const_and_precisions.o itm_constants.o itm_types.o
beamdata.o green_func_p.o const_and_precisions.o
# Alternative search paths
vpath %.f90 src
@ -28,14 +28,12 @@ $(LIBFILE): $(OBJ)
ar -rv $@ $^
# Dependencies on modules
main.o: gray_main.o
gray_main.o: gray-externals.o itm_types.o
gray-externals.o: green_func_p.o reflections.o scatterspl.o beamdata.o
main.o: const_and_precisions.o
gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o
green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o
itm_constants.o: itm_types.o
scatterspl.o: itm_types.o
beamdata.o: itm_types.o
scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o
# General object compilation command
%.o: %.f90

72
src/beamdata.f90 Normal file
View File

@ -0,0 +1,72 @@
module beamdata
use const_and_precisions, only : r8
implicit none
integer, parameter :: jmx=31,kmx=36,nmx=8000
integer, save :: nrayr,nrayth,nstep
real(r8), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav
real(r8), dimension(:,:,:), allocatable, save :: pdjki,currj,didst
integer, dimension(:,:), allocatable, save :: iiv,iop,iow
real(r8), dimension(:,:), allocatable, save :: tau1v
real(r8), dimension(:), allocatable, save :: q
real(r8), dimension(:,:,:), allocatable, save :: yyrfl !(:,:,6)
real(r8), dimension(:,:,:), allocatable, save :: ywrk,ypwrk !(6,:,:)
real(r8), dimension(:,:,:), allocatable, save :: xc,xco,du1,du1o !(3,:,:)
real(r8), dimension(:,:,:), allocatable, save :: gri,dgrad2v !(3,:,:)
real(r8), dimension(:,:,:,:), allocatable, save :: ggri !(3,3,:,:)
real(r8), dimension(:,:), allocatable, save :: grad2
real(r8), dimension(:), allocatable, save :: dffiu,ddffiu
contains
subroutine alloc_beam(ierr)
implicit none
integer, intent(out) :: ierr
call dealloc_beam
allocate(psjki(nrayr,nrayth,nstep), ppabs(nrayr,nrayth,nstep), &
pdjki(nrayr,nrayth,nstep), ccci(nrayr,nrayth,nstep), &
currj(nrayr,nrayth,nstep), didst(nrayr,nrayth,nstep), &
tauv(nrayr,nrayth,nstep), alphav(nrayr,nrayth,nstep), &
iiv(nrayr,nrayth), iop(nrayr,nrayth), &
iow(nrayr,nrayth), tau1v(nrayr,nrayth), &
q(nrayr), yyrfl(nrayr,nrayth,6), &
ywrk(6,nrayr,nrayth), ypwrk(6,nrayr,nrayth), &
xc(3,nrayr,nrayth), xco(3,nrayr,nrayth), &
du1(3,nrayr,nrayth), du1o(3,nrayr,nrayth), &
gri(3,nrayr,nrayth), dgrad2v(3,nrayr,nrayth), &
ggri(3,3,nrayr,nrayth), grad2(nrayr,nrayth), &
dffiu(nrayr), ddffiu(nrayr), stat=ierr)
if (ierr/=0) call dealloc_beam
end subroutine alloc_beam
subroutine dealloc_beam
implicit none
if (allocated(psjki)) deallocate(psjki)
if (allocated(ppabs)) deallocate(ppabs)
if (allocated(pdjki)) deallocate(pdjki)
if (allocated(ccci)) deallocate(ccci)
if (allocated(currj)) deallocate(currj)
if (allocated(didst)) deallocate(didst)
if (allocated(tauv)) deallocate(tauv)
if (allocated(alphav)) deallocate(alphav)
if (allocated(iiv)) deallocate(iiv)
if (allocated(iop)) deallocate(iop)
if (allocated(iow)) deallocate(iow)
if (allocated(tau1v)) deallocate(tau1v)
if (allocated(q)) deallocate(q)
if (allocated(yyrfl)) deallocate(yyrfl)
if (allocated(ywrk)) deallocate(ywrk)
if (allocated(ypwrk)) deallocate(ypwrk)
if (allocated(xc)) deallocate(xc)
if (allocated(xco)) deallocate(xco)
if (allocated(du1)) deallocate(du1)
if (allocated(du1o)) deallocate(du1o)
if (allocated(gri)) deallocate(gri)
if (allocated(dgrad2v)) deallocate(dgrad2v)
if (allocated(ggri)) deallocate(ggri)
if (allocated(grad2)) deallocate(grad2)
if (allocated(dffiu)) deallocate(dffiu)
if (allocated(ddffiu)) deallocate(ddffiu)
end subroutine dealloc_beam
end module beamdata

40
src/const_and_precisions.f90
View File

@ -1,18 +1,22 @@
!########################################################################!
MODULE const_and_precisions
use itm_types, only : wp_ => r8
use itm_constants, only : pi => itm_pi, e_ => itm_qe, me_ => itm_me, c_ => itm_c
!########################################################################!
IMPLICIT NONE
PUBLIC
!------------------------------------------------------------------------
! common precisions
!------------------------------------------------------------------------
! INTEGER, PARAMETER :: sp_ = 4 ! single precision
! INTEGER, PARAMETER :: dp_ = 8 ! double precision
! INTEGER, PARAMETER :: wp_ = dp_ ! work-precision
! INTEGER, PARAMETER :: odep_ = dp_ ! ODE-solver precision
! INTEGER, PARAMETER :: I1 = SELECTED_INT_KIND (2) ! Integer*1
! INTEGER, PARAMETER :: I2 = SELECTED_INT_KIND (4) ! Integer*2
INTEGER, PARAMETER :: I4 = SELECTED_INT_KIND (9) ! Integer*4
! INTEGER, PARAMETER :: I8 = SELECTED_INT_KIND (18) ! Integer*8
! INTEGER, PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37) ! Real*4
INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8
INTEGER, PARAMETER :: wp_ = R8 ! work-precision
! INTEGER, PARAMETER :: odep_ = R8 ! ODE-solver precision
! INTEGER, PARAMETER :: xp_ = wp_ ! for ext. modules if necessary
!------------------------------------------------------------------------
! precisions which are in use in CONFIG_yat
@ -28,7 +32,7 @@
!========================================================================
REAL(wp_), PARAMETER :: zero = 0.0_wp_
REAL(wp_), PARAMETER :: unit = 1.0_wp_
! REAL(wp_), PARAMETER :: pi = 3.141592653589793_wp_
real (kind = R8), parameter :: pi = 3.141592653589793238462643383280_R8
! REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
! REAL(wp_), PARAMETER :: sqrt_2 = 1.414213562373095_wp_
! REAL(wp_), PARAMETER :: rad = pi/180.0_wp_
@ -47,7 +51,7 @@
!========================================================================
! Computer constants
!========================================================================
REAL(wp_), PARAMETER :: comp_eps = EPSILON(unit)
REAL(kind = R8), PARAMETER :: comp_eps = EPSILON(unit)
! REAL(wp_), PARAMETER :: comp_eps2 = comp_eps**2
! REAL(wp_), PARAMETER :: comp_tiny = TINY(unit)
! REAL(wp_), PARAMETER :: comp_huge = HUGE(unit)
@ -65,16 +69,24 @@
!========================================================================
! Physical constants (SI)
!========================================================================
! REAL(wp_), PARAMETER :: e_ = 1.602176487d-19 ! [C]
! REAL(wp_), PARAMETER :: me_ = 9.10938215d-31 ! [kg]
! REAL(wp_), PARAMETER :: mp_ = 1.672621637d-27 ! [kg]
real (kind = R8), parameter :: e_ = 1.602176487e-19_R8 ! [C]
real (kind = R8), parameter :: me_ = 9.10938215e-31_R8 ! electron mass [kg]
! real (kind = R8), parameter :: mp_ = 1.672621637e-27_R8 ! proton mass [kg]
! REAL(wp_), PARAMETER :: rmpe_ = mp_/me_
! REAL(wp_), PARAMETER :: c_ = 2.99792458d+08 ! [m/s]
! REAL(wp_), PARAMETER :: eps0_ = 8.854187817d-12 ! [F/m]
! real (kind = R8), parameter :: md_ = 3.34358320e-27_R8 ! deuteron mass, kg
! real (kind = R8), parameter :: mt_ = 5.00735588e-27_R8 ! triton mass, kg
! real (kind = R8), parameter :: ma_ = 6.64465620e-27_R8 ! alpha mass, kg
! real (kind = R8), parameter :: amu_ = 1.660538782e-27_R8 ! amu, kg
real (kind = R8), parameter :: c_ = 2.99792458e8_R8 ! speed of light [m/s]
real (kind = R8), parameter :: mu0_ = 4.0e-7_R8 * pi
! real (kind = R8), parameter :: eps0_ = 1.0_R8 / (mu0_ * c_ * c_) ! [F/m]
! real (kind = R8), parameter :: avogr_ = 6.02214179e23_R8
! real (kind = R8), parameter :: KBolt_ = 1.3806504e-23_R8
!------------------------------------------------------------------------
! Useful definitions
!------------------------------------------------------------------------
REAL(wp_), PARAMETER :: keV_ = 1000*e_ ! [J]
real (kind = R8), parameter :: ev_ = e_ ! [J]
REAL(wp_), PARAMETER :: keV_ = 1000*ev_ ! [J]
REAL(wp_), PARAMETER :: mc2_SI = me_*c_**2 ! [J]
REAL(wp_), PARAMETER :: mc2_ = mc2_SI/keV_ ! [keV]
! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s]

167
src/gray-externals.f
View File

@ -103,7 +103,7 @@ c
c
c
subroutine after_onestep(i,istop)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : psjki,ppabs,ccci,iiv,tauv,
. iop,iow,tau1v,yyrfl,nrayr,nrayth
implicit real*8 (a-h,o-z)
@ -284,7 +284,7 @@ c
c
c
subroutine print_output(i,j,k)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : ywrk,psjki,tauv,alphav,pdjki,
. currj,didst,q,tau1v,nrayr!,nrayth
implicit real*8 (a-h,o-z)
@ -447,7 +447,7 @@ c
. rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te,
. dne, zeff, qsf, powin)
use green_func_p, only:Setup_SpitzFunc
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : nrayr,nrayth,nstep
implicit real*8 (a-h,o-z)
integer ijetto, mr, mz, nrho, nbnd
@ -732,10 +732,12 @@ c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
anx0c=(anr0c*x00-anphi0c*y00)/r00
any0c=(anr0c*y00+anphi0c*x00)/r00
c
print*,' input file read'
! call myflush
c
c read data for Te , ne , Zeff from file if iprof>0
c
if (iprof.eq.1) then
c nprof=98
c lprf=length(filenmprf)
@ -746,6 +748,8 @@ c read profiles from input arguments
call profdata(nrho, psijet, te, dne, zeff)
c close(nprof)
end if
print*,' profiles fitted'
! call myflush
c
c read equilibrium data from file if iequil=2
c
@ -758,6 +762,8 @@ c . status= 'unknown', unit=neqdsk)
call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd,
. rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf)
c close(neqdsk)
print*,' equilibrium fitted'
! call myflush
c print density, temperature, safecty factor, toroidal current dens
c versus psi, rhop, rhot
@ -841,7 +847,7 @@ c
c
c
subroutine surf_anal
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
implicit real*8(a-h,o-z)
common/parban/b0,rr0m,zr0m,rpam
common/parbres/bres
@ -994,7 +1000,7 @@ c
c
subroutine equidata(ijetto,mr,mz,r,z,psin2d,psiaxis,psiedge,
. rax,zax,nbnd,rbnd,zbnd,mrho,psijet,fpjet,qjet)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use reflections, only : inside
implicit real*8 (a-h,o-z)
integer ijetto,mr,mz,nbnd,mrho
@ -1100,7 +1106,7 @@ c psi function
psia0=psiedge-psiaxis
psia=psia0*factb
sgniphi=sign(1.0d0,-psia0)
c
cc
c do j=1,nz
c do i=1,nr
c write(620,2021) rv(i),zv(j),psin(i,j)
@ -1120,15 +1126,15 @@ c psi(i,j)=psin(i,j)*psia
enddo
enddo
iopt=0
call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
call dierckx_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,
. kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
. wrk,lwrk,iwrk,liwrk,ier)
call dierckx_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)
end if
nsrt=nsr
nszt=nsz
@ -1178,10 +1184,10 @@ c if ier=-1 data are fitted using sspl=0
nsrt=nsr
nszt=nsz
end if
c
c re-evaluate psi on the grid using the spline (only for debug)
c
c call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
cc
cc re-evaluate psi on the grid using the spline (only for debug and cniteq)
cc
c call dierckx_bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
c
c do j=1,nz
@ -1195,31 +1201,31 @@ c write(619,*) ' '
c enddo
c
c2021 format(5(1x,e16.9))
c
nur=1
nuz=0
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier)
call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,
. nz,ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier)
c
nur=0
nuz=1
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier)
call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,
. nz,ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier)
c
nur=2
nuz=0
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier)
call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,
. nz,ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier)
c
nur=0
nuz=2
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier)
call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,
. nz,ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier)
c
nur=1
nuz=1
call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
. ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier)
call dierckx_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
@ -1235,10 +1241,10 @@ c
xb=0.0d0
xe=1.0d0
ssfp=0.01d0
call curfit(iopt,nrho,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft,
. tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier)
call dierckx_curfit(iopt,nrho,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,nrho,ier)
call dierckx_splev(tfp,nsft,cfp,3,psinr,fpoli,nrho,ier)
fpolas=fpoli(nrho)
c
c no limiter shape provided yet
@ -1406,6 +1412,7 @@ c
c
c
subroutine points_ox(rz,zz,rf,zf,psinvf,info)
use const_and_precisions, only : comp_eps
implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
@ -1413,7 +1420,7 @@ c
common/psival/psinv
xvec(1)=rz
xvec(2)=zz
tol = sqrt(dpmpar(1))
tol = sqrt(comp_eps)
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then
write(*,'(a,i2,a,2f8.4)') ' info subr points_ox =',info,
@ -1448,6 +1455,7 @@ c
c
c
subroutine points_tgo(rz,zz,rf,zf,psin,info)
use const_and_precisions, only : comp_eps
implicit real*8 (a-h,o-z)
parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
@ -1456,7 +1464,7 @@ c
h=psin
xvec(1)=rz
xvec(2)=zz
tol = sqrt(dpmpar(1))
tol = sqrt(comp_eps)
call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then
end if
@ -1611,7 +1619,7 @@ c
integer mrho
real*8 psijet(mrho), te(mrho), dne(mrho), zeff(mrho)
c
parameter(npmx=250,npest=npmx+4)
parameter(npmx=501,npest=npmx+4)
dimension psrad(npmx),terad(npmx),derad(npmx),zfc(npmx)
dimension ct(npmx,4),cz(npmx,4)
parameter(lwrkf=npmx*4+npest*16)
@ -1672,17 +1680,17 @@ c
kspl=3
sspl=.001d0
c
call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
. tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
call dierckx_curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,
. nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
c
call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
call dierckx_splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
nu=1
call gsplder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
call dierckx_splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
dnpp=densi(npp)
ddnpp=ddensi(npp)
c
nu=2
call gsplder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
call dierckx_splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
d2dnpp=d2densi(npp)
if(derad(npp).eq.0.0d0) then
@ -1823,10 +1831,10 @@ c spline interpolation of rhopol versus rhotor
xe=1.0d0
ss=0.00001
kspl=3
call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp,
. trp,crp,rp,wrkp,lwrkp,iwrkp,ier)
call dierckx_curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,
. nsrp,trp,crp,rp,wrkp,lwrkp,iwrkp,ier)
c write(*,*) ier
call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier)
call dierckx_splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier)
do i=1,nnr
write(644,*) rhop(i),rhot(i),rhopi(i)
end do
@ -1842,7 +1850,7 @@ c write(*,*) ier
common/coffrn/nsrp
common/coffrp/crp
rrs(1)=rhot
call splev(trp,nsrp,crp,3,rrs,ffspl,1,ier)
call dierckx_splev(trp,nsrp,crp,3,rrs,ffspl,1,ier)
frhopol=ffspl(1)
return
end
@ -2039,7 +2047,7 @@ c
c
c
subroutine contours_psi(h,np,rcn,zcn,ipr)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z)
parameter(mest=4,kspl=3)
parameter(nnw=501,nnh=501)
@ -2073,10 +2081,11 @@ c
do ic=2,np
zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0
iopt=1
call sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
if(ier.gt.0) write(*,*) ' sprofil =',ier
call dierckx_sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,
. ier)
if(ier.gt.0) write(*,*) ' dierckx_sprofil =',ier
val=h*psiant+psinop
call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
call dierckx_sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
if (zeroc(1).gt.rwallm) then
rcn(ic)=zeroc(1)
zcn(ic)=zc
@ -2104,13 +2113,13 @@ c
c
c
subroutine flux_average
use itm_constants, only : pi=>itm_pi, itm_mu0
use const_and_precisions, only : pi, mu0=>mu0_
implicit real*8 (a-h,o-z)
real*8 lam
parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41)
parameter(zero=0.0d0,one=1.0d0)
parameter(ccj=1.0d0/itm_mu0)
parameter(ccj=1.0d0/mu0)
parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
. njest*nnintp+nlest+54)
@ -2448,7 +2457,7 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
iopt=0
s=0.0d0
call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam,
call dierckx_regrid(iopt,nintp,rpstab,nlam,alam,ffhlam,
. zero,one,zero,one,ksp,ksp,s,
. njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier)
njpt=njp
@ -3115,7 +3124,7 @@ c
c
c
subroutine plas_deriv(xx,yy,zz)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z)
dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3)
dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3)
@ -3464,10 +3473,10 @@ c
c
if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then
rrs(1)=psinv
call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
call dierckx_splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
fpolv=ffspl(1)
nu=1
call gsplder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
call dierckx_splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
dfpolv=ffspl(1)
ffpv=fpolv*dfpolv/psia
end if
@ -3497,9 +3506,9 @@ c
c
subroutine tor_curr(rpsim,zpsim,ajphi)
use itm_constants, only : pi=>itm_pi, itm_mu0
use const_and_precisions, only : pi, mu0=>mu0_
implicit real*8 (a-h,o-z)
parameter(ccj=1.0d0/itm_mu0)
parameter(ccj=1.0d0/mu0)
common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
call equinum(rpsim,zpsim)
bzz= dpsidr/rpsim
@ -3525,10 +3534,11 @@ c
c
iopt=1
zc=zmaxis
call sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
if(ier.gt.0) write(*,*) ' sprofil =',ier
call dierckx_sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,
. ier)
if(ier.gt.0) write(*,*) ' dierckx_sprofil =',ier
val=h*psiant+psinop
call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
call dierckx_sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
r1=zeroc(1)
r2=zeroc(2)
return
@ -3558,7 +3568,7 @@ c
c
subroutine density(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250,npest=npmx+4)
parameter(npmx=501,npest=npmx+4)
dimension xxs(1),ffs(1)
dimension tfn(npest),cfn(npest),wrkfd(npest)
c
@ -3599,11 +3609,11 @@ c
else
xxs(1)=arg
ier=0
call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
call dierckx_splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
dens=ffs(1)
nu=1
ier=0
call gsplder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
call dierckx_splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
ddens=ffs(1)
if(ier.gt.0) write(*,*) ier
if(abs(dens).lt.1.0d-10) dens=0.0d0
@ -3618,7 +3628,7 @@ c
c
function temperature(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250)
parameter(npmx=501)
dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4)
c
common/parqte/te0,dte0,alt1,alt2
@ -3646,7 +3656,7 @@ c
c
function fzeff(arg)
implicit real*8 (a-h,o-z)
parameter(npmx=250)
parameter(npmx=501)
dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4)
c
common/iipr/iprof
@ -3673,7 +3683,7 @@ c
c beam tracing initial conditions igrad=1
c
subroutine ic_gb
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,
. gri,ggri
@ -3943,7 +3953,7 @@ c
c ray tracing initial conditions igrad=0
c
subroutine ic_rt
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,
. gri,ggri
@ -4088,7 +4098,7 @@ c
subroutine ic_rt2
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,grad2,dgrad2v,
. gri,ggri,yyrfl
@ -4266,7 +4276,7 @@ c
c
c
subroutine pabs_curr(i,j,k)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst,
. ccci,q,tau1v
implicit real*8 (a-h,o-z)
@ -5845,7 +5855,7 @@ c
c
c
subroutine pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
use beamdata, only : psjki,iiv,ppabs,ccci,pdjki,
. nrayr,nrayth,nstep
implicit real*8(a-h,o-z)
@ -5958,8 +5968,11 @@ c radial coordinate of i-(i+1) interval mid point
idecr=-1
is=1
do i=1,ii
if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i))
if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i)))
if(ipec.eq.0) then
xxi(i)=abs(psjki(j,k,i))
else
xxi(i)=sqrt(abs(psjki(j,k,i)))
end if
if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then
ypt(i)=ppabs(j,k,i)
yamp(i)=ccci(j,k,i)
@ -5983,15 +5996,16 @@ c radial coordinate of i-(i+1) interval mid point
else
if(xxi(i).gt.rtbc) exit
if(xxi(i).lt.xxi(i-1)) then
isev(is)=i !!!!!!!!!! it should be isev(is)=i-1
! isev(is)=i !!!!!!!!!! it should be isev(is)=i-1
isev(is)=i-1
is=is+1
idecr=-1
end if
end if
end if
end do
c
isev(is)=i-1
c
ppa1=0.0d0
cci1=0.0d0
do iis=1,is-1
@ -6010,11 +6024,12 @@ c
iind=1
end if
do ind=ind1,ind2,iind !!!!!!!!!! is it safe?
iind=ind !!!!!!!!!! iind reused in the loop!
if (idecr.eq.-1) iind=ind-1
rt1=rtab1(iind)
! iind=ind !!!!!!!!!! iind reused in the loop!
indi=ind !!!!!!!!!! iind reused in the loop!
if (idecr.eq.-1) indi=ind-1
rt1=rtab1(indi)
call locatex(xxi,iise,iise0,iise,rt1,itb1)
if(itb1.gt.iise0.and.itb1.lt.iise) then
if(itb1.ge.iise0.and.itb1.lt.iise) then
call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),
. ypt(itb1+1),rt1,ppa2)
call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),
@ -6370,7 +6385,7 @@ c
end
subroutine pol_limit(qq,uu,vv)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
implicit none
integer*4 ipolc
real*8 bv(3),anv(3)
@ -6468,7 +6483,7 @@ c
subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl)
use itm_constants, only : pi=>itm_pi
use const_and_precisions, only : pi
implicit none
integer*4 ivac,irfl
real*8 anv(3),xv(3),xvrfl(3)

14
src/gray_main.f90
View File

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

277
src/grayl.f
View File

@ -1660,11 +1660,11 @@ c
subroutine calcei(arg,result,int)
c----------------------------------------------------------------------
c
c this fortran 77 packet computes the exponential integrals ei(x),
c e1(x), and exp(-x)*ei(x) for real arguments x where
c this fortran 77 packet computes the exponential integrals eint(x),
c e1(x), and exp(-x)*eint(x) for real arguments x where
c
c integral (from t=-infinity to t=x) (exp(t)/t), x > 0,
c ei(x) =
c eint(x) =
c -integral (from t=-x to t=infinity) (exp(t)/t), x < 0,
c
c and where the first integral is a principal value integral.
@ -1672,14 +1672,14 @@ c the packet contains three function type subprograms: ei, eone,
c and expei; and one subroutine type subprogram: calcei. the
c calling statements for the primary entries are
c
c y = ei(x), where x .ne. 0,
c y = eint(x), where x .ne. 0,
c
c y = eone(x), where x .gt. 0,
c and
c y = expei(x), where x .ne. 0,
c
c and where the entry points correspond to the functions ei(x),
c e1(x), and exp(-x)*ei(x), respectively. the routine calcei
c and where the entry points correspond to the functions eint(x),
c e1(x), and exp(-x)*eint(x), respectively. the routine calcei
c is intended for internal packet use only, all computations within
c the packet being concentrated in this routine. the function
c subprograms invoke calcei with the fortran statement
@ -1689,9 +1689,9 @@ c
c function parameters for calcei
c call arg result int
c
c ei(x) x .ne. 0 ei(x) 1
c eone(x) x .gt. 0 -ei(-x) 2
c expei(x) x .ne. 0 exp(-x)*ei(x) 3
c eint(x) x .ne. 0 eint(x) 1
c eone(x) x .gt. 0 -eint(-x) 2
c expei(x) x .ne. 0 exp(-x)*eint(x) 3
c----------------------------------------------------------------------
integer i,int
double precision
@ -1936,11 +1936,11 @@ c----------------------------------------------------------------------
return
c---------- last line of calcei ----------
end
function ei(x)
function eint(x)
c--------------------------------------------------------------------
c
c this function program computes approximate values for the
c exponential integral ei(x), where x is real.
c exponential integral eint(x), where x is real.
c
c author: w. j. cody
c
@ -1948,11 +1948,11 @@ c latest modification: january 12, 1988
c
c--------------------------------------------------------------------
integer int
double precision ei, x, result
double precision eint, x, result
c--------------------------------------------------------------------
int = 1
call calcei(x,result,int)
ei = result
eint = result
return
c---------- last line of ei ----------
end
@ -1960,7 +1960,7 @@ c---------- last line of ei ----------
c--------------------------------------------------------------------
c
c this function program computes approximate values for the
c function exp(-x) * ei(x), where ei(x) is the exponential
c function exp(-x) * eint(x), where eint(x) is the exponential
c integral, and x is real.
c
c author: w. j. cody
@ -2004,11 +2004,11 @@ c
subroutine calcei3(arg,result)
c----------------------------------------------------------------------
c
c this fortran 77 packet computes the exponential integrals ei(x),
c e1(x), and exp(-x)*ei(x) for real arguments x where
c this fortran 77 packet computes the exponential integrals eint(x),
c e1(x), and exp(-x)*eint(x) for real arguments x where
c
c integral (from t=-infinity to t=x) (exp(t)/t), x > 0,
c ei(x) =
c eint(x) =
c -integral (from t=-x to t=infinity) (exp(t)/t), x < 0,
c
c and where the first integral is a principal value integral.
@ -2016,14 +2016,14 @@ c the packet contains three function type subprograms: ei, eone,
c and expei; and one subroutine type subprogram: calcei. the
c calling statements for the primary entries are
c
c y = ei(x), where x .ne. 0,
c y = eint(x), where x .ne. 0,
c
c y = eone(x), where x .gt. 0,
c and
c y = expei(x), where x .ne. 0,
c
c and where the entry points correspond to the functions ei(x),
c e1(x), and exp(-x)*ei(x), respectively. the routine calcei
c and where the entry points correspond to the functions eint(x),
c e1(x), and exp(-x)*eint(x), respectively. the routine calcei
c is intended for internal packet use only, all computations within
c the packet being concentrated in this routine. the function
c subprograms invoke calcei with the fortran statement
@ -2033,9 +2033,9 @@ c
c function parameters for calcei
c call arg result int
c
c ei(x) x .ne. 0 ei(x) 1
c eone(x) x .gt. 0 -ei(-x) 2
c expei(x) x .ne. 0 exp(-x)*ei(x) 3
c eint(x) x .ne. 0 eint(x) 1
c eone(x) x .gt. 0 -eint(-x) 2
c expei(x) x .ne. 0 exp(-x)*eint(x) 3
c----------------------------------------------------------------------
integer i,int
double precision
@ -4677,7 +4677,8 @@ c
c
c
c
subroutine surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
subroutine dierckx_surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,
* nxest,nyest,
* nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
c given the set of data points (x(i),y(i),z(i)) and the set of positive
c numbers w(i),i=1,...,m, subroutine surfit determines a smooth bivar-
@ -5072,8 +5073,8 @@ c we partition the working space and determine the spline approximation
lby = lbx+nek
lsx = lby+nek
lsy = lsx+m*km1
call fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
* eps,tol,maxit,nest,km1,km2,ib1,ib3,ncest,nrint,nreg,nx,tx,
call dierckx_fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,
* nyest,eps,tol,maxit,nest,km1,km2,ib1,ib3,ncest,nrint,nreg,nx,tx,
* ny,ty,c,fp,wrk1(1),wrk1(lfp),wrk1(lco),wrk1(lf),wrk1(lff),
* wrk1(la),wrk1(lq),wrk1(lbx),wrk1(lby),wrk1(lsx),wrk1(lsy),
* wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier)
@ -5081,8 +5082,8 @@ c we partition the working space and determine the spline approximation
end
subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest,
* nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest,
subroutine dierckx_fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,
* nxest,nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest,
* nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h,
* index,nummer,wrk,lwrk,ier)
c ..
@ -5106,7 +5107,7 @@ c ..local scalars..
c ..local arrays..
real*8 hx(6),hy(6)
c ..function references..
real*8 abs,fprati,sqrt
real*8 abs,dierckx_fprati,sqrt
integer min0
c ..subroutine references..
c fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota
@ -5242,7 +5243,7 @@ c a minimal bandwidth.
ky1 = ky+1
130 iband = iband1+1
c arrange the data points according to the panel they belong to.
call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg)
call dierckx_fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg)
c find ncof, the number of b-spline coefficients.
nk1x = nx-kx1
nk1y = ny-ky1
@ -5274,9 +5275,9 @@ c fetch a new data point.
wi = w(in)
zi = z(in)*wi
c evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in).
call fpbspl(tx,nx,kx,x(in),l1,hx)
call dierckx_fpbspl(tx,nx,kx,x(in),l1,hx)
c evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in).
call fpbspl(ty,ny,ky,y(in),l2,hy)
call dierckx_fpbspl(ty,ny,ky,y(in),l2,hy)
c store the value of these b-splines in spx and spy respectively.
do 160 i=1,kx1
spx(in,i) = hx(i)
@ -5307,16 +5308,16 @@ c rotate the row into triangle by givens transformations .
piv = h(i)
if(piv.eq.0.) go to 220
c calculate the parameters of the givens transformation.
call fpgivs(piv,a(irot,1),cos,sin)
call dierckx_fpgivs(piv,a(irot,1),cos,sin)
c apply that transformation to the right hand side.
call fprota(cos,sin,zi,f(irot))
call dierckx_fprota(cos,sin,zi,f(irot))
if(i.eq.iband) go to 230
c apply that transformation to the left hand side.
i2 = 1
i3 = i+1
do 210 j=i3,iband
i2 = i2+1
call fprota(cos,sin,h(j),a(irot,i2))
call dierckx_fprota(cos,sin,h(j),a(irot,i2))
210 continue
220 continue
c add the contribution of the row to the sum of squares of residual
@ -5339,7 +5340,7 @@ c check whether the observation matrix is rank deficient.
if(a(i,1).le.sigma) go to 280
270 continue
c backward substitution in case of full rank.
call fpback(a,f,ncof,iband,c,nc)
call dierckx_fpback(a,f,ncof,iband,c,nc)
rank = ncof
do 275 i=1,ncof
q(i,1) = a(i,1)/dmax
@ -5357,7 +5358,7 @@ c check whether there is sufficient working space
lf =1
lh = lf+ncof
la = lh+iband
call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la),
call dierckx_fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la),
* wrk(lf),wrk(lh))
do 295 i=1,ncof
q(i,1) = q(i,1)/dmax
@ -5489,13 +5490,13 @@ c test whether there are interior knots in the x-direction.
if(nk1x.eq.kx1) go to 440
c evaluate the discotinuity jumps of the kx-th order derivative of
c the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1.
call fpdisc(tx,nx,kx2,bx,nmax)
call dierckx_fpdisc(tx,nx,kx2,bx,nmax)
440 ky2 = ky1 + 1
c test whether there are interior knots in the y-direction.
if(nk1y.eq.ky1) go to 450
c evaluate the discontinuity jumps of the ky-th order derivative of
c the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1.
call fpdisc(ty,ny,ky2,by,nmax)
call dierckx_fpdisc(ty,ny,ky2,by,nmax)
c initial value for p.
450 p1 = 0.
f1 = fp0-s
@ -5549,14 +5550,14 @@ c square roots.
i2 = min0(iband1,ncof-irot)
if(piv.eq.0.) if(i2) 550,550,520
c calculate the parameters of the givens transformation.
call fpgivs(piv,q(irot,1),cos,sin)
call dierckx_fpgivs(piv,q(irot,1),cos,sin)
c apply that givens transformation to the right hand side.
call fprota(cos,sin,zi,ff(irot))
call dierckx_fprota(cos,sin,zi,ff(irot))
if(i2.eq.0) go to 550
c apply that givens transformation to the left hand side.
do 510 l=1,i2
l1 = l+1
call fprota(cos,sin,h(l1),q(irot,l1))
call dierckx_fprota(cos,sin,h(l1),q(irot,l1))
510 continue
520 do 530 l=1,i2
h(l) = h(l+1)
@ -5589,14 +5590,14 @@ c rotate the new row into triangle by givens transformations .
i2 = min0(iband3,ncof-irot)
if(piv.eq.0.) if(i2) 630,630,600
c calculate the parameters of the givens transformation.
call fpgivs(piv,q(irot,1),cos,sin)
call dierckx_fpgivs(piv,q(irot,1),cos,sin)
c apply that givens transformation to the right hand side.
call fprota(cos,sin,zi,ff(irot))
call dierckx_fprota(cos,sin,zi,ff(irot))
if(i2.eq.0) go to 630
c apply that givens transformation to the left hand side.
do 590 l=1,i2
l1 = l+1
call fprota(cos,sin,h(l1),q(irot,l1))
call dierckx_fprota(cos,sin,h(l1),q(irot,l1))
590 continue
600 do 610 l=1,i2
h(l) = h(l+1)
@ -5617,7 +5618,7 @@ c check whether the matrix is rank deficient.
if(q(i,1).le.sigma) go to 670
660 continue
c backward substitution in case of full rank.
call fpback(q,ff,ncof,iband4,c,nc)
call dierckx_fpback(q,ff,ncof,iband4,c,nc)
rank = ncof
go to 675
c in case of rank deficiency, find the minimum norm solution.
@ -5626,7 +5627,7 @@ c in case of rank deficiency, find the minimum norm solution.
lf = 1
lh = lf+ncof
la = lh+iband4
call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la),
call dierckx_fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la),
* wrk(lf),wrk(lh))
675 do 680 i=1,ncof
q(i,1) = q(i,1)/dmax
@ -5687,7 +5688,7 @@ c test whether the iteration process proceeds as theoretically
c expected.
760 if(f2.ge.f1 .or. f2.le.f3) go to 800
c find the new value of p.
p = fprati(p1,f1,p2,f2,p3,f3)
p = dierckx_fprati(p1,f1,p2,f2,p3,f3)
770 continue
c error codes and messages.
780 ier = lwest
@ -5748,7 +5749,7 @@ c if not, interchange x and y once more.
940 return
end
subroutine fprank(a,f,n,m,na,tol,c,sq,rank,aa,ff,h)
subroutine dierckx_fprank(a,f,n,m,na,tol,c,sq,rank,aa,ff,h)
c subroutine fprank finds the minimum norm solution of a least-
c squares problem in case of rank deficiency.
c
@ -5803,12 +5804,12 @@ c the rank deficiency is increased by one.
i2 = min0(n-ii,m1)
piv = h(1)
if(piv.eq.0.) go to 30
call fpgivs(piv,a(ii,1),cos,sin)
call fprota(cos,sin,yi,f(ii))
call dierckx_fpgivs(piv,a(ii,1),cos,sin)
call dierckx_fprota(cos,sin,yi,f(ii))
if(i2.eq.0) go to 70
do 20 j=1,i2
j1 = j+1
call fprota(cos,sin,h(j1),a(ii,j1))
call dierckx_fprota(cos,sin,h(j1),a(ii,j1))
h(j) = h(j1)
20 continue
go to 50
@ -5894,13 +5895,13 @@ c rotate this column into aa by givens transformations.
h(j2) = h(j3)
150 continue
go to 180
160 call fpgivs(piv,aa(jj,1),cos,sin)
160 call dierckx_fpgivs(piv,aa(jj,1),cos,sin)
if(j1.eq.0) go to 200
kk = jj
do 170 j2=1,j1
j3 = j2+1
kk = kk-1
call fprota(cos,sin,h(j3),aa(kk,j3))
call dierckx_fprota(cos,sin,h(j3),aa(kk,j3))
h(j2) = h(j3)
170 continue
180 jj = jj-1
@ -5984,7 +5985,8 @@ c to zero the small diagonal elements of matrix (a).
return
end
subroutine fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg)
subroutine dierckx_fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,
* nreg)
c subroutine fporde sorts the data points (x(i),y(i)),i=1,2,...,m
c according to the panel tx(l)<=x<tx(l+1),ty(k)<=y<ty(k+1), they belong
c to. for each panel a stack is constructed containing the numbers
@ -6033,8 +6035,8 @@ c ..
subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,
* iwrk,kwrk,ier)
subroutine dierckx_bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,
* lwrk,iwrk,kwrk,ier)
c subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
c ,my a bivariate spline s(x,y) of degrees kx and ky, given in the
c b-spline representation.
@ -6128,12 +6130,12 @@ c are invalid control is immediately repassed to the calling program.
50 continue
60 ier = 0
iw = mx*(kx+1)+1
call fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),wrk(iw),
* iwrk(1),iwrk(mx+1))
call dierckx_fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),
* wrk(iw),iwrk(1),iwrk(mx+1))
100 return
end
c
subroutine fpback(a,z,n,k,c,nest)
subroutine dierckx_fpback(a,z,n,k,c,nest)
c subroutine fpback calculates the solution of the system of
c equations a*c = z with a a n x n upper triangular matrix
c of bandwidth k.
@ -6166,7 +6168,8 @@ c ..
end
c
subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly)
subroutine dierckx_fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,
* lx,ly)
c ..scalar arguments..
integer nx,ny,kx,ky,mx,my
c ..array arguments..
@ -6195,7 +6198,7 @@ c ..
l = l1
l1 = l+1
go to 10
20 call fpbspl(tx,nx,kx,arg,l,h)
20 call dierckx_fpbspl(tx,nx,kx,arg,l,h)
lx(i) = l-kx1
do 30 j=1,kx1
wx(i,j) = h(j)
@ -6215,7 +6218,7 @@ c ..
l = l1
l1 = l+1
go to 50
60 call fpbspl(ty,ny,ky,arg,l,h)
60 call dierckx_fpbspl(ty,ny,ky,arg,l,h)
ly(i) = l-ky1
do 70 j=1,ky1
wy(i,j) = h(j)
@ -6245,7 +6248,7 @@ c ..
return
end
subroutine fpbspl(t,n,k,x,l,h)
subroutine dierckx_fpbspl(t,n,k,x,l,h)
c subroutine fpbspl evaluates the (k+1) non-zero b-splines of
c degree k at t(l) <= x < t(l+1) using the stable recurrence
c relation of de boor and cox.
@ -6279,7 +6282,7 @@ c ..
end
c
subroutine fpchec(x,m,t,n,k,ier)
subroutine dierckx_fpchec(x,m,t,n,k,ier)
c subroutine fpchec verifies the number and the position of the knots
c t(j),j=1,2,...,n of a spline of degree k, in relation to the number
c and the position of the data points x(i),i=1,2,...,m. if all of the
@ -6342,7 +6345,7 @@ c check condition no 5
80 return
end
c
subroutine fpdisc(t,n,k2,b,nest)
subroutine dierckx_fpdisc(t,n,k2,b,nest)
c subroutine fpdisc calculates the discontinuity jumps of the kth
c derivative of the b-splines of degree k at the knots t(k+2)..t(n-k-1)
c ..scalar arguments..
@ -6387,7 +6390,7 @@ c ..
end
c
subroutine fpgivs(piv,ww,cos,sin)
subroutine dierckx_fpgivs(piv,ww,cos,sin)
c subroutine fpgivs calculates the parameters of a givens
c transformation .
c ..
@ -6405,9 +6408,9 @@ c ..
ww = dd
return
end
subroutine fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,
* ty,ny,p,c,nc,fp,fpx,fpy,mm,mynx,kx1,kx2,ky1,ky2,spx,spy,right,q,
* ax,ay,bx,by,nrx,nry)
subroutine dierckx_fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,
* ky,tx,nx,ty,ny,p,c,nc,fp,fpx,fpy,mm,mynx,kx1,kx2,ky1,ky2,spx,spy,
* right,q,ax,ay,bx,by,nrx,nry)
c ..
c ..scalar arguments..
real*8 p,fp
@ -6470,7 +6473,7 @@ c ion problem in the x-direction.
l1 = l+1
number = number+1
go to 10
20 call fpbspl(tx,nx,kx,arg,l,h)
20 call dierckx_fpbspl(tx,nx,kx,arg,l,h)
do 30 i=1,kx1
spx(it,i) = h(i)
30 continue
@ -6491,7 +6494,7 @@ c ion problem in the y-direction.
l1 = l+1
number = number+1
go to 60
70 call fpbspl(ty,ny,ky,arg,l,h)
70 call dierckx_fpbspl(ty,ny,ky,arg,l,h)
do 80 i=1,ky1
spy(it,i) = h(i)
80 continue
@ -6501,11 +6504,11 @@ c ion problem in the y-direction.
100 if(p.le.0.0d0) go to 120
c calculate the non-zero elements of the matrix (bx).
if(ifbx.ne.0 .or. nx.eq.2*kx1) go to 110
call fpdisc(tx,nx,kx2,bx,nx)
call dierckx_fpdisc(tx,nx,kx2,bx,nx)
ifbx = 1
c calculate the non-zero elements of the matrix (by).
110 if(ifby.ne.0 .or. ny.eq.2*ky1) go to 120
call fpdisc(ty,ny,ky2,by,ny)
call dierckx_fpdisc(ty,ny,ky2,by,ny)
ifby = 1
c reduce the matrix (ax) to upper triangular form (rx) using givens
c rotations. apply the same transformations to the rows of matrix q
@ -6557,12 +6560,12 @@ c rotate the new row of matrix (ax) into triangle.
piv = h(i)
if(piv.eq.0.0d0) go to 240
c calculate the parameters of the givens transformation.
call fpgivs(piv,ax(irot,1),cos,sin)
call dierckx_fpgivs(piv,ax(irot,1),cos,sin)
c apply that transformation to the rows of matrix q.
iq = (irot-1)*my
do 220 j=1,my
iq = iq+1
call fprota(cos,sin,right(j),q(iq))
call dierckx_fprota(cos,sin,right(j),q(iq))
220 continue
c apply that transformation to the columns of (ax).
if(i.eq.ibandx) go to 250
@ -6570,7 +6573,7 @@ c apply that transformation to the columns of (ax).
i3 = i+1
do 230 j=i3,ibandx
i2 = i2+1
call fprota(cos,sin,h(j),ax(irot,i2))
call dierckx_fprota(cos,sin,h(j),ax(irot,i2))
230 continue
240 continue
250 if(nrold.eq.number) go to 270
@ -6627,11 +6630,11 @@ c rotate the new row of matrix (ay) into triangle.
piv = h(i)
if(piv.eq.0.0d0) go to 390
c calculate the parameters of the givens transformation.
call fpgivs(piv,ay(irot,1),cos,sin)
call dierckx_fpgivs(piv,ay(irot,1),cos,sin)
c apply that transformation to the colums of matrix g.
ic = irot
do 370 j=1,nk1x
call fprota(cos,sin,right(j),c(ic))
call dierckx_fprota(cos,sin,right(j),c(ic))
ic = ic+nk1y
370 continue
c apply that transformation to the columns of matrix (ay).
@ -6640,7 +6643,7 @@ c apply that transformation to the columns of matrix (ay).
i3 = i+1
do 380 j=i3,ibandy
i2 = i2+1
call fprota(cos,sin,h(j),ay(irot,i2))
call dierckx_fprota(cos,sin,h(j),ay(irot,i2))
380 continue
390 continue
400 if(nrold.eq.number) go to 420
@ -6652,7 +6655,7 @@ c solution of the linear system (ry) c (rx)' = h.
c first step: solve the system (ry) (c1) = h.
k = 1
do 450 i=1,nk1x
call fpback(ay,c(k),nk1y,ibandy,c(k),ny)
call dierckx_fpback(ay,c(k),nk1y,ibandy,c(k),ny)
k = k+nk1y
450 continue
c second step: solve the system c (rx)' = (c1).
@ -6664,7 +6667,7 @@ c second step: solve the system c (rx)' = (c1).
right(i) = c(l)
l = l+nk1y
460 continue
call fpback(ax,right,nk1x,ibandx,right,nx)
call dierckx_fpback(ax,right,nk1x,ibandx,right,nx)
l = k
do 470 i=1,nk1x
c(l) = right(i)
@ -6731,7 +6734,7 @@ c adjust the different parameters.
return
end
c
subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
subroutine dierckx_fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
c subroutine fpknot locates an additional knot for a spline of degree
c k and adjusts the corresponding parameters,i.e.
c t : the position of the knots.
@ -6797,7 +6800,7 @@ c adjust the different parameters.
end
c
subroutine fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,
subroutine dierckx_fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,
* nxest,nyest,tol,maxit,nc,nx,tx,ny,ty,c,fp,fp0,fpold,reducx,
* reducy,fpintx,fpinty,lastdi,nplusx,nplusy,nrx,nry,nrdatx,nrdaty,
* wrk,lwrk,ier)
@ -6817,7 +6820,7 @@ c ..local scalars
* nk1x,nk1y,nmaxx,nmaxy,nminx,nminy,nplx,nply,npl1,nrintx,
* nrinty,nxe,nxk,nye
c
real*8 fprati
real*8 dierckx_fprati
c ..subroutine references..
c fpgrre,fpknot
@ -7006,10 +7009,10 @@ c of squared residuals fpintx(j),j=1,2,...,nx-2*kx-1 (fpinty(j),j=1,2,
c ...,ny-2*ky-1) for the data points having their absciss (ordinate)-
c value belonging to that interval.
c fp gives the total sum of squared residuals.
call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty,
* ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx),
* wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby),
* nrx,nry)
call dierckx_fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,
* nx,ty,ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,
* wrk(lsx),wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),
* wrk(lby),nrx,nry)
if(ier.eq.(-2)) fp0 = fp
c test whether the least-squares spline is an acceptable solution.
if(iopt.lt.0) go to 440
@ -7055,7 +7058,7 @@ c addition in the x-direction.
ifsx = 0
do 220 l=1,nplusx
c add a new knot in the x-direction
call fpknot(x,mx,tx,nx,fpintx,nrdatx,nrintx,nxest,1)
call dierckx_fpknot(x,mx,tx,nx,fpintx,nrdatx,nrintx,nxest,1)
c test whether we cannot further increase the number of knots in the
c x-direction.
if(nx.eq.nxe) go to 250
@ -7068,7 +7071,7 @@ c addition in the y-direction.
ifsy = 0
do 240 l=1,nplusy
c add a new knot in the y-direction.
call fpknot(y,my,ty,ny,fpinty,nrdaty,nrinty,nyest,1)
call dierckx_fpknot(y,my,ty,ny,fpinty,nrdaty,nrinty,nyest,1)
c test whether we cannot further increase the number of knots in the
c y-direction.
if(ny.eq.nye) go to 250
@ -7107,10 +7110,10 @@ c iteration process to find the root of f(p)=s.
do 350 iter = 1,maxit
c find the smoothing spline sp(x,y) and the corresponding sum of
c squared residuals fp.
call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty,
* ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx),
* wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby),
* nrx,nry)
call dierckx_fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,
* nx,ty,ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,
* wrk(lsx),wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),
* wrk(lby),nrx,nry)
c test whether the approximation sp(x,y) is an acceptable solution.
fpms = fp-s
if(abs(fpms).lt.acc) go to 440
@ -7143,7 +7146,7 @@ c expected.
330 if(f2.gt.0.0d0) ich1 = 1
340 if(f2.ge.f1 .or. f2.le.f3) go to 410
c find the new value of p.
p = fprati(p1,f1,p2,f2,p3,f3)
p = dierckx_fprati(p1,f1,p2,f2,p3,f3)
350 continue
c error codes and messages.
400 ier = 3
@ -7157,7 +7160,7 @@ c error codes and messages.
440 return
end
c
subroutine fprota(cos,sin,a,b)
subroutine dierckx_fprota(cos,sin,a,b)
c subroutine fprota applies a givens rotation to a and b.
c ..
c ..scalar arguments..
@ -7174,7 +7177,7 @@ c ..
c
c
c
double precision function fprati(p1,f1,p2,f2,p3,f3)
double precision function dierckx_fprati(p1,f1,p2,f2,p3,f3)
c given three points (p1,f1),(p2,f2) and (p3,f3), function fprati
c gives the value of p such that the rational interpolating function
c of the form r(p) = (u*p+v)/(p+w) equals zero at p.
@ -7200,11 +7203,11 @@ c adjust the value of p1,f1,p3 and f3 such that f1 > 0 and f3 < 0.
go to 40
30 p3 = p2
f3 = f2
40 fprati = p
40 dierckx_fprati = p
return
end
c
subroutine regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,
subroutine dierckx_regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,
* nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
c given the set of values z(i,j) on the rectangular grid (x(i),y(j)),
c i=1,...,mx;j=1,...,my, subroutine regrid determines a smooth bivar-
@ -7523,7 +7526,7 @@ c are invalid, control is immediately repassed to the calling program.
tx(j) = xe
j = j-1
30 continue
call fpchec(x,mx,tx,nx,kx,ier)
call dierckx_fpchec(x,mx,tx,nx,kx,ier)
if(ier.ne.0) go to 70
if(ny.lt.nminy .or. ny.gt.nyest) go to 70
j = ny
@ -7532,7 +7535,7 @@ c are invalid, control is immediately repassed to the calling program.
ty(j) = ye
j = j-1
40 continue
call fpchec(y,my,ty,ny,ky,ier)
call dierckx_fpchec(y,my,ty,ny,ky,ier)
if(ier) 70,60,70
50 if(s.lt.0.0d0) go to 70
if(s.eq.0.0d0 .and. (nxest.lt.(mx+kx1) .or. nyest.lt.(my+ky1)) )
@ -7547,14 +7550,14 @@ c we partition the working space and determine the spline approximation
knry = knrx+mx
kndx = knry+my
kndy = kndx+nxest
call fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
* tol,maxit,nc,nx,tx,ny,ty,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),
call dierckx_fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,nxest,
* nyest,tol,maxit,nc,nx,tx,ny,ty,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),
* wrk(lfpx),wrk(lfpy),iwrk(1),iwrk(2),iwrk(3),iwrk(knrx),
* iwrk(knry),iwrk(kndx),iwrk(kndy),wrk(lww),jwrk,ier)
70 return
end
c
subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,
subroutine dierckx_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,
* wrk,lwrk,iwrk,kwrk,ier)
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
@ -7724,13 +7727,13 @@ c we calculate the b-spline coefficients of this spline
c we partition the working space and evaluate the partial derivative
300 iwx = 1+nxx*nyy
iwy = iwx+mx*(kx1-nux)
call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky,
* x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1))
call dierckx_fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,
* kky,x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1))
400 return
end
subroutine coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,
subroutine dierckx_coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,
* wrk,lwrk,ier)
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
@ -7882,7 +7885,7 @@ c we calculate the b-spline coefficients of this spline
c
c
c
subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,
subroutine dierckx_curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,
* wrk,lwrk,iwrk,ier)
c given the set of data points (x(i),y(i)) and the set of positive
c numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
@ -8126,7 +8129,7 @@ c are invalid, control is immediately repassed to the calling program.
t(j) = xe
j = j-1
20 continue
call fpchec(x,m,t,n,k,ier)
call dierckx_fpchec(x,m,t,n,k,ier)
if(ier) 50,40,50
30 if(s.lt.0.0d0) go to 50
if(s.eq.0.0d0 .and. nest.lt.(m+k1)) go to 50
@ -8138,15 +8141,15 @@ c we partition the working space and determine the spline approximation.
ib = ia+nest*k1
ig = ib+nest*k2
iq = ig+nest*k2
call fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,
* wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
call dierckx_fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,
* t,c,fp,wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
50 return
end
c
c
c
subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,
* n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
subroutine dierckx_fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,
* k1,k2,n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
c ..
c ..scalar arguments..
real*8 xb,xe,s,tol,fp
@ -8163,7 +8166,7 @@ c ..local scalars..
c ..local arrays..
real*8 h(7)
c ..function references
real*8 abs,fprati
real*8 abs,dierckx_fprati
integer max0,min0
c ..subroutine references..
c fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
@ -8274,7 +8277,7 @@ c search for knot interval t(l) <= xi < t(l+1).
l = l+1
go to 85
c evaluate the (k+1) non-zero b-splines at xi and store them in q.
90 call fpbspl(t,n,k,xi,l,h)
90 call dierckx_fpbspl(t,n,k,xi,l,h)
do 95 i=1,k1
q(it,i) = h(i)
h(i) = h(i)*wi
@ -8286,16 +8289,16 @@ c rotate the new row of the observation matrix into triangle.
piv = h(i)
if(piv.eq.0.0d0) go to 110
c calculate the parameters of the givens transformation.
call fpgivs(piv,a(j,1),cos,sin)
call dierckx_fpgivs(piv,a(j,1),cos,sin)
c transformations to right hand side.
call fprota(cos,sin,yi,z(j))
call dierckx_fprota(cos,sin,yi,z(j))
if(i.eq.k1) go to 120
i2 = 1
i3 = i+1
do 100 i1 = i3,k1
i2 = i2+1
c transformations to left hand side.
call fprota(cos,sin,h(i1),a(j,i2))
call dierckx_fprota(cos,sin,h(i1),a(j,i2))
100 continue
110 continue
c add contribution of this row to the sum of squares of residual
@ -8307,7 +8310,7 @@ c right hand sides.
fpint(n-1) = fpold
nrdata(n) = nplus
c backward substitution to obtain the b-spline coefficients.
call fpback(a,z,nk1,k1,c,nest)
call dierckx_fpback(a,z,nk1,k1,c,nest)
c test whether the approximation sinf(x) is an acceptable solution.
if(iopt.lt.0) go to 440
fpms = fp-s
@ -8358,7 +8361,7 @@ c t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
fpint(nrint) = fpart
do 190 l=1,nplus
c add a new knot.
call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
call dierckx_fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
c if n=nmax we locate the knots as for interpolation.
if(n.eq.nmax) go to 10
c test whether we cannot further increase the number of knots.
@ -8392,7 +8395,7 @@ c guaranteed by taking f1>0 and f3<0. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c evaluate the discontinuity jump of the kth derivative of the
c b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
call fpdisc(t,n,k2,b,nest)
call dierckx_fpdisc(t,n,k2,b,nest)
c initial value for p.
p1 = 0.0d0
f1 = fp0-s
@ -8427,23 +8430,23 @@ c the row of matrix b is rotated into triangle by givens transformation
do 290 j=it,nk1
piv = h(1)
c calculate the parameters of the givens transformation.
call fpgivs(piv,g(j,1),cos,sin)
call dierckx_fpgivs(piv,g(j,1),cos,sin)
c transformations to right hand side.
call fprota(cos,sin,yi,c(j))
call dierckx_fprota(cos,sin,yi,c(j))
if(j.eq.nk1) go to 300
i2 = k1
if(j.gt.n8) i2 = nk1-j
do 280 i=1,i2
c transformations to left hand side.
i1 = i+1
call fprota(cos,sin,h(i1),g(j,i1))
call dierckx_fprota(cos,sin,h(i1),g(j,i1))
h(i) = h(i1)
280 continue
h(i2+1) = 0.0d0
290 continue
300 continue
c backward substitution to obtain the b-spline coefficients.
call fpback(g,c,nk1,k2,c,nest)
call dierckx_fpback(g,c,nk1,k2,c,nest)
c computation of f(p).
fp = 0.0d0
l = k2
@ -8489,7 +8492,7 @@ c test whether the iteration process proceeds as theoretically
c expected.
350 if(f2.ge.f1 .or. f2.le.f3) go to 410
c find the new value for p.
p = fprati(p1,f1,p2,f2,p3,f3)
p = dierckx_fprati(p1,f1,p2,f2,p3,f3)
360 continue
c error codes and messages.
400 ier = 3
@ -8504,7 +8507,7 @@ c error codes and messages.
c
c
c
subroutine gsplder(t,n,c,k,nu,x,y,m,wrk,ier)
subroutine dierckx_splder(t,n,c,k,nu,x,y,m,wrk,ier)
c subroutine gsplder evaluates in a number of points x(i),i=1,2,...,m
c the derivative of order nu of a spline s(x) of degree k,given in
c its b-spline representation.
@ -8629,7 +8632,7 @@ c search for knot interval t(l) <= arg < t(l+1)
l1 = l+1
go to 140
c evaluate the non-zero b-splines of degree k-nu at arg.
150 call fpbspl(t,n,kk,arg,l,h)
150 call dierckx_fpbspl(t,n,kk,arg,l,h)
c find the value of the derivative at x=arg.
sp = 0.0d0
ll = l-k1
@ -8644,7 +8647,7 @@ c find the value of the derivative at x=arg.
c
c
c
subroutine splev(t,n,c,k,x,y,m,ier)
subroutine dierckx_splev(t,n,c,k,x,y,m,ier)
c subroutine splev evaluates in a number of points x(i),i=1,2,...,m
c a spline s(x) of degree k, given in its b-spline representation.
c
@ -8727,7 +8730,7 @@ c search for knot interval t(l) <= arg < t(l+1)
l1 = l+1
go to 40
c evaluate the non-zero b-splines at arg.
50 call fpbspl(t,n,k,arg,l,h)
50 call dierckx_fpbspl(t,n,k,arg,l,h)
c find the value of s(x) at x=arg.
sp = 0.0d0
ll = l-k1
@ -8742,7 +8745,7 @@ c find the value of s(x) at x=arg.
c
c
c
subroutine sproota(val,t,n,c,zero,mest,m,ier)
subroutine dierckx_sproota(val,t,n,c,zero,mest,m,ier)
c subroutine sproot finds the zeros of a cubic spline s(x),which is
c given in its normalized b-spline representation.
c
@ -8885,7 +8888,7 @@ c t(l) <= x <= t(l+1).
* z3.and.z4).or.nz0.and.(z1.and.(nz3.or.nz2.and.z4).or.z2.and.
* nz3.and.nz4))))go to 200
c find the zeros of ql(y).
100 call fpcuro(a3,a2,a1,a0,y,j)
100 call dierckx_fpcuro(a3,a2,a1,a0,y,j)
if(j.eq.0) go to 200
c find which zeros of pl(x) are zeros of s(x).
do 150 i=1,j
@ -8926,7 +8929,7 @@ c the zeros of s(x) are arranged in increasing order.
end
c
subroutine sprofil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
subroutine dierckx_sprofil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
c if iopt=0 subroutine sprofil calculates the b-spline coefficients of
c the univariate spline f(y) = s(u,y) with s(x,y) a bivariate spline of
c degrees kx and ky, given in the b-spline representation.
@ -9005,7 +9008,7 @@ c the b-splinecoefficients of f(y) = s(u,y).
l = l1
l1 = l+1
go to 110
120 call fpbspl(tx,nx,kx,u,l,h)
120 call dierckx_fpbspl(tx,nx,kx,u,l,h)
m0 = (l-kx1)*nky1+1
do 140 i=1,nky1
m = m0
@ -9028,7 +9031,7 @@ c the b-splinecoefficients of g(x) = s(x,u).
l = l1
l1 = l+1
go to 210
220 call fpbspl(ty,ny,ky,u,l,h)
220 call dierckx_fpbspl(ty,ny,ky,u,l,h)
m0 = l-ky
do 240 i=1,nkx1
m = m0
@ -9043,7 +9046,7 @@ c the b-splinecoefficients of g(x) = s(x,u).
300 return
end
c
subroutine fpcuro(a,b,c,d,x,n)
subroutine dierckx_fpcuro(a,b,c,d,x,n)
c subroutine fpcuro finds the real zeros of a cubic polynomial
c p(x) = a*x**3+b*x**2+c*x+d.
c

32
src/itm_constants.f90
View File

@ -1,32 +0,0 @@
!> Module implementing the ITM physics constants
!>
!> Source:
!> based on SOLPS b2mod_constants.F
!> '09/12/07 xpb : source CODATA 2006 (http://www.nist.gov/)'
!> pulled from ets r100
!>
!> \author David Coster
!>
!> \version "$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $"
module itm_constants
use itm_types
real (kind = R8), parameter :: itm_pi = 3.141592653589793238462643383280_R8
real (kind = R8), parameter :: itm_c = 2.99792458e8_R8 ! speed of light, m/s
real (kind = R8), parameter :: itm_me = 9.10938215e-31_R8 ! electron mass, kg
real (kind = R8), parameter :: itm_mp = 1.672621637e-27_R8 ! proton mass, kg
real (kind = R8), parameter :: itm_md = 3.34358320e-27_R8 ! deuteron mass, kg
real (kind = R8), parameter :: itm_mt = 5.00735588e-27_R8 ! triton mass, kg
real (kind = R8), parameter :: itm_ma = 6.64465620e-27_R8 ! alpha mass, kg
real (kind = R8), parameter :: itm_amu = 1.660538782e-27_R8 ! amu, kg
real (kind = R8), parameter :: itm_ev = 1.602176487e-19_R8
real (kind = R8), parameter :: itm_qe = itm_ev
real (kind = R8), parameter :: itm_mu0 = 4.0e-7_R8 * itm_pi
real (kind = R8), parameter :: itm_eps0 = 1.0_R8 / (itm_mu0 * itm_c * itm_c)
real (kind = R8), parameter :: itm_avogr = 6.02214179e23_R8
real (kind = R8), parameter :: itm_KBolt = 1.3806504e-23_R8
character (len=64), parameter :: itm_constants_version = '$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $'
end module itm_constants

50
src/itm_types.f90
View File

@ -1,50 +0,0 @@
!> Module implementing the ITM basic types
!>
!> Source:
!> based on SOLPS b2mod_types.F
!> pulled from ets r100 and extended with input from C. Konz, T. Ribeiro & B. Scott
!>
!> \author David Coster
!>
!> \version "$Id: itm_types.f90 144 2010-10-07 09:26:24Z konz $"
module itm_types
INTEGER, PARAMETER :: ITM_I1 = SELECTED_INT_KIND (2) ! Integer*1
INTEGER, PARAMETER :: ITM_I2 = SELECTED_INT_KIND (4) ! Integer*2
INTEGER, PARAMETER :: ITM_I4 = SELECTED_INT_KIND (9) ! Integer*4
INTEGER, PARAMETER :: ITM_I8 = SELECTED_INT_KIND (18) ! Integer*8
INTEGER, PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37) ! Real*4
INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8
INTEGER, PARAMETER :: itm_int_invalid = -999999999
REAL(R8), PARAMETER :: itm_r8_invalid = -9.0D40
interface itm_is_valid
module procedure itm_is_valid_int4, itm_is_valid_int8, itm_is_valid_real8
end interface
contains
logical function itm_is_valid_int4(in_int)
implicit none
integer(ITM_I4) in_int
itm_is_valid_int4 = in_int .ne. itm_int_invalid
return
end function itm_is_valid_int4
logical function itm_is_valid_int8(in_int)
implicit none
integer(ITM_I8) in_int
itm_is_valid_int8 = in_int .ne. itm_int_invalid
return
end function itm_is_valid_int8
logical function itm_is_valid_real8(in_real)
implicit none
real(R8) in_real
itm_is_valid_real8 = abs(in_real - itm_r8_invalid) .gt. abs(itm_r8_invalid) * 1.0e-15_R8
return
end function itm_is_valid_real8
end module itm_types

286
src/main.f90 Normal file
View File

@ -0,0 +1,286 @@
program main
use const_and_precisions, only : r8, pi
implicit none
integer, parameter :: cocosout=3
integer :: u, ierr, cocos, exp2pi, exp2piout, idum, i, j, ipsinorm
logical :: ispsinorm
logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout
character(len=255) :: desc, eqdskfile, prffile
character(len=*), parameter :: fmt2000='(a48,3i4)',fmt2020='(5e16.9)',fmt2022='(2i5)'
integer :: nr, nz, nbnd, nprf
real(r8) :: deltar, deltaz, r0, r1, zmid, rax, zax, psiax, psib, bref, ipla, &
dummy, z1, dr, dz, dpsin, bref0, ipla0, qb0, psifact, psii, powin, pec, icd
real(r8), dimension(:,:), allocatable :: psi2d
real(r8), dimension(:), allocatable :: r, z, psi1d, q, fdia, rbnd, zbnd, &
psiprf, te, dne, zeff, cte, cne, czeff, teeq, dneeq, zeffeq, dpdv, jcd
real(r8), external :: fspli
! === read filenames and input power from gray.data ===
u = get_free_unit()
open(u,file='gray.data',status= 'unknown')
read(u,*)
read(u,*)
read(u,*) powin
do i=1,13
read(u,*)
end do
read(u,*) desc
eqdskfile=trim(desc)//'.eqdsk'
read(u,*) ipsinorm
ispsinorm=(ipsinorm==1)
read(u,*)
read(u,*) desc
prffile=trim(desc)//'.prf'
read(u,*)
read(u,*) dummy,dummy,cocos
close(u)
! === read EQDSK ===
u = get_free_unit()
open(unit=u,file=eqdskfile,status='OLD',action='READ',iostat=ierr)
if(ierr/=0) then
write(*,*) 'Cannot open file '//trim(eqdskfile)
stop
end if
! read grid size and allocate arrays
read (u,fmt2000) desc,idum,nr,nz
allocate(psi2d(nr,nz), r(nr), z(nz), psi1d(nr), q(nr), fdia(nr), stat=ierr)
if (ierr/=0) then
close(u)
call free_allocs
write(*,*) 'cannot allocate arrays for equilibrium data'
stop
end if
! read 0D fields
read (u,fmt2020) deltar, deltaz, r0, r1, zmid
read (u,fmt2020) rax, zax, psiax, psib, bref
read (u,fmt2020) ipla, (dummy,i=1,4)
read (u,fmt2020) (dummy,i=1,5)
! read 1D-2D fields
read (u,fmt2020) (fdia(i),i=1,nr)
read (u,fmt2020) (dummy,i=1,nr)
read (u,fmt2020) (dummy,i=1,nr)
read (u,fmt2020) (dummy,i=1,nr)
read (u,fmt2020) ((psi2d(i,j),i=1,nr),j=1,nz)
read (u,fmt2020) (q(i),i=1,nr)
! read boundary size and allocate arrays
read (u,fmt2022) nbnd
if (nbnd>0) then
allocate(rbnd(nbnd), zbnd(nbnd), stat=ierr)
if (ierr/=0) then
close(u)
call free_allocs
write(*,*) 'cannot allocate arrays for boundary data'
stop
end if
! read boundary shape
read (u,fmt2020) (rbnd(i), zbnd(i),i=1,nbnd)
end if
close(u)
! normalize psi2d
if (.not.ispsinorm) psi2d = (psi2d - psiax)/(psib - psiax)
! interpret cocos numbers
call decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout)
! check sign consistency
ipla0=ipla
bref0=bref
if (psiincr) then
ipla=sign(ipla,psib - psiax)
else
ipla=sign(ipla,psiax - psib)
end if
bref=sign(bref,fdia(nr))
if (ipla/ipla0<0 .or. bref/bref0<0) then
write(*,*) 'Warning: sign inconsistency in Ipla/psi or Bref/Fdia'
end if
qb0=q(nr)
if (qpos) then
q=sign(q,ipla*bref)
else
q=sign(q,-ipla*bref)
end if
if (q(nr)/qb0<0) then
write(*,*) 'Warning: sign inconsistency found among q, Ipla and Bref'
end if
! convert cocosin to cocosout
if (phiccw.neqv.phiccwout) then
! opposite direction of toroidal angle phi in cocosin and cocosout
bref=-bref
ipla=-ipla
fdia=-fdia
end if
if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) then
! psi and Ip signs don't change accordingly
psib=-psib
psiax=-psiax
end if
if (qpos .neqv. qposout) q=-q
! q has opposite sign for given sign of Bphi*Ip
if (exp2pi/=exp2piout) then
! convert Wb to Wb/rad or viceversa
psifact=(2._r8*pi)**(exp2piout-exp2pi)
psib=psib*psifact
psiax=psiax*psifact
end if
! fill equi-spaced R, z, psi arrays
z1 = zmid - deltaz/2._r8
dr = deltar/(nr-1)
dz = deltaz/(nz-1)
dpsin = 1._r8/(nr-1)
do i=1,nr
r(i) = r1 + (i-1)*dr
psi1d(i) = (i-1)*dpsin
end do
do j=1,nz
z(j) = z1 + (j-1)*dz
end do
! === read profiles ===
u = get_free_unit()
open(unit=u,file=prffile,status='OLD',action='READ',iostat=ierr)
if(ierr/=0) then
write(*,*) 'Cannot open file '//trim(prffile)
stop
end if
read(u,*) nprf
if (nprf>0) then
allocate(psiprf(nprf),te(nprf),dne(nprf),zeff(nprf), stat=ierr)
if (ierr==0) allocate(cte(4*nprf),cne(4*nprf),czeff(4*nprf), stat=ierr)
if (ierr/=0) then
close(u)
call free_allocs
write(*,*) 'cannot allocate arrays for input 1d profiles'
stop
end if
do i=1,nprf
read(u,*) psiprf(i),te(i),dne(i),zeff(i)
end do
else
write(*,*) 'no data for 1d profiles'
stop
end if
close(u)
! spline interpolation for resampling on uniform grid
call difcsg(psiprf,te ,nprf,0,cte ,ierr)
if (ierr==0) call difcsg(psiprf,dne ,nprf,0,cne ,ierr)
if (ierr==0) call difcsg(psiprf,zeff,nprf,0,czeff,ierr)
if (ierr/=0) then
call free_allocs
write(*,*) 'error in input 1d profiles interpolation'
stop
end if
allocate(teeq(nr),dneeq(nr),zeffeq(nr),dpdv(nr),jcd(nr), stat=ierr)
if (ierr/=0) then
call free_allocs
write(*,*) 'cannot allocate arrays for resampled and output 1d profiles'
stop
end if
do i=1,nr
psii=psi1d(i)
call vlocate(psiprf,nprf,psii,j)
j=max(1,min(j,nprf-1))
dpsin=psii-psiprf(j)
teeq(i) =fspli(cte, nprf,j,dpsin)
dneeq(i) =fspli(cne, nprf,j,dpsin)
zeffeq(i)=fspli(czeff,nprf,j,dpsin)
end do
! convert keV to eV, 10^19 m^-3 to m^-3, and MW to W
teeq=teeq*1.e3_r8
dneeq=dneeq*1.e19_r8
powin=powin*1.e6_r8
! === call GRAY subroutine ===
call gray_main(1, nr, nz, r, z, psi2d, psiax, psib, &
rax, zax, nbnd, rbnd, zbnd, nr, psi1d, fdia, teeq, dneeq, zeffeq, q, &
powin, dpdv, jcd, pec, icd, ierr)
call free_allocs
contains
function get_free_unit(umin,umax) result(i)
implicit none
integer :: i
integer, intent(in), optional :: umin, umax
integer, parameter :: max_allowed = 99
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) then
i=-2 ! cannot inquire unit i
exit
end if
if (ex.and..not.op) exit ! unit i exists and is not open
i = i + 1
end do
end function get_free_unit
subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
implicit none
integer, intent(in) :: cocos
integer, intent(out) :: exp2pi
logical, intent(out) :: phiccw,psiincr,qpos
integer :: cocosm10,cocosm4
cocosm10=mod(cocos,10)
cocosm4=mod(cocosm10,4)
! cocos>10 psi in Wb, cocos<10 psi in Wb/rad
exp2pi=(cocos-cocosm10)/10
! cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW
phiccw=(mod(cocosm10,2)==1)
! cocos mod 10 = 1,2,5,6: psi increasing with positive Ip
psiincr=(cocosm4==1 .or. cocosm4==2)
! cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip
qpos=(cocosm10<3 .or. cocosm10>6)
end subroutine decode_cocos
subroutine free_allocs
implicit none
if(allocated(psi2d)) deallocate(psi2d)
if(allocated(r)) deallocate(r)
if(allocated(z)) deallocate(z)
if(allocated(psi1d)) deallocate(psi1d)
if(allocated(q)) deallocate(q)
if(allocated(fdia)) deallocate(fdia)
if(allocated(rbnd)) deallocate(rbnd)
if(allocated(zbnd)) deallocate(zbnd)
if(allocated(psiprf)) deallocate(psiprf)
if(allocated(te)) deallocate(te)
if(allocated(dne)) deallocate(dne)
if(allocated(zeff)) deallocate(zeff)
if(allocated(cte)) deallocate(cte)
if(allocated(cne)) deallocate(cne)
if(allocated(czeff)) deallocate(czeff)
if(allocated(teeq)) deallocate(teeq)
if(allocated(dneeq)) deallocate(dneeq)
if(allocated(zeffeq)) deallocate(zeffeq)
if(allocated(dpdv)) deallocate(dpdv)
if(allocated(jcd)) deallocate(jcd)
end subroutine free_allocs
end program main

7
src/scatterspl.f90
View File

@ -1,6 +1,6 @@
subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
tx,nknt_x,ty,nknt_y,coeff,ierr)
use itm_types, only : r8
use const_and_precisions, only : r8, comp_eps
implicit none
integer :: n
real(r8), dimension(n) :: x, y, z
@ -15,7 +15,6 @@ subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
integer :: ierr
integer :: iopt
real(r8), parameter :: eps_r8=epsilon(1._r8)
real(r8) :: resid
integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest
real(r8), dimension(:), allocatable :: tx_tmp,ty_tmp
@ -38,8 +37,8 @@ subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk))
iopt=0
call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl,sspl, &
nxest,nyest,ne,eps_r8,nknt_x,tx_tmp,nknt_y,ty_tmp, &
call dierckx_surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, &
sspl,nxest,nyest,ne,comp_eps,nknt_x,tx_tmp,nknt_y,ty_tmp, &
coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr)
tx(1:nknt_x)=tx_tmp(1:nknt_x)
ty(1:nknt_y)=ty_tmp(1:nknt_y)