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)