From 074f3313550b997257739ed58129cb045bd89032 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Mon, 16 Jun 2014 16:41:43 +0000 Subject: [PATCH] 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. --- Makefile | 12 +- Makefile.standalone | 16 +- src/beamdata.f90 | 72 +++++++++ src/const_and_precisions.f90 | 40 +++-- src/gray-externals.f | 167 ++++++++++---------- src/gray_main.f90 | 14 +- src/grayl.f | 277 ++++++++++++++++----------------- src/itm_constants.f90 | 32 ---- src/itm_types.f90 | 50 ------ src/main.f90 | 286 +++++++++++++++++++++++++++++++++++ src/scatterspl.f90 | 7 +- 11 files changed, 643 insertions(+), 330 deletions(-) create mode 100644 src/beamdata.f90 delete mode 100644 src/itm_constants.f90 delete mode 100644 src/itm_types.f90 create mode 100644 src/main.f90 diff --git a/Makefile b/Makefile index 4e3ba22..97f3803 100644 --- a/Makefile +++ b/Makefile @@ -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 # --------------------------------------------------------------- diff --git a/Makefile.standalone b/Makefile.standalone index 6359b46..5a58e4d 100644 --- a/Makefile.standalone +++ b/Makefile.standalone @@ -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 diff --git a/src/beamdata.f90 b/src/beamdata.f90 new file mode 100644 index 0000000..ad96107 --- /dev/null +++ b/src/beamdata.f90 @@ -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 diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index 26ce054..226d214 100755 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -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] diff --git a/src/gray-externals.f b/src/gray-externals.f index 43b8846..5d679f1 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -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) diff --git a/src/gray_main.f90 b/src/gray_main.f90 index 75a7752..1b31b4d 100644 --- a/src/gray_main.f90 +++ b/src/gray_main.f90 @@ -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 diff --git a/src/grayl.f b/src/grayl.f index 761950a..4ea42f5 100644 --- a/src/grayl.f +++ b/src/grayl.f @@ -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 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 diff --git a/src/itm_constants.f90 b/src/itm_constants.f90 deleted file mode 100644 index a617670..0000000 --- a/src/itm_constants.f90 +++ /dev/null @@ -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 diff --git a/src/itm_types.f90 b/src/itm_types.f90 deleted file mode 100644 index 8a16580..0000000 --- a/src/itm_types.f90 +++ /dev/null @@ -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 diff --git a/src/main.f90 b/src/main.f90 new file mode 100644 index 0000000..f5b71ef --- /dev/null +++ b/src/main.f90 @@ -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 \ No newline at end of file diff --git a/src/scatterspl.f90 b/src/scatterspl.f90 index 025ab65..b8ce1d0 100644 --- a/src/scatterspl.f90 +++ b/src/scatterspl.f90 @@ -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)