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

made beam arrays allocatable and max # of knots in spline for scattered data reduced to n/4 to save memory

Browse Source
This commit is contained in:
Lorenzo Figini 2014-06-13 15:12:36 +00:00
parent b6d8dd5a6f
commit 9ca1ccd817
4 changed files with 106 additions and 234 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Makefile
View File

@ -97,11 +97,12 @@ clean:
# Dependencies # Dependencies
# ------------ # ------------
gray_main.o: gray-externals.o itm_types.o gray_main.o: gray-externals.o itm_types.o
gray-externals.o: green_func_p.o reflections.o scatterspl.o gray-externals.o: green_func_p.o reflections.o scatterspl.o beamdata.o
green_func_p.o: const_and_precisions.o green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o const_and_precisions.o: itm_types.o itm_constants.o
itm_constants.o: itm_types.o itm_constants.o: itm_types.o
scatterspl.o: itm_types.o scatterspl.o: itm_types.o
beamdata.o: itm_types.o
# Special rule to preprocess source file and include svn revision # Special rule to preprocess source file and include svn revision
# --------------------------------------------------------------- # ---------------------------------------------------------------

5
Makefile.standalone
View File

@ -4,7 +4,7 @@ LIBFILE=lib$(EXE).a
# Objects list # Objects list
OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \ OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.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 itm_constants.o itm_types.o
# Alternative search paths # Alternative search paths
vpath %.f90 src vpath %.f90 src
@ -30,11 +30,12 @@ $(LIBFILE): $(OBJ)
# Dependencies on modules # Dependencies on modules
main.o: gray_main.o main.o: gray_main.o
gray_main.o: gray-externals.o itm_types.o gray_main.o: gray-externals.o itm_types.o
gray-externals.o: green_func_p.o reflections.o scatterspl.o gray-externals.o: green_func_p.o reflections.o scatterspl.o beamdata.o
green_func_p.o: const_and_precisions.o green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o const_and_precisions.o: itm_types.o itm_constants.o
itm_constants.o: itm_types.o itm_constants.o: itm_types.o
scatterspl.o: itm_types.o scatterspl.o: itm_types.o
beamdata.o: itm_types.o
# General object compilation command # General object compilation command
%.o: %.f90 %.o: %.f90

325
src/gray-externals.f
View File

@ -1,6 +1,7 @@
subroutine gray_integration subroutine gray_integration
use beamdata, only : nstep
implicit none implicit none
integer istep,nstep,istop,index_rt integer istep,istop,index_rt
real*8 st,dst,strfl11 real*8 st,dst,strfl11
integer i integer i
real*8 st0 real*8 st0
@ -8,7 +9,6 @@
common/ss/st common/ss/st
common/dsds/dst common/dsds/dst
common/istep/istep common/istep/istep
common/nstep/nstep
common/istop/istop common/istop/istop
common/strfl11/strfl11 common/strfl11/strfl11
common/index_rt/index_rt common/index_rt/index_rt
@ -38,6 +38,7 @@ c ray integration: end
subroutine after_gray_integration(rhopin,nrho,dpdvout,ajcdout) subroutine after_gray_integration(rhopin,nrho,dpdvout,ajcdout)
use beamdata, only : nrayr
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(zero=0.0d0) parameter(zero=0.0d0)
character*255 filenmeqq,filenmprf,filenmbm character*255 filenmeqq,filenmprf,filenmbm
@ -48,7 +49,6 @@ c
common/ibeam/ibeam common/ibeam/ibeam
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
common/nray/nrayr,nrayth
common/iieq/iequil common/iieq/iequil
common/iipr/iprof common/iipr/iprof
common/index_rt/index_rt common/index_rt/index_rt
@ -104,27 +104,16 @@ c
c c
subroutine after_onestep(i,istop) subroutine after_onestep(i,istop)
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : psjki,ppabs,ccci,iiv,tauv,
. iop,iow,tau1v,yyrfl,nrayr,nrayth
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
parameter(taucr=12.0d0,cvdr=pi/180.0d0) parameter(taucr=12.0d0,cvdr=pi/180.0d0)
dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension iop(jmx,kmx),iow(jmx,kmx)
dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
dimension xv(3),anv(3),xvrfl(3),anvrfl(3) dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
common/pcjki/ppabs,ccci
common/atjki/tauv,alphav
common/tau1v/tau1v
c
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
common/nray/nrayr,nrayth
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/istgr/istpr,istpl common/istgr/istpr,istpl
c c
common/iiv/iiv
common/iov/iop,iow
common/psjki/psjki
common/psival/psinv common/psival/psinv
common/psinv11/psinv11 common/psinv11/psinv11
common/ierr/ierr common/ierr/ierr
@ -138,7 +127,6 @@ c
common/ipol/ipolc common/ipol/ipolc
common/iovmin/iopmin,iowmin common/iovmin/iopmin,iowmin
common/densbnd/psdbnd common/densbnd/psdbnd
common/yyrfl/yyrfl
common/powrfl/powrfl common/powrfl/powrfl
common/dstvac/dstvac common/dstvac/dstvac
common/strfl11/strfl11 common/strfl11/strfl11
@ -297,32 +285,19 @@ c
c c
subroutine print_output(i,j,k) subroutine print_output(i,j,k)
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : ywrk,psjki,tauv,alphav,pdjki,
. currj,didst,q,tau1v,nrayr!,nrayth
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36,nmx=8000)
parameter(taucr=12.0d0) parameter(taucr=12.0d0)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension q(jmx),tau1v(jmx,kmx)
complex*16 ex,ey,ez complex*16 ex,ey,ez
c c
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nhn/nhn common/nhn/nhn
common/iokh/iohkawa common/iokh/iohkawa
common/p0/p0mw common/p0/p0mw
common/tau1v/tau1v
common/qw/q
common/index_rt/index_rt common/index_rt/index_rt
common/ss/st common/ss/st
common/nray/nrayr,nrayth
common/istgr/istpr,istpl common/istgr/istpr,istpl
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/iieq/iequil common/iieq/iequil
@ -340,8 +315,6 @@ c
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
common/nprw/anprr,anpri common/nprw/anprr,anpri
c
common/wrk/ywrk,ypwrk
c c
x=ywrk(1,j,k) x=ywrk(1,j,k)
y=ywrk(2,j,k) y=ywrk(2,j,k)
@ -436,7 +409,7 @@ c
c c
subroutine prfile subroutine prfile
implicit none implicit none
integer*4 index_rt integer index_rt
common/index_rt/index_rt common/index_rt/index_rt
If(index_rt.eq.1) then If(index_rt.eq.1) then
@ -475,6 +448,7 @@ c
. dne, zeff, qsf, powin) . dne, zeff, qsf, powin)
use green_func_p, only:Setup_SpitzFunc use green_func_p, only:Setup_SpitzFunc
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : nrayr,nrayth,nstep
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
integer ijetto, mr, mz, nrho, nbnd integer ijetto, mr, mz, nrho, nbnd
real*8 r(mr), z(mz), psin(mr,mz) real*8 r(mr), z(mz), psin(mr,mz)
@ -489,13 +463,12 @@ c
character*255 filenmeqq,filenmprf,filenmbm character*255 filenmeqq,filenmprf,filenmbm
parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
parameter(cvdr=pi/180.0d0) parameter(cvdr=pi/180.0d0)
parameter(nmx=8000,nbb=5000) parameter(nbb=5000)
real*8 rlim(nbb),zlim(nbb) real*8 rlim(nbb),zlim(nbb)
c c
common/xgcn/xgcn common/xgcn/xgcn
common/ipec/ipec,nnd common/ipec/ipec,nnd
common/nstep/nstep
common/ibeam/ibeam common/ibeam/ibeam
common/ist/istpr0,istpl0 common/ist/istpr0,istpl0
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
@ -504,7 +477,6 @@ c
c c
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
c c
common/nray/nrayr,nrayth
common/rwmax/rwmax common/rwmax/rwmax
common/dsds/dst common/dsds/dst
common/igrad/igrad common/igrad/igrad
@ -802,9 +774,9 @@ c set simple limiter as two cylindrical walls at rwallm and r00
rlim(3)=rlim(2) rlim(3)=rlim(2)
rlim(4)=rlim(1) rlim(4)=rlim(1)
rlim(5)=rlim(1) rlim(5)=rlim(1)
zlim(1)=-dst*nmx*1.d-2 zlim(1)=-dst*nstep*2.d-2
zlim(2)=zlim(1) zlim(2)=zlim(1)
zlim(3)=dst*nmx*1.d-2 zlim(3)=dst*nstep*2.d-2
zlim(4)=zlim(3) zlim(4)=zlim(3)
zlim(5)=zlim(1) zlim(5)=zlim(1)
ipass=abs(ipass) ipass=abs(ipass)
@ -1041,8 +1013,8 @@ c
parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3) parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3)
dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh) dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh)
dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk)
parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh) c parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh)
dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) c dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
dimension rv1d(nnw*nnh),zv1d(nnw*nnh),wpsi(nnw*nnh) dimension rv1d(nnw*nnh),zv1d(nnw*nnh),wpsi(nnw*nnh)
parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh)
parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh)
@ -1174,29 +1146,32 @@ c valid data only inside boundary
end do end do
end do end do
c add boundary to the fit c add boundary to the fit
if (ij+nbnd.le.nnw*nnh) then c
do i=1,nbnd c *** TEMPORARILY COMMENTED TO TRY REDUCING MEMORY CONSUMPTION - BEGIN ***
ij=ij+1 c if (ij+nbnd.le.nnw*nnh) then
rv1d(ij)=rbnd(i) c do i=1,nbnd
zv1d(ij)=zbnd(i) c ij=ij+1
fvpsi(ij)=1.0d0 c rv1d(ij)=rbnd(i)
wpsi(ij)=1.0d0 c zv1d(ij)=zbnd(i)
end do c fvpsi(ij)=1.0d0
end if c wpsi(ij)=1.0d0
c end do
c end if
c *** TEMPORARILY COMMENTED TO TRY REDUCING MEMORY CONSUMPTION - END ***
c
nrz=ij nrz=ij
c c
c fit as a scattered set of points c fit as a scattered set of points
c c
nsr=nr+4 nsr=nr/4+4
nsz=nz+4 nsz=nz/4+4
call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
. rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier)
!!!!!!!!!!!! Sistemare anche dipendenze Makefile scatterspl:...
c if ier=-1 data are fitted using sspl=0 c if ier=-1 data are fitted using sspl=0
if(ier.eq.-1) then if(ier.eq.-1) then
sspl=0.0d0 sspl=0.0d0
nsr=nr+4 nsr=nr/4+4
nsz=nz+4 nsz=nz/4+4
call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
. rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier)
end if end if
@ -1530,7 +1505,7 @@ c
rhot=0.0d0 rhot=0.0d0
call density(psin) call density(psin)
call tor_curr_psi(eps,ajphi) call tor_curr_psi(eps,ajphi)
te=temp(psin) te=temperature(psin)
qq=qpsi(1) qq=qpsi(1)
c c
write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6 write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6
@ -1541,7 +1516,7 @@ c
rhop=sqrt(psin) rhop=sqrt(psin)
c c
call density(psin) call density(psin)
te=temp(psin) te=temperature(psin)
c c
call vlocate(psinr,nrho,psin,ips) call vlocate(psinr,nrho,psin,ips)
ips=min(max(ips,1),nrho-1) ips=min(max(ips,1),nrho-1)
@ -2518,33 +2493,23 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda
end end
subroutine vectinit subroutine vectinit
use beamdata, only : jmx,kmx,nmx,psjki,tauv,alphav,pdjki,ppabs,
. currj,didst,ccci,iiv,iop,iow,tau1v,
. nrayr,nrayth,nstep,alloc_beam
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx)
parameter(tmax=5,npts=500) parameter(tmax=5,npts=500)
real*8 ttv(npts+1),extv(npts+1) real*8 ttv(npts+1),extv(npts+1)
common/ttex/ttv,extv common/ttex/ttv,extv
c c
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
common/iiv/iiv c
common/iov/iop,iow common/ierr/ierr
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nray/nrayr,nrayth
common/nstep/nstep
common/tau1v/tau1v
c c
if(nstep.gt.nmx) nstep=nmx if(nstep.gt.nmx) nstep=nmx
if(nrayr.gt.jmx) nrayr=jmx if(nrayr.gt.jmx) nrayr=jmx
if(nrayth.gt.kmx) nrayth=kmx if(nrayth.gt.kmx) nrayth=kmx
call alloc_beam(ierr)
if (ierr.ne.0) return
c c
do i=1,nstep do i=1,nstep
do k=1,nrayth do k=1,nrayth
@ -2576,26 +2541,20 @@ c
return return
end end
subroutine vectfree
use beamdata, only : dealloc_beam
implicit none
c
call dealloc_beam
c
return
end
subroutine vectinit2 subroutine vectinit2
use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,
. didst,ccci,iiv,iop,iow,nrayr,nrayth,nstep
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx)
common/iiv/iiv
common/iov/iop,iow
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/nray/nrayr,nrayth
common/nstep/nstep
do i=1,nstep do i=1,nstep
do k=1,nrayth do k=1,nrayth
@ -2641,16 +2600,9 @@ c
c c
c c
subroutine updatepos subroutine updatepos
implicit real*8 (a-h,o-z) use beamdata, only : xc,xco,du1,du1o,ywrk,nrayr,nrayth
parameter(jmx=31,kmx=36)
dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c c
common/nray/nrayr,nrayth implicit real*8 (a-h,o-z)
common/grco/xco,du1o
common/grc/xc,du1
common/wrk/ywrk,ypwrk
c c
do j=1,nrayr do j=1,nrayr
do k=1,nrayth do k=1,nrayth
@ -2681,26 +2633,12 @@ c
c c
c c
subroutine gradi subroutine gradi
use beamdata, only : dffiu,ddffiu,xc,xco,du1,du1o,gri,ggri,
. dgrad2v,grad2,nrayr,nrayth
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36)
dimension dffiu(jmx),ddffiu(jmx)
dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
dimension dgrad2v(3,jmx,kmx)
dimension grad2(jmx,kmx)
dimension dxv1(3),dxv2(3),dxv3(3),dgu(3) dimension dxv1(3),dxv2(3),dxv3(3),dgu(3)
dimension dgg1(3),dgg2(3),dgg3(3) dimension dgg1(3),dgg2(3),dgg3(3)
dimension df1(3),df2(3),df3(3) dimension df1(3),df2(3),df3(3)
c
common/nray/nrayr,nrayth
common/fi/dffiu,ddffiu
common/grco/xco,du1o
common/grc/xc,du1
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
c c
c compute grad u and grad(S_I) c compute grad u and grad(S_I)
c c
@ -2890,23 +2828,16 @@ c
c Runge-Kutta integrator c Runge-Kutta integrator
c c
subroutine rkint4 subroutine rkint4
use beamdata, only : nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v,
. gri,ggri
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36) parameter(ndim=6)
dimension y(ndim),yy(ndim) dimension y(ndim),yy(ndim)
dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim) dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim)
dimension dgr2(3),dgr(3),ddgr(3,3) dimension dgr2(3),dgr(3),ddgr(3,3)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nray/nrayr,nrayth
common/dsds/dst common/dsds/dst
c c
common/wrk/ywrk,ypwrk
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/gr/gr2 common/gr/gr2
common/dgr/dgr2,dgr,ddgr common/dgr/dgr2,dgr,ddgr
common/igrad/igrad common/igrad/igrad
@ -2961,21 +2892,14 @@ c
c c
c c
subroutine gwork(j,k) subroutine gwork(j,k)
use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36) parameter(ndim=6)
dimension yy(ndim),yyp(ndim) dimension yy(ndim),yyp(ndim)
dimension dgr2(3),dgr(3),ddgr(3,3) dimension dgr2(3),dgr(3),ddgr(3,3)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/igrad/igrad common/igrad/igrad
c c
common/wrk/ywrk,ypwrk
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/gr/gr2 common/gr/gr2
common/dgr/dgr2,dgr,ddgr common/dgr/dgr2,dgr,ddgr
c c
@ -3011,8 +2935,7 @@ c
c c
subroutine fwork(y,dery) subroutine fwork(y,dery)
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6) dimension y(6),dery(6)
dimension y(ndim),dery(ndim)
dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3) dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
dimension derdxv(3),danpldxv(3),derdnv(3) dimension derdxv(3),danpldxv(3),derdnv(3)
dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3) dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
@ -3693,7 +3616,7 @@ c
c c
c c
c c
function temp(arg) function temperature(arg)
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(npmx=250) parameter(npmx=250)
dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4) dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4)
@ -3704,17 +3627,17 @@ c
common/coffte/ct common/coffte/ct
common/eqnn/nr,nz,nrho,npp,nintp common/eqnn/nr,nz,nrho,npp,nintp
c c
temp=0.0d0 temperature=0.0d0
if(arg.ge.1.0d0.or.arg.lt.0.0d0) return if(arg.ge.1.0d0.or.arg.lt.0.0d0) return
if(iprof.eq.0) then if(iprof.eq.0) then
proft=(1.0d0-arg**alt1)**alt2 proft=(1.0d0-arg**alt1)**alt2
temp=(te0-dte0)*proft+dte0 temperature=(te0-dte0)*proft+dte0
else else
call vlocate(psrad,npp,arg,k) call vlocate(psrad,npp,arg,k)
if(k.eq.0) k=1 if(k.eq.0) k=1
if(k.eq.npp) k=npp-1 if(k.eq.npp) k=npp-1
dps=arg-psrad(k) dps=arg-psrad(k)
temp=fspli(ct,npmx,k,dps) temperature=fspli(ct,npmx,k,dps)
endif endif
return return
end end
@ -3751,15 +3674,12 @@ c beam tracing initial conditions igrad=1
c c
subroutine ic_gb subroutine ic_gb
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,
. gri,ggri
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3) parameter(zero=0.0d0,izero=0,one=1.0d0)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(cvdr=pi/180.0d0) parameter(cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy
complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy
complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
@ -3767,19 +3687,11 @@ c
external catand external catand
c parameter(ui=(0.0d0,1.0d0)) c parameter(ui=(0.0d0,1.0d0))
c c
common/nray/nrayr,nrayth
common/rwmax/rwmax common/rwmax/rwmax
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c common/anic/anx0c,any0c,anz0c
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/fi/dffiu,ddffiu
common/mirr/x00,y00,z00 common/mirr/x00,y00,z00
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
c c
ui=(0.0d0,1.0d0) ui=(0.0d0,1.0d0)
@ -4032,28 +3944,17 @@ c ray tracing initial conditions igrad=0
c c
subroutine ic_rt subroutine ic_rt
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,
. gri,ggri
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3) parameter(zero=0.0d0,izero=0,one=1.0d0)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(cvdr=pi/180.0d0) parameter(cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nray/nrayr,nrayth
common/rwmax/rwmax common/rwmax/rwmax
common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
common/anic/anx0c,any0c,anz0c common/anic/anx0c,any0c,anz0c
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/fi/dffiu,ddffiu
common/mirr/x00,y00,z00 common/mirr/x00,y00,z00
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
c c
csth=anz0c csth=anz0c
@ -4188,26 +4089,14 @@ c
subroutine ic_rt2 subroutine ic_rt2
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
. xc0=>xc,du10=>du1,grad2,dgrad2v,
. gri,ggri,yyrfl
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3) parameter(zero=0.0d0,izero=0,one=1.0d0)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(cvdr=pi/180.0d0) parameter(cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension yyrfl(jmx,kmx,ndim)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
dimension dffiu(jmx),ddffiu(jmx)
dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c c
common/nray/nrayr,nrayth
common/wrk/ywrk0,ypwrk0
common/grc/xc0,du10
common/grad2jk/grad2
common/dgrad2jk/dgrad2v
common/gradjk/gri
common/ggradjk/ggri
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
common/yyrfl/yyrfl
do j=1,nrayr do j=1,nrayr
do k=1,nrayth do k=1,nrayth
@ -4287,11 +4176,8 @@ c
c ray power weigth coefficient q(j) c ray power weigth coefficient q(j)
c c
subroutine pweigth subroutine pweigth
use beamdata, only : nrayr,nrayth,q
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(jmx=31)
dimension q(jmx)
common/qw/q
common/nray/nrayr,nrayth
common/rwmax/rwmax common/rwmax/rwmax
c c
dr=1.0d0 dr=1.0d0
@ -4381,22 +4267,10 @@ c
c c
subroutine pabs_curr(i,j,k) subroutine pabs_curr(i,j,k)
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst,
. ccci,q,tau1v
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000)
dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension q(jmx),tau1v(jmx,kmx)
c c
common/psjki/psjki
common/atjki/tauv,alphav
common/dpjjki/pdjki,currj
common/pcjki/ppabs,ccci
common/dcjki/didst
common/tau1v/tau1v
c
common/qw/q
common/p0/p0mw common/p0/p0mw
c c
common/dsds/dst common/dsds/dst
@ -4491,7 +4365,7 @@ c
akim=0.0d0 akim=0.0d0
effjcd=0.0d0 effjcd=0.0d0
c c
tekev=temp(psinv) tekev=temperature(psinv)
amu=mc2/tekev amu=mc2/tekev
zeff=fzeff(psinv) zeff=fzeff(psinv)
c c
@ -5895,12 +5769,9 @@ c gg=F(u)/u with F(u) as in Cohen paper
subroutine projxyzt(iproj,nfile) subroutine projxyzt(iproj,nfile)
use beamdata, only : ywrk,nrayr,nrayth
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c c
common/nray/nrayr,nrayth
common/wrk/ywrk,ypwrk
common/psinv11/psinv11 common/psinv11/psinv11
common/istep/istep common/istep/istep
common/ss/st common/ss/st
@ -5975,35 +5846,23 @@ c
c c
subroutine pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) subroutine pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout)
use itm_constants, only : pi=>itm_pi use itm_constants, only : pi=>itm_pi
use beamdata, only : psjki,iiv,ppabs,ccci,pdjki,
. nrayr,nrayth,nstep
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000)
parameter(rtbc=1.0d0) parameter(rtbc=1.0d0)
dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho)
dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx) real*8, dimension(:), allocatable :: xxi,ypt,yamp
dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) real*8, dimension(:), allocatable :: rtab,rhotv,rtab1,ajphiv,
dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx) . dpdv,dvolt,darea,ratjav,ratjbv,ratjplv,ajplv,ajcdav,
dimension xxi(nmx),ypt(nmx),yamp(nmx) . ajcdbv,pins,currins,fi
dimension rtab(nndmx),rhotv(nndmx)
dimension rtab1(0:nndmx)
dimension ajphiv(nndmx),dpdv(nndmx)
dimension dvolt(nndmx),darea(nndmx)
dimension ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx)
dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx)
dimension pins(nndmx),currins(nndmx),fi(nndmx)
parameter(llmx=21) parameter(llmx=21)
dimension isev(llmx) dimension isev(llmx)
c c
common/nray/nrayr,nrayth
common/istep/istep common/istep/istep
common/dsds/dst common/dsds/dst
common/ipec/ipec,nnd common/ipec/ipec,nnd
common/ieccd/ieccd common/ieccd/ieccd
common/index_rt/index_rt common/index_rt/index_rt
c
common/iiv/iiv
common/psjki/psjki
common/pcjki/ppabs,ccci
common/dpjjki/pdjki,currj
c c
common/cent/btrcen,rcen common/cent/btrcen,rcen
common/angles/alpha0,beta0 common/angles/alpha0,beta0
@ -6014,6 +5873,7 @@ c
c c
common/polcof/psipol,chipol common/polcof/psipol,chipol
c c
allocate(xxi(nstep),ypt(nstep),yamp(nstep))
stf=istep*dst stf=istep*dst
if (ipec.ge.0) then if (ipec.ge.0) then
nd=nnd nd=nnd
@ -6021,6 +5881,10 @@ c
! ipec=-1 uses rho_pol input grid ! ipec=-1 uses rho_pol input grid
nd=nrho nd=nrho
end if end if
allocate(rtab(nd),rhotv(nd),rtab1(0:nd),ajphiv(nd),
. dpdv(nd),dvolt(nd),darea(nd),ratjav(nd),
. ratjbv(nd),ratjplv(nd),ajplv(nd),ajcdav(nd),
. ajcdbv(nd),pins(nd),currins(nd),fi(nd))
if(pabstot.gt.0.0d0) then if(pabstot.gt.0.0d0) then
do ll=1,llmx do ll=1,llmx
@ -6090,7 +5954,7 @@ c radial coordinate of i-(i+1) interval mid point
do k=1,kkk do k=1,kkk
ise0=0 ise0=0
ii=iiv(j,k) ii=iiv(j,k)
if (ii.lt.nmx.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 if (ii.lt.nstep.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1
idecr=-1 idecr=-1
is=1 is=1
do i=1,ii do i=1,ii
@ -6336,6 +6200,9 @@ c fill output arguments only if ipec=-1
end do end do
end if end if
deallocate(xxi,ypt,yamp)
deallocate(rtab,rhotv,rtab1,ajphiv,dpdv,dvolt,darea,ratjav,ratjbv,
. ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi)
return return
49 format(i5,20(1x,e12.5)) 49 format(i5,20(1x,e12.5))
99 format(30(1x,e12.5)) 99 format(30(1x,e12.5))

7
src/gray_main.f90
View File

@ -38,8 +38,10 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, & call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, &
nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin) nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin)
call vectinit call vectinit
if(igrad.eq.0) call ic_rt if(iercom.eq.0) then
if(igrad.gt.0) call ic_gb if(igrad.eq.0) call ic_rt
if(igrad.gt.0) call ic_gb
end if
if(iercom.gt.0) then if(iercom.gt.0) then
ierr=iercom ierr=iercom
write(*,*) ' IERR = ', ierr write(*,*) ' IERR = ', ierr
@ -93,6 +95,7 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
end if end if
write(*,*) write(*,*)
write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8 write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8
call vectfree
pec=pabstott*1.e6_r8 pec=pabstott*1.e6_r8
icd=currtott*1.e6_r8 icd=currtott*1.e6_r8
dpdv=dpdv*1.e6_r8 dpdv=dpdv*1.e6_r8