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