From 79af5f0b4b14827224974dd28f03e5120e7d342d Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Mon, 8 Jun 2015 14:06:03 +0000 Subject: [PATCH] added modules graydata_flags, graydata_par, graydata_anequil, interp_eqprof, magsurf_data in nocommon branch --- Makefile | 11 +- src/gray.f | 1132 +++++++++++++++++++++++++--------------------------- 2 files changed, 543 insertions(+), 600 deletions(-) diff --git a/Makefile b/Makefile index 0ffe52b..6bccb99 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,8 @@ EXE=gray # Objects list MAINOBJ=gray.o OTHOBJ=dispersion.o calcei_mod.o dqagmv.o grayl.o reflections.o green_func_p.o \ - const_and_precisions.o + const_and_precisions.o graydata_flags.o graydata_par.o graydata_anequil.o \ + magsurf_data.o interp_eqprof.o # Alternative search paths vpath %.f90 src @@ -23,10 +24,16 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules -gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions.o +gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions.o \ + graydata_flags.o graydata_par.o graydata_anequil.o magsurf_data.o interp_eqprof.o green_func_p.o: const_and_precisions.o reflections.o: const_and_precisions.o dispersion.o: calcei_mod.o dqagmv.o +graydata_flags.o: const_and_precisions.o +graydata_par.o: const_and_precisions.o +graydata_anequil.o: const_and_precisions.o +magsurf_data.o: const_and_precisions.o +interp_eqprof.o: const_and_precisions.o # General object compilation command %.o: %.f90 diff --git a/src/gray.f b/src/gray.f index 84853d9..2a2f908 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1,12 +1,13 @@ + use graydata_flags, only : ipass,igrad + implicit real*8 (a-h,o-z) + common/ierr/ierr - common/igrad/igrad common/mode/sox common/p0/p0mw common/powrfl/powrfl common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot - common/ipass/ipass c read data plus initialization @@ -67,15 +68,16 @@ c second pass into plasma subroutine gray_integration + use graydata_flags, only : dst + implicit none integer istep,istop,index_rt integer nstep - real*8 st,dst,strfl11 + real*8 st,strfl11 integer i real*8 st0 common/ss/st - common/dsds/dst common/istep/istep common/nstep/nstep common/istop/istop @@ -108,20 +110,16 @@ c ray integration: end subroutine after_gray_integration use const_and_precisions, only : zero + use graydata_flags, only : ibeam,iwarm,ilarm,iequil,iprof, + . filenmeqq,filenmprf,filenmbm + use graydata_anequil, only : dens0,te0 + implicit real*8 (a-h,o-z) - character*255 filenmeqq,filenmprf,filenmbm c common/ss/st - 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 common/taumnx/taumn,taumx,pabstot,currtot - common/pardens/dens0,aln1,aln2 - common/parqte/te0,dte0,alt1,alt2 c c print all ray positions in local reference system c @@ -168,9 +166,13 @@ c subroutine after_onestep(i,istop) use const_and_precisions, only : pi use reflections, only : inside + use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad + use graydata_par, only : psipol0,chipol0 + use interp_eqprof, only : zbmin,zbmax,rlim,zlim,nlim + implicit real*8 (a-h,o-z) complex*16 ext,eyt,extr,eytr,exin2,eyin2 - parameter(jmx=31,kmx=36,nmx=8000,nbb=5000) + parameter(jmx=31,kmx=36,nmx=8000) parameter(taucr=12.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) @@ -180,7 +182,6 @@ c dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3) dimension xwcl(3),xvjk(3,jmx,kmx),anvjk(3,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) - real*8 rlim(nbb),zlim(nbb) logical ins_pl,inside_plasma external inside_plasma c @@ -188,9 +189,7 @@ c 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 @@ -204,19 +203,13 @@ c common/xv/xv common/anv/anv c - common/pol0/psipol0,chipol0 common/polcof/psipol,chipol common/evt/ext,eyt common/yyrfl/yyrfl common/powrfl/powrfl common/strfl11/strfl11 - common/dsds/dst common/index_rt/index_rt - common/ipass/ipass - common/bound/zbmin,zbmax - common/limiter/rlim,zlim,nlim c - common/igrad/igrad common/wrk/ywrk,ypwrk c pabstot=0.0d0 @@ -385,6 +378,7 @@ c first pass in multi-pass simulation is stopped when at least one c ray has reflected and all rays are directed away from c reflection point, or when no reflection has occurred and c central ray re-enters the plasma + if((ipass.eq.1 .and. ((iopmin.gt.1) .or. . (taumn.lt.1d+30.and.taumn.gt.taucr))) . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. @@ -458,6 +452,9 @@ c c subroutine print_output(i,j,k) use const_and_precisions, only : pi + use graydata_flags, only : iequil,istpl0 + use graydata_par, only : psdbnd + implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0) @@ -485,8 +482,6 @@ c common/ss/st common/nray/nrayr,nrayth common/istgr/istpr,istpl - common/ist/istpr0,istpl0 - common/iieq/iequil c common/parpl/brr,bphi,bzz,ajphi common/btot/btot @@ -495,7 +490,6 @@ c common/dens/dens,ddens common/tete/tekev common/absor/alpha,effjcd,akim,tau0 - common/densbnd/psdbnd common/epolar/ex,ey,ez c common/nplr/anpl,anpr @@ -633,67 +627,33 @@ c use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_, * cvdr=>degree,pi use green_func_p, only:Setup_SpitzFunc + use graydata_flags + use graydata_par + use graydata_anequil + use interp_eqprof, only : rmxm,rlim,zlim,nlim,zbmin,zbmax, + . btrcen,rcen,alloc_lim + implicit real*8 (a-h,o-z) character*8 wdat character*10 wtim - character*255 filenmeqq,filenmprf,filenmbm - parameter(nmx=8000,nbb=5000) - real*8 rlim(nbb),zlim(nbb) + parameter(nmx=8000) c common/xgcn/xgcn - common/ipec/ipec,nnd common/nstep/nstep - common/ibeam/ibeam - common/ist/istpr0,istpl0 - common/warm/iwarm,ilarm - common/ieccd/ieccd - common/idst/idst - common/imx/imx -c - common/filesn/filenmeqq,filenmprf,filenmbm c common/nray/nrayr,nrayth - common/rwmax/rwmax - common/dsds/dst - common/igrad/igrad - common/ipass/ipass - common/rwallm/rwallm - common/limiter/rlim,zlim,nlim - common/iieq/iequil - common/icocos/icocos - common/ixp/ixp - common/ipsn/ipsinorm - common/sspl/sspl - common/iipr/iprof - common/factb/factb - common/facttn/factt,factn - common/cent/btrcen,rcen c common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 - common/pol0/psipol0,chipol0 - common/ipol/ipol -c - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/pardens/dens0,aln1,aln2 - common/parban/b0,rr0m,zr0m,rpam - common/parqq/q0,qa,alq - common/parqte/te0,dte0,alt1,alt2 - common/zz/Zeffan - common/bound/zbmin,zbmax c common/parbres/bres - common/densbnd/psdbnd - common/nfile/neqdsk,nprof - common/sgnib/sgnbphi,sgniphi common/p0/p0mw c common/mode/sox common/angles/alpha0,beta0 - common/scal/iscal c open(2,file='gray.data',status= 'unknown') c @@ -948,6 +908,8 @@ c versus psi, rhop, rhot if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 + call alloc_lim(ierr) + if (ierr.ne.0) stop rlim(1)=rwallm rlim(2)=max(r00*1.d-2,rmxm) rlim(3)=rlim(2) @@ -1020,22 +982,33 @@ c c subroutine surf_anal use const_and_precisions, only : pi + use graydata_anequil, only : b0,rr0m,zr0m,rpam + use magsurf_data, only : npsi,npoints,psicon,rcon,zcon, + . alloc_surf_anal + implicit real*8(a-h,o-z) - common/parban/b0,rr0m,zr0m,rpam common/parbres/bres + + npsi=10 + npoints=101 + call alloc_surf_anal(ierr) + if (ierr.ne.0) stop c c print circular magnetic surfaces iequil=1 c write(71,*) '#i psin R z' dal=2.0d-2*pi drn=0.1d0 - do i=1,10 + do i=1,npsi rni=i*drn - do k=1,101 + psicon(i)=rni + do k=1,npoints drrm=rpam*rni*cos((k-1)*dal) dzzm=rpam*rni*sin((k-1)*dal) rrm=rr0m+drrm zzm=zr0m+dzzm + rcon(i,k)=rrm + zcon(i,k)=zrm write(71,111) i,rni,rrm,zzm end do write(71,*) ' ' @@ -1080,8 +1053,9 @@ c c c subroutine read_beams + use graydata_flags, only : filenmbm, ibeam + implicit real*8(a-h,o-z) - character*255 filenmeqq,filenmprf,filenmbm parameter(nstrmx=50) c dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) @@ -1094,8 +1068,6 @@ c dimension phi1v(nstrmx),phi2v(nstrmx) dimension cphi1(nstrmx,4),cphi2(nstrmx,4) c - common/ibeam/ibeam - common/filesn/filenmeqq,filenmprf,filenmbm common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 common/angles/alpha0,beta0 @@ -1193,62 +1165,25 @@ c c subroutine equidata use const_and_precisions, only : pi + use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk + use graydata_par, only : sgnbphi,sgniphi,factb + use interp_eqprof, only : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz, + . psia,psiant,psinop,btrcen,rcen,btaxis,rmaxis,zmaxis,rmnm, + . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,fpolas,phitedge,rrtor,rup,zup, + . rlw,zlw,tr,tz,tfp,cc=>cceq,cfp,cc01=>cceq01,cc10=>cceq10, + . cc20=>cceq20,cc02=>cceq02,cc11=>cceq11,psinr,qpsi,rv,zv,psin, + . psi,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd + implicit real*8 (a-h,o-z) - parameter(nnw=501,nnh=501) - parameter(nbb=5000) c parameter(np=100) + parameter(kspl=3) character*48 stringa - dimension fpol(nnw),pres(nnw),qpsi(nnw) - dimension ffprim(nnw),pprim(nnw) - dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) - dimension rbbbs(nbb),zbbbs(nbb) - dimension rlim(nbb),zlim(nbb) c dimension rcon(2*np+1),zcon(2*np+1) + integer, dimension(:), allocatable :: iwrk,iwrkf,iwrkbsp + real*8, dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol, + . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim c - parameter(nrest=nnw+4,nzest=nnh+4) - parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54) - 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) - 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(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh) - dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) - parameter(lwrkf=nnw*4+nrest*16) - dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw) - dimension fpoli(nnw) -c - common/pareq1/psia - common/pareq1t/psiant,psinop - common/cent/btrcen,rcen - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/psinr/psinr - common/qpsi/qpsi - common/cpsin/rv,zv,psin - common/cpsi/psi - common/eqnn/nr,nz,npp,nintp - common/ipsn/ipsinorm - common/sspl/sspl - common/nfile/neqdsk,nprof - common/bound/zbmin,zbmax - common/sgnib/sgnbphi,sgniphi - common/factb/factb - common/ixp/ixp - common/icocos/icocos - - common/coffeqt/tr,tz - common/coffeqtp/tfp - common/coffeq/cc - common/coffeqd/cc01,cc10,cc20,cc02,cc11 - common/coffeqn/nsrt,nszt,nsft - common/cofffp/cfp - common/fpas/fpolas common/rhotsx/rhotsx - common/phitedge/phitedge - common/cnt/rup,zup,rlw,zlw - common/limiter/rlim,zlim,nlim c c read from file eqdsk c (see http://fusion.gat.com/efit/g_eqdsk.html) @@ -1256,10 +1191,46 @@ c c spline interpolation of psi and derivatives c if(icocos.gt.0) then - read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz + read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz else - read (neqdsk,*) nr,nz + read (neqdsk,*) nr,nz end if +c + call alloc_equilvec(ierr) + if (ierr.ne.0) stop + + nrest=nr+4 + nzest=nz+4 + lwrk=4*(nr+nz)+11*(nrest+nzest)+nrest*nz+nzest+54 + liwrk=nz+nr+nrest+nzest+3 + lw10=nr*3+nz*4+nr*nz + lw01=nr*4+nz*3+nr*nz + lw20=nr*2+nz*4+nr*nz + lw02=nr*4+nz*2+nr*nz + lw11=nr*3+nz*3+nr*nz + lwrkf=nr*4+nrest*16 + lwrkbsp=4*(nr+nz) + liwrkbsp=nr+nz + + if(allocated(fvpsi)) deallocate(fvpsi) + if(allocated(pres)) deallocate(pres) + if(allocated(ffprim)) deallocate(ffprim) + if(allocated(pprim)) deallocate(pprim) + if(allocated(ffvpsi)) deallocate(ffvpsi) + if(allocated(fpol)) deallocate(fpol) + if(allocated(fpoli)) deallocate(fpoli) + if(allocated(wrkf)) deallocate(wrkf) + if(allocated(iwrkf)) deallocate(iwrkf) + if(allocated(wf)) deallocate(wf) + if(allocated(wrk)) deallocate(wrk) + if(allocated(iwrk)) deallocate(iwrk) + if(allocated(wrkbsp)) deallocate(wrkbsp) + if(allocated(iwrkbsp)) deallocate(iwrkbsp) + + allocate(fvpsi(nr*nz),pres(nr),ffprim(nr),pprim(nr), + . ffvpsi(nr*nz),fpol(nr),fpoli(nr),wrkf(lwrkf),iwrkf(nrest), + . wf(nr),wrk(lwrk),iwrk(liwrk),wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)) +c if(ipsinorm.eq.0) then read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid read (neqdsk,2020) rmaxis,zmaxis,psiaxis,psiedge,btrcen @@ -1458,6 +1429,10 @@ c c read plasma boundaries from eqdsk file c read (neqdsk,*) nbbbs,nlim + + call alloc_bnd(ierr) + if (ierr.ne.0) stop + if(nbbbs.gt.0) then if(ipsinorm.eq.1) . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) @@ -1598,7 +1573,7 @@ c compute flux surface averaged quantities zlw=zmaxis-(zmaxis-zbmin)/10.0d0 call flux_average c ipr=1 -c call contours_psi(1.0d0,np,rcon,zcon,ipr) +c call contours_psi(1.0d0,rcon,zcon,ipr) c do ii=1,2*np+1 c write(52,*) rcon(ii), zcon(ii) c end do @@ -1630,6 +1605,8 @@ c locate btot=bres c call bfield_res c + deallocate(iwrk,iwrkf,iwrkbsp,fpoli,fvpsi,ffvpsi,fpol, + . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim) return end c @@ -1659,11 +1636,12 @@ c c c subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) + use interp_eqprof, only: psia + implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz - common/pareq1/psia select case(iflag) case(1) @@ -1706,13 +1684,14 @@ c c c subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) + use interp_eqprof, only: psia + implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/psival/psinv common/cnpsi/h - common/pareq1/psia select case(iflag) case(1) @@ -1737,13 +1716,11 @@ c c c subroutine print_prof + use interp_eqprof, only : psinr,qpsi,nr + implicit real*8 (a-h,o-z) - parameter(nnw=501,eps=1.d-4) - dimension psinr(nnw),qpsi(nnw) + parameter(eps=1.d-4) c - common/psinr/psinr - common/qpsi/qpsi - common/eqnn/nr,nz,npp,nintp common/dens/dens,ddens c write(55,*) ' #psi rhot ne Te q Jphi' @@ -1803,40 +1780,34 @@ c c c subroutine surfq(qval,psival) + use magsurf_data, only : npoints + use interp_eqprof, only : nr,psinr,qpsi + implicit real*8 (a-h,o-z) - parameter(nnw=501) - parameter(ncnt=100,ncntt=2*ncnt+1) - dimension psinr(nnw),qpsi(nnw) - dimension rcon(ncntt),zcon(ncntt) -c - common/psinr/psinr - common/qpsi/qpsi - common/eqnn/nr,nz,npp,nintp + dimension rcn(npoints),zcn(npoints) c c locate psi surface for q=qval c + ncnt=(npoints-1)/2 call locate(qpsi,nr,qval,i1) call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1), . qval,psival) ipr=1 - call contours_psi(psival,ncnt,rcon,zcon,ipr) + call contours_psi(psival,rcn,zcn,ipr) return end c c c subroutine bfield_res + use interp_eqprof, only : rv,zv,nr,nz,btotal + implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) - dimension rv(nnw),zv(nnh),psin(nnw,nnh) - dimension btotal(nnw,nnh) parameter(icmx=2002) dimension rrcb(icmx),zzcb(icmx),ncpts(10) c - common/cpsin/rv,zv,psin - common/eqnn/nr,nz,npp,nintp common/parbres/bres - common/btt/btotal c c Btotal on psi grid c @@ -1875,28 +1846,15 @@ c c c subroutine profdata + use graydata_flags, only : iprof,iscal,nprof + use graydata_par, only : psdbnd,factb,factt,factn + use interp_eqprof, only : nsfd,npp,psnpp,aa,bb,cc, + . psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec + implicit real*8 (a-h,o-z) - parameter(npmx=250,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) - dimension tfn(npest),cfn(npest),wrkf(lwrkf),iwrkf(npest),wf(npmx) - dimension densi(npest),ddensi(npest),d2densi(npest),wrkfd(npest) -c - common/densbnd/psdbnd - common/denspp/psnpp,aa,bb,cc - common/eqnn/nr,nz,npp,nintp - common/iipr/iprof - common/nfile/neqdsk,nprof - common/crad/psrad,derad,terad,zfc - common/coffte/ct - common/coffz/cz - common/factb/factb - common/coffdt/tfn - common/coffdnst/nsfd - common/cofffn/cfn - common/scal/iscal - common/facttn/factt,factn + integer, dimension(:), allocatable :: iwrkf + real*8, dimension(:), allocatable :: wrkf,wf,wrkfd, + . densi,ddensi,d2densi c c read plasma profiles from file if iprof>0 c @@ -1912,6 +1870,24 @@ c c if (iprof.gt.0) then read(nprof,*) npp + + call alloc_profvec(ierr) + if (ierr.ne.0) stop + + npest=npp+4 + lwrkf=npp*4+npest*16 + + if(allocated(wrkf)) deallocate(wrkf) + if(allocated(iwrkf)) deallocate(iwrkf) + if(allocated(wf)) deallocate(wf) + if(allocated(wrkfd)) deallocate(wrkfd) + if(allocated(densi)) deallocate(densi) + if(allocated(ddensi)) deallocate(ddensi) + if(allocated(d2densi)) deallocate(d2densi) + + allocate(wrkf(lwrkf),iwrkf(npest),wf(npp),densi(npest), + . wrkfd(npest),ddensi(npest),d2densi(npest)) + read(nprof,*) psrad0,terad0,derad0,zfc0 if(psrad0.ne.0.0d0) psrad0=0.0d0 psrad(1)=psrad0 @@ -1931,31 +1907,31 @@ c c spline approximation of temperature and Zeff c iopt=0 - call difcsn(psrad,terad,npmx,npp,iopt,ct,ier) + call difcsn(psrad,terad,npp,npp,iopt,ct,ier) c iopt=0 - call difcsn(psrad,zfc,npmx,npp,iopt,cz,ier) + call difcsn(psrad,zfc,npp,npp,iopt,cz,ier) c c spline approximation of density c - iopt=0 - xb=0.0d0 - xe=psrad(npp) - kspl=3 - sspl=.001d0 + iopt=0 + xb=0.0d0 + xe=psrad(npp) + 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) c - call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd, - . tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier) + call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier) + nu=1 + call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier) + dnpp=densi(npp) + ddnpp=ddensi(npp) c - call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier) - nu=1 - call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier) - dnpp=densi(npp) - ddnpp=ddensi(npp) -c - nu=2 - call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier) - d2dnpp=d2densi(npp) + nu=2 + call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier) + d2dnpp=d2densi(npp) if(derad(npp).eq.0.0d0) then psdbnd=psrad(npp) @@ -1975,52 +1951,56 @@ c c end if c + deallocate(iwrkf,wrkf,wf,densi,ddensi,d2densi,wrkfd) + return end c c c - subroutine rhotor(nr) + subroutine rhotor(nnr) + use interp_eqprof, only : nr,psinr,qpsi,crhot,cq + implicit real*8(a-h,o-z) - parameter(nnw=501) - dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4) - common/psinr/psinr - common/qpsi/qpsi + real*8, dimension(:), allocatable :: rhotnr common/rhotsx/rhotsx - common/crhot/crhot - common/cq/cq +c + if(allocated(cq)) deallocate(cq) + if(allocated(crhot)) deallocate(crhot) + allocate(rhotnr(nr),cq(nr,4),crhot(nr,4)) c c normalized toroidal rho : ~ Integral q(psi) dpsi c iopt=0 - call difcsn(psinr,qpsi,nnw,nr,iopt,cq,ier) + call difcsn(psinr,qpsi,nr,nnr,iopt,cq,ier) c rhotnr(1)=0.0d0 - do k=1,nr-1 + do k=1,nnr-1 dx=psinr(k+1)-psinr(k) drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0 . +dx*cq(k,4)/4.0d0))) rhotnr(k+1)=rhotnr(k)+drhot end do - rhotsx=rhotnr(nr) + rhotsx=rhotnr(nnr) do k=1,nr - rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr)) + rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nnr)) end do c c spline interpolation of rhotor c iopt=0 - call difcs(psinr,rhotnr,nr,iopt,crhot,ier) + call difcsn(psinr,rhotnr,nr,nnr,iopt,crhot,ier) + + deallocate(rhotnr) + return end function fq_eq(psi) + use interp_eqprof, only : psinr,nr,cq + implicit real*8(a-h,o-z) - parameter(nnw=501) - dimension psinr(nnw),cq(nnw,4) - common/psinr/psinr - common/eqnn/nr,nz,npp,nintp - common/cq/cq + irt=int((nr-1)*psi+1) if(irt.eq.0) irt=1 if(irt.eq.nr) irt=nr-1 @@ -2030,12 +2010,9 @@ c end function frhotor_eq(psi) + use interp_eqprof, only : psinr,nr,crhot + implicit real*8(a-h,o-z) - parameter(nnw=501) - dimension psinr(nnw),crhot(nnw,4) - common/psinr/psinr - common/eqnn/nr,nz,npp,nintp - common/crhot/crhot c irt=int((nr-1)*psi+1) c if(irt.eq.0) irt=1 @@ -2047,28 +2024,26 @@ c if(irt.eq.nr) irt=nr-1 end function frhotor(psi) + use graydata_flags, only : iequil + implicit real*8(a-h,o-z) - common/iieq/iequil if(iequil.eq.2) frhotor=frhotor_eq(psi) if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) return end function frhotor_av(psi) + use magsurf_data, only : rpstab, crhotq, npsi + implicit real*8(a-h,o-z) - parameter(nnintp=101) - dimension rpstab(nnintp),crhotq(nnintp,4) - common/pstab/rpstab - common/eqnn/nr,nz,npp,nintp - common/crhotq/crhotq rpsi=sqrt(psi) - ip=int((nintp-1)*rpsi+1) + ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 -c if(ip.eq.nintp) ip=nintp-1 - ip=min(max(1,ip),nintp-1) +c if(ip.eq.npsi) ip=npsi-1 + ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) - frhotor_av=spli(crhotq,nintp,ip,dps) + frhotor_av=spli(crhotq,npsi,ip,dps) return end @@ -2124,23 +2099,24 @@ c end do end subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi) + use interp_eqprof, only : dr,dz,nr,nz,rqgrid=>rv,zqgrid=>zv,psin, + . btotal + implicit real*8 (a-h,o-z) c c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c - parameter(nnw=501,nnh=501,icmx=2002,nna=nnw*nnh) - dimension a(nna),ja(3,2),lx(1000),npts(10) + parameter(icmx=2002) + dimension ja(3,2),lx(1000),npts(10) dimension rcon(icmx),zcon(icmx) - dimension rqgrid(nnw),zqgrid(nnh),psin(nnw,nnh),btotal(nnw,nnh) -c - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/cpsin/rqgrid,zqgrid,psin - common/btt/btotal - common/eqnn/nr,nz,npp,nintp + real*8, dimension(:), allocatable :: a c data px/0.5d0/ c + if(allocated(a)) deallocate(a) + allocate(a(nr*nz)) +c if(ichoi.eq.0) then do j=1,nz do i=1,nr @@ -2308,28 +2284,30 @@ c npts(ncon) = icount - iclast iclast = icount 50 continue +c + deallocate(a) c return end c c c - subroutine contours_psi(h,np,rcn,zcn,ipr) + subroutine contours_psi(h,rcn,zcn,ipr) + use graydata_par, only : rwallm + use magsurf_data, only : npoints + use interp_eqprof, only : psia,psiant,psinop,nsr=>nsrt,nsz=>nszt, + . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) parameter(mest=4,kspl=3) - parameter(nnw=501,nnh=501) - parameter(nrest=nnw+4,nzest=nnh+4) - dimension rcn(2*np+1),zcn(2*np+1) - dimension cc(nnw*nnh),tr(nrest),tz(nzest) - dimension czc(nrest),zeroc(mest) + dimension zeroc(mest) + dimension rcn(npoints),zcn(npoints) + real*8, dimension(:), allocatable :: czc c - common/pareq1t/psiant,psinop - common/coffeqn/nsr,nsz,nsft - common/coffeq/cc - common/coffeqt/tr,tz - common/cnt/rup,zup,rlw,zlw - common/rwallm/rwallm + np=(npoints-1)/2 + if(allocated(czc)) deallocate(czc) + allocate(czc(nr+4)) c ra=rup rb=rlw @@ -2341,8 +2319,8 @@ c th=pi/dble(np) rcn(1)=rlw zcn(1)=zlw - rcn(2*np+1)=rlw - zcn(2*np+1)=zlw + rcn(npoints)=rlw + zcn(npoints)=zlw rcn(np+1)=rup zcn(np+1)=zup do ic=2,np @@ -2355,22 +2333,25 @@ c if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) zcn(ic)=zc - rcn(2*np+2-ic)=zeroc(2) - zcn(2*np+2-ic)=zc + rcn(npoints+1-ic)=zeroc(2) + zcn(npoints+1-ic)=zc else rcn(ic)=zeroc(2) zcn(ic)=zc - rcn(2*np+2-ic)=zeroc(3) - zcn(2*np+2-ic)=zc + rcn(npoints+1-ic)=zeroc(3) + zcn(npoints+1-ic)=zc end if end do if (ipr.gt.0) then - do ii=1,2*np+1 + do ii=1,npoints write(71,111) ii,h,rcn(ii),zcn(ii) end do write(71,*) write(71,*) end if + + deallocate(czc) + return 111 format(i6,12(1x,e12.5)) end @@ -2379,61 +2360,48 @@ c c subroutine flux_average use const_and_precisions, only : zero,one,pi,ccj=>mu0inv + use magsurf_data + use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge + implicit real*8 (a-h,o-z) real*8 lam + real*8, dimension(:), allocatable :: rctemp,zctemp - parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) + parameter(nnintp=101,ncnt=100,nlam=41) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54) parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) - - dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) - dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) - dimension dadrhotv(nnintp),cdadrhot(nnintp,4) - dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4) - dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) - dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) - dimension rcon(2*ncnt+1),zcon(2*ncnt+1) +c + dimension dadrhotv(nnintp) + dimension dvdrhotv(nnintp) dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) - dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4) - dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4) - dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4) - dimension cfc(nnintp,4),crhotq(nnintp,4) - dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp) + dimension vratjpl(nnintp) dimension alam(nlam),fhlam(nnintp,nlam) dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) - dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)) dimension iwrk(kwrk),wrk(lwrk) - dimension ch01(lw01),weights(nlam) + dimension weights(nlam) c - common/cent/btrcen,rcen - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/eqnn/nr,nz,npp,nintp common/fpol/fpolv,ffpv c - common/pstab/rpstab - common/flav/vvol,rri,rbav,bmxpsi,bmnpsi - common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc - common/cratj/cratja,cratjb,cratjpl - common/crhotq/crhotq - common/phitedge/phitedge - common/cdadrhot/cdadrhot - common/cdvdrhot/cdvdrhot + npsi=nnintp + ninpr=(npsi-1)/10 + npoints = 2*ncnt+1 c - common/coffvl/ch - common/coffdvl/ch01 - common/coffvlt/tjp,tlm - common/coffvln/njpt,nlmt + call alloc_surfvec(ierr) + if(allocated(tjp)) deallocate(tjp) + if(allocated(tlm)) deallocate(tlm) + if(allocated(ch)) deallocate(ch) + if(allocated(ch01)) deallocate(ch01) + allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), + . ch01(lw01),rctemp(npoints),zctemp(npoints),stat=ierr) + if (ierr.ne.0) return c computation of flux surface averaged quantities write(71,*)' #i psin R z' - nintp=nnintp - ninpr=(nintp-1)/10 - dlam=1.0d0/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam @@ -2460,7 +2428,11 @@ c computation of flux surface averaged quantities qq=1.0d0 fc=1.0d0 + psicon(1)=0.0d0 + rcon(1,:)=0.0d0 + zcon(1,:)=0.0d0 pstab(1)=0.0d0 + rhot_eq(1)=0.0d0 rpstab(1)=0.0d0 rhotqv(1)=0.0d0 vcurrp(1)=0.0d0 @@ -2479,19 +2451,25 @@ c computation of flux surface averaged quantities rhot2q=0.0d0 qqv(1)=1.0d0 + dadrhotv(1)=0.0d0 + dvdrhotv(1)=0.0d0 + C rup=rmaxis C rlw=rmaxis C zup=(zbmax+zmaxis)/2.0d0 C zlw=(zmaxis+zbmin)/2.0d0 - do jp=2,nintp - height=dble(jp-1)/dble(nintp-1) - if(jp.eq.nintp) height=0.9999d0 + do jp=2,npsi + height=dble(jp-1)/dble(npsi-1) + psicon(jp)=height + if(jp.eq.npsi) height=0.9999d0 ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 height2=height*height - call contours_psi(height2,ncnt,rcon,zcon,ipr) + call contours_psi(height2,rctemp,zctemp,ipr) + rcon(jp,:) = rctemp + zcon(jp,:) = zctemp c r2iav=0.0d0 anorm=0.0d0 @@ -2505,37 +2483,37 @@ c bmmx=-1.0d+30 bmmn=1.0d+30 - call tor_curr(rcon(1),zcon(1),ajphi0) - call bpol(rcon(1),zcon(1),brr,bzz) + call tor_curr(rctemp(1),zctemp(1),ajphi0) + call bpol(rctemp(1),zctemp(1),brr,bzz) call equinum_fpol(0) - bphi=fpolv/rcon(1) + bphi=fpolv/rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 - rpsim0=rcon(1) + rpsim0=rctemp(1) - do inc=1,ncntt-1 + do inc=1,npoints-1 inc1=inc+1 - dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2) - dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2) - dlp=sqrt((rcon(inc1)-rcon(inc))**2+ - . (zcon(inc1)-zcon(inc))**2) - drc=(rcon(inc1)-rcon(inc)) + dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) + dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) + dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ + . (zctemp(inc1)-zctemp(inc))**2) + drc=(rctemp(inc1)-rctemp(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5d0*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) - rzp=rcon(inc1)*zcon(inc1) - rz=rcon(inc)*zcon(inc) + rzp=rctemp(inc1)*zctemp(inc1) + rz=rctemp(inc)*zctemp(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 - rpsim=rcon(inc1) - zpsim=zcon(inc1) + rpsim=rctemp(inc1) + zpsim=zctemp(inc1) call bpol(rpsim,zpsim,brr,bzz) call tor_curr(rpsim,zpsim,ajphi) call equinum_fpol(0) @@ -2612,13 +2590,11 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = . *dadpsi/pi dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) . *dvdpsi/pi -c -c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) - . /dble(nintp-1) + . /dble(npsi-1) rhotqv(jp)=sqrt(rhot2q) c computation of fraction of circulating/trapped fraction fc, ft @@ -2633,7 +2609,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ rl2=1.0d0-lam*bv(1)/bmmx rl0=0.d0 if(rl2.gt.0) rl0=sqrt(rl2) - do inc=1,ncntt-1 + do inc=1,npoints-1 rl2=1.0d0-lam*bv(inc+1)/bmmx rl=0.0d0 if(rl2.gt.0) rl=sqrt(rl2) @@ -2670,58 +2646,58 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) - do jp=1,nintp - rhotqv(jp)=rhotqv(jp)/rhotqv(nintp) - if(jp.eq.nintp) then - rhotqv(jp)=1.0d0 - rpstab(jp)=1.0d0 - pstab(jp)=1.0d0 + do jp=1,npsi + rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) + if(jp.eq.npsi) then + rhotqv(jp)=1.0d0 + rpstab(jp)=1.0d0 + pstab(jp)=1.0d0 end if - rhot_eq=frhotor_eq(pstab(jp)) - write(56,99) pstab(jp),rhot_eq,rhotqv(jp), + rhot_eq(jp)=frhotor_eq(pstab(jp)) + write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), - . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) - . ,vratja(jp),vratjb(jp) + . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), + . vratja(jp),vratjb(jp) end do -c rarea=sqrt(varea(nintp)/pi) +c rarea=sqrt(varea(npsi)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 - call difcs(rpstab,vvol,nintp,iopt,cvol,ier) + call difcs(rpstab,vvol,npsi,iopt,cvol,ier) iopt=0 - call difcs(rpstab,rbav,nintp,iopt,crbav,ier) + call difcs(rpstab,rbav,npsi,iopt,crbav,ier) iopt=0 - call difcs(rpstab,rri,nintp,iopt,crri,ier) + call difcs(rpstab,rri,npsi,iopt,crri,ier) iopt=0 - call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier) + call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) iopt=0 - call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier) + call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) iopt=0 - call difcs(rpstab,vratja,nintp,iopt,cratja,ier) + call difcs(rpstab,vratja,npsi,iopt,cratja,ier) iopt=0 - call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier) + call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) iopt=0 - call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier) + call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) iopt=0 - call difcs(rpstab,varea,nintp,iopt,carea,ier) + call difcs(rpstab,varea,npsi,iopt,carea,ier) iopt=0 - call difcs(rpstab,ffc,nintp,iopt,cfc,ier) + call difcs(rpstab,ffc,npsi,iopt,cfc,ier) iopt=0 - call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier) + call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) iopt=0 - call difcs(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier) + call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) iopt=0 - call difcs(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier) + call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0d0 - call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam, + call regrid(iopt,npsi,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 @@ -2736,40 +2712,41 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda end function fdadrhot(rpsi) + use magsurf_data, only : rpstab,cdadrhot,npsi + implicit real*8(a-h,o-z) - parameter(nnintp=101) - dimension rpstab(nnintp),cdadrhot(nnintp,4) - common/pstab/rpstab - common/eqnn/nr,nz,npp,nintp - common/cdadrhot/cdadrhot - ip=int((nintp-1)*rpsi+1) + + ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 -c if(ip.eq.nintp) ip=nintp-1 - ip=min(max(1,ip),nintp-1) +c if(ip.eq.npsi) ip=npsi-1 + ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) - fdadrhot=spli(cdadrhot,nintp,ip,dps) + fdadrhot=spli(cdadrhot,npsi,ip,dps) return end function fdvdrhot(rpsi) + use magsurf_data, only : rpstab,cdvdrhot,npsi + implicit real*8(a-h,o-z) - parameter(nnintp=101) - dimension rpstab(nnintp),cdvdrhot(nnintp,4) - common/pstab/rpstab - common/eqnn/nr,nz,npp,nintp - common/cdvdrhot/cdvdrhot - ip=int((nintp-1)*rpsi+1) - ip=min(max(1,ip),nintp-1) + + ip=int((npsi-1)*rpsi+1) + ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) - fdvdrhot=spli(cdvdrhot,nintp,ip,dps) + fdvdrhot=spli(cdvdrhot,npsi,ip,dps) return end subroutine flux_average_an + use graydata_anequil, only : q0,qa,alq,b0,rr0m,zr0m,rpam + use magsurf_data + use interp_eqprof, only : btrcen + implicit real*8 (a-h,o-z) real*8 lam + real*8, dimension(:), allocatable :: rctemp,zctemp - parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) + parameter(nnintp=101,ncnt=100,nlam=41) parameter(zero=0.0d0,one=1.0d0) parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) @@ -2778,49 +2755,34 @@ c if(ip.eq.nintp) ip=nintp-1 parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) - dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) - dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) - dimension dadrhotv(nnintp),cdadrhot(nnintp,4) - dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4) - dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) - dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) - dimension rcon(2*ncnt+1),zcon(2*ncnt+1) + dimension dadrhotv(nnintp) + dimension dvdrhotv(nnintp) dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) - dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4) - dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4) - dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4) - dimension cfc(nnintp,4),crhotq(nnintp,4) - dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp) + dimension vratjpl(nnintp) dimension alam(nlam),fhlam(nnintp,nlam) dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) - dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)) dimension iwrk(kwrk),wrk(lwrk) - dimension ch01(lw01),weights(nlam) + dimension weights(nlam) c - common/cent/btrcen,rcen - common/parban/b0,rr0m,zr0m,rpam - common/parqq/q0,qa,alq - common/eqnn/nr,nz,npp,nintp common/fpol/fpolv,ffpv c common/psival/psin common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz c - common/pstab/rpstab - common/flav/vvol,rri,rbav,bmxpsi,bmnpsi - common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc - common/cratj/cratja,cratjb,cratjpl - common/crhotq/crhotq - common/phitedge/phitedge - common/cdadrhot/cdadrhot - common/cdvdrhot/cdvdrhot + npsi=nnintp + ninpr=(npsi-1)/10 + npoints = 2*ncnt+1 +c + call alloc_surfvec(ierr) + if(allocated(tjp)) deallocate(tjp) + if(allocated(tlm)) deallocate(tlm) + if(allocated(ch)) deallocate(ch) + if(allocated(ch01)) deallocate(ch01) + allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), + . ch01(lw01),rctemp(npoints),zctemp(npoints),stat=ierr) + if (ierr.ne.0) return c - common/coffvl/ch - common/coffdvl/ch01 - common/coffvlt/tjp,tlm - common/coffvln/njpt,nlmt - c computation of flux surface averaged quantities rmaxis=rr0m @@ -2832,9 +2794,6 @@ c computation of flux surface averaged quantities write(71,*)' #i psin R z' - nintp=nnintp - ninpr=(nintp-1)/10 - dlam=1.0d0/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam @@ -2861,7 +2820,11 @@ c computation of flux surface averaged quantities qq=1.0d0 fc=1.0d0 + psicon(1)=0.0d0 + rcon(1,:)=0.0d0 + zcon(1,:)=0.0d0 pstab(1)=0.0d0 + rhot_eq(1)=0.0d0 rpstab(1)=0.0d0 rhotqv(1)=0.0d0 vcurrp(1)=0.0d0 @@ -2887,14 +2850,17 @@ C rlw=rmaxis C zup=(zbmax+zmaxis)/2.0d0 C zlw=(zmaxis+zbmin)/2.0d0 - do jp=2,nintp - height=dble(jp-1)/dble(nintp-1) - if(jp.eq.nintp) height=0.9999d0 + do jp=2,npsi + height=dble(jp-1)/dble(npsi-1) + psicon(jp)=height + if(jp.eq.npsi) height=0.9999d0 height2=height*height ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 - call contours_psi_an(height2,ncnt,rcon,zcon,ipr) + call contours_psi_an(height2,rctemp,zctemp,ipr) + rcon(jp,:) = rctemp + zcon(jp,:) = zctemp c r2iav=0.0d0 anorm=0.0d0 @@ -2908,40 +2874,40 @@ c bmmx=-1.0d+30 bmmn=1.0d+30 - call equian(rcon(1),zcon(1)) - dbvcdc13=-ddpsidzz/rcon(1) - dbvcdc31= ddpsidrr/rcon(1)-bzz/rcon(1) + call equian(rctemp(1),zctemp(1)) + dbvcdc13=-ddpsidzz/rctemp(1) + dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) ajphi=ccj*(dbvcdc13-dbvcdc31) - brr=-dpsidz/rcon(1) - bzz= dpsidr/rcon(1) - bphi=fpolv/rcon(1) + brr=-dpsidz/rctemp(1) + bzz= dpsidr/rctemp(1) + bphi=fpolv/rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 - rpsim0=rcon(1) + rpsim0=rctemp(1) - do inc=1,ncntt-1 + do inc=1,npoints-1 inc1=inc+1 - dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2) - dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2) - dlp=sqrt((rcon(inc1)-rcon(inc))**2+ - . (zcon(inc1)-zcon(inc))**2) - drc=(rcon(inc1)-rcon(inc)) + dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) + dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) + dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ + . (zctemp(inc1)-zctemp(inc))**2) + drc=(rctemp(inc1)-rctemp(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5d0*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) - rzp=rcon(inc1)*zcon(inc1) - rz=rcon(inc)*zcon(inc) + rzp=rctemp(inc1)*zctemp(inc1) + rz=rctemp(inc)*zctemp(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 - rpsim=rcon(inc1) - zpsim=zcon(inc1) + rpsim=rctemp(inc1) + zpsim=zctemp(inc1) call equian(rpsim,zpsim) brr=-dpsidz/rpsim bzz= dpsidr/rpsim @@ -3029,7 +2995,7 @@ c dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) - . /dble(nintp-1) + . /dble(npsi-1) c print*,jp,rhot2q,qqv(jp),rpstab(jp),qqv(jp-1),rpstab(jp-1) rhotqv(jp)=sqrt(rhot2q) c rhotqv(jp)=rn @@ -3049,7 +3015,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ rl2=1.0d0-lam*bv(1)/bmmx rl0=0.d0 if(rl2.gt.0) rl0=sqrt(rl2) - do inc=1,ncntt-1 + do inc=1,npoints-1 rl2=1.0d0-lam*bv(inc+1)/bmmx rl=0.0d0 if(rl2.gt.0) rl=sqrt(rl2) @@ -3086,58 +3052,58 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) - do jp=1,nintp - rhotqv(jp)=rhotqv(jp)/rhotqv(nintp) - if(jp.eq.nintp) then + do jp=1,npsi + rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) + if(jp.eq.npsi) then rhotqv(jp)=1.0d0 rpstab(jp)=1.0d0 pstab(jp)=1.0d0 end if - rhot_eq=frhotor_an(sqrt(pstab(jp))) - write(56,99) pstab(jp),rhot_eq,rhotqv(jp), + rhot_eq(jp)=frhotor_an(sqrt(pstab(jp))) + write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), - . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) - . ,vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) + . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), + . vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) end do -c rarea=sqrt(varea(nintp)/pi) +c rarea=sqrt(varea(npsi)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 - call difcs(rpstab,vvol,nintp,iopt,cvol,ier) + call difcs(rpstab,vvol,npsi,iopt,cvol,ier) iopt=0 - call difcs(rpstab,rbav,nintp,iopt,crbav,ier) + call difcs(rpstab,rbav,npsi,iopt,crbav,ier) iopt=0 - call difcs(rpstab,rri,nintp,iopt,crri,ier) + call difcs(rpstab,rri,npsi,iopt,crri,ier) iopt=0 - call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier) + call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) iopt=0 - call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier) + call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) iopt=0 - call difcs(rpstab,vratja,nintp,iopt,cratja,ier) + call difcs(rpstab,vratja,npsi,iopt,cratja,ier) iopt=0 - call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier) + call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) iopt=0 - call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier) + call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) iopt=0 - call difcs(rpstab,varea,nintp,iopt,carea,ier) + call difcs(rpstab,varea,npsi,iopt,carea,ier) iopt=0 - call difcs(rpstab,ffc,nintp,iopt,cfc,ier) + call difcs(rpstab,ffc,npsi,iopt,cfc,ier) iopt=0 - call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier) + call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) iopt=0 - call difcs(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier) + call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) iopt=0 - call difcs(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier) + call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0d0 - call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam, + call regrid(iopt,npsi,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 @@ -3153,6 +3119,10 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda subroutine rhopol_an + use graydata_par, only : sgniphi + use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam + use interp_eqprof, only : psia + implicit real*8(a-h,o-z) parameter(nnr=101,nrest=nnr+4) parameter(lwrk=nnr*4+nrest*16) @@ -3167,10 +3137,6 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda common/coffrptt/trot common/coffrnt/nsrot common/coffrpt/crot - common/parqq/q0,qa,alq - common/parban/b0,rr0m,zr0m,rpam - common/pareq1/psia - common/sgnib/sgnbphi,sgniphi dr=1.0d0/dble(nnr-1) @@ -3238,16 +3204,20 @@ c spline interpolation of rhotor versus rhopol end - subroutine contours_psi_an(h,np,rcn,zcn,ipr) + subroutine contours_psi_an(h,rcn,zcn,ipr) + use graydata_anequil, only : rr0m,zr0m,rpam + use magsurf_data, only : npoints + implicit real*8 (a-h,o-z) parameter(pi=3.14159265358979d0) parameter(mest=4,kspl=3) - dimension rcn(2*np+1),zcn(2*np+1) - common/parban/b0,rr0m,zr0m,rpam + dimension rcn(npoints),zcn(npoints) + + np=(npoints-1)/2 th=pi/dble(np) rn=frhotor_an(sqrt(h)) - do ic=1,2*np+1 + do ic=1,npoints zcn(ic)=zr0m+rpam*rn*sin(th*(ic-1)) rcn(ic)=rr0m+rpam*rn*cos(th*(ic-1)) @@ -3266,6 +3236,8 @@ c spline interpolation of rhotor versus rhopol subroutine vectinit use dispersion, only: expinit + use graydata_flags, only : iwarm + implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) dimension psjki(jmx,kmx,nmx) @@ -3275,7 +3247,6 @@ c spline interpolation of rhotor versus rhopol dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) dimension istore(jmx,kmx),anwcl(3),xwcl(3) c - common/warm/iwarm,ilarm common/iiv/iiv common/iov/iop,iow,ihcd,istore common/psjki/psjki @@ -3635,6 +3606,8 @@ c c Runge-Kutta integrator c subroutine rkint4 + use graydata_flags, only : dst,igrad + implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36) dimension y(ndim),yy(ndim) @@ -3645,7 +3618,6 @@ c 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 @@ -3654,7 +3626,6 @@ c common/ggradjk/ggri common/gr/gr2 common/dgr/dgr2,dgr,ddgr - common/igrad/igrad c h=dst hh=h*0.5d0 @@ -3706,6 +3677,8 @@ c c c subroutine gwork(j,k) + use graydata_flags, only : igrad + implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36) dimension yy(ndim),yyp(ndim) @@ -3713,8 +3686,6 @@ c 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 @@ -3755,6 +3726,8 @@ c c c subroutine fwork(y,dery) + use graydata_flags, only : idst,igrad + implicit real*8 (a-h,o-z) parameter(ndim=6) dimension y(ndim),dery(ndim) @@ -3769,7 +3742,6 @@ c common/nplr/anpl,anpr common/bb/bv common/dbb/derbv - common/igrad/igrad common/xgxg/xg common/ygyg/yg common/dxgyg/derxg,deryg @@ -3778,7 +3750,6 @@ c common/ierr/ierr common/anv/anv common/xv/xv - common/idst/idst c xx=y(1) yy=y(2) @@ -3938,6 +3909,9 @@ c c subroutine plas_deriv(xx,yy,zz) use const_and_precisions, only : pi + use graydata_flags, only : iequil + use graydata_par, only : sgnbphi + 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) @@ -3951,12 +3925,10 @@ c common/dxgdps/dxgdpsi common/ygyg/yg common/dxgyg/derxg,deryg - common/iieq/iequil common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/fpol/fpolv,ffpv common/psival/psinv - common/sgnib/sgnbphi,sgniphi xg=0.0d0 yg=9.9d1 @@ -4105,12 +4077,13 @@ c c c subroutine equian(rrm,zzm) + use graydata_par, only : sgnbphi,sgniphi + use graydata_anequil, only : b0,rr0m,zr0m,rpam,q0,qa,alq + use interp_eqprof, only : psia + implicit real*8 (a-h,o-z) c - common/parqq/q0,qa,alq - common/parban/b0,rr0m,zr0m,rpam common/psival/psinv - common/pareq1/psia common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/fpol/fpolv,ffpv @@ -4118,7 +4091,6 @@ c common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddens - common/sgnib/sgnbphi,sgniphi dpsidrp=0.0d0 d2psidrp=0.0d0 @@ -4171,22 +4143,16 @@ c subroutine equinum_psi(rpsim,zpsim) + use interp_eqprof, only : rmnm,rmxm,zmnm,zmxm,psiant,psinop, + . tr,tz,ccspl=>cceq,nsrt,nszt + implicit real*8 (a-h,o-z) - parameter(nnw=501,nnh=501) - parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=8,liwrk=2) - dimension ccspl(nnw*nnh) - dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) + dimension wrk(lwrk),iwrk(liwrk) dimension rrs(1),zzs(1),ffspl(1) parameter(nrs=1,nzs=1) c common/psival/psinv - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/pareq1t/psiant,psinop -c - common/coffeqt/tr,tz - common/coffeq/ccspl - common/coffeqn/nsrt,nszt,nsft c psinv=-1.0d0 c @@ -4207,18 +4173,20 @@ c end subroutine equinum_derpsi(rpsim,zpsim,iderpsi) + use interp_eqprof, only : nr,nz,rmnm,rmxm,zmnm,zmxm,cc01=>cceq01, + . cc10=>cceq10,cc20=>cceq20,cc02=>cceq02,cc11=>cceq11 + implicit real*8 (a-h,o-z) - parameter(nnw=501,nnh=501) - 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(lw11=nnw*3+nnh*3+nnw*nnh) - dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) integer*4 iderpsi c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/coffeqd/cc01,cc10,cc20,cc02,cc11 +c + lw10=nr*3+nz*4+nr*nz + lw01=nr*4+nz*3+nr*nz + lw20=nr*2+nz*4+nr*nz + lw02=nr*4+nz*2+nr*nz + lw11=nr*3+nz*3+nr*nz c dpsidr=0.0d0 dpsidz=0.0d0 @@ -4287,20 +4255,15 @@ c here lengths are measured in meters end subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw) + use interp_eqprof, only : nr,nz,psia,tr,tz,nsrt,nszt + implicit real*8 (a-h,o-z) - parameter(nnw=501,nnh=501) - parameter(nrest=nnw+4,nzest=nnh+4) parameter(liwrk=2) parameter(nrs=1,nzs=1) - dimension tr(nrest),tz(nzest),iwrk(liwrk) + dimension iwrk(liwrk) dimension rrs(1),zzs(1),ffspl(1) dimension cc(lw) - common/pareq1/psia - common/coffeqt/tr,tz - common/coffeqn/nsrt,nszt,nsft - common/eqnn/nr,nz,npp,nintp - rrs(1)=rpsim zzs(1)=zpsim nsr=nsrt @@ -4317,20 +4280,15 @@ c here lengths are measured in meters end subroutine equinum_fpol(iderfpol) + use interp_eqprof, only : nr,psia,tfp,nsft,cfp,fpolas + implicit real*8 (a-h,o-z) - parameter(nnw=501) - parameter(nrest=nnw+4) dimension rrs(1),ffspl(1) - dimension tfp(nrest),cfp(nrest),wrkfd(nrest) + dimension wrkfd(nr+4) integer*4 iderfpol c common/psival/psinv common/fpol/fpolv,ffpv - common/pareq1/psia - common/coffeqtp/tfp - common/coffeqn/nsrt,nszt,nsft - common/cofffp/cfp - common/fpas/fpolas fpolv=fpolas dfpolv=0.0d0 @@ -4378,8 +4336,10 @@ c end subroutine tor_curr_psi(h,ajphi) + use interp_eqprof, only : zmaxis + implicit real*8 (a-h,o-z) - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz + call psi_raxis(h,r1,r2) call tor_curr(r2,zmaxis,ajphi) return @@ -4399,18 +4359,16 @@ c end c subroutine psi_raxis(h,r1,r2) + use interp_eqprof, only : psiant,psinop,zmaxis,nsr=>nsrt, + . nsz=>nszt,cc=>cceq,tr,tz,nr + implicit real*8 (a-h,o-z) parameter(mest=4,kspl=3) - parameter(nnw=501,nnh=501) - parameter(nrest=nnw+4,nzest=nnh+4) - dimension cc(nnw*nnh),tr(nrest),tz(nzest) - dimension czc(nrest),zeroc(mest) + dimension zeroc(mest) + real*8, dimension(:), allocatable :: czc c - common/pareq1t/psiant,psinop - common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz - common/coffeqn/nsr,nsz,nsft - common/coffeq/cc - common/coffeqt/tr,tz + if(allocated(czc)) deallocate(czc) + allocate(czc(nr+4)) c iopt=1 zc=zmaxis @@ -4420,15 +4378,19 @@ c call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) + + deallocate(czc) + return end c c c subroutine sub_xg_derxg + use interp_eqprof, only : psia + implicit real*8 (a-h,o-z) common/psival/psinv - common/pareq1/psia common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn @@ -4446,19 +4408,17 @@ c c c subroutine density(arg) + use graydata_flags, only : iprof + use graydata_par, only : psdbnd + use graydata_anequil, only : dens0,aln1,aln2 + use interp_eqprof, only : psnpp,aad=>aa,bbd=>bb,ccd=>cc,tfn,nsfd, + . cfn,npp + implicit real*8 (a-h,o-z) - parameter(npmx=250,npest=npmx+4) dimension xxs(1),ffs(1) - dimension tfn(npest),cfn(npest),wrkfd(npest) + dimension wrkfd(npp+4) c - common/densbnd/psdbnd - common/denspp/psnpp,aad,bbd,ccd - common/iipr/iprof - common/pardens/dens0,aln1,aln2 common/dens/dens,ddens - common/coffdt/tfn - common/coffdnst/nsfd - common/cofffn/cfn c c computation of density [10^19 m^-3] and derivative wrt psi c @@ -4506,15 +4466,11 @@ c c c function temp(arg) + use graydata_flags, only : iprof + use graydata_anequil, only : te0,dte0,alt1,alt2 + use interp_eqprof, only : psrad,ct,npp + implicit real*8 (a-h,o-z) - parameter(npmx=250) - dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4) -c - common/parqte/te0,dte0,alt1,alt2 - common/iipr/iprof - common/crad/psrad,derad,terad,zfc - common/coffte/ct - common/eqnn/nr,nz,npp,nintp c temp=0.0d0 if(arg.ge.1.0d0.or.arg.lt.0.0d0) return @@ -4526,7 +4482,7 @@ c if(k.eq.0) k=1 if(k.eq.npp) k=npp-1 dps=arg-psrad(k) - temp=spli(ct,npmx,k,dps) + temp=spli(ct,npp,k,dps) endif return end @@ -4534,15 +4490,11 @@ c c c function fzeff(arg) + use graydata_flags, only : iprof + use graydata_anequil, only : zeffan + use interp_eqprof, only : psrad,cz,npp + implicit real*8 (a-h,o-z) - parameter(npmx=250) - dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4) -c - common/iipr/iprof - common/crad/psrad,derad,terad,zfc - common/coffz/cz - common/eqnn/nr,nz,npp,nintp - common/zz/Zeffan c fzeff=1 ps=arg @@ -4553,7 +4505,7 @@ c call locate(psrad,npp,ps,k) k=max(1,min(k,npp-1)) dps=ps-psrad(k) - fzeff=spli(cz,npmx,k,dps) + fzeff=spli(cz,npp,k,dps) endif return end @@ -4562,6 +4514,9 @@ c beam tracing initial conditions igrad=1 c subroutine ic_gb use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree + use graydata_flags, only : ipol + use graydata_par, only : rwmax,psipol0,chipol0 + implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36) @@ -4579,7 +4534,6 @@ c external catand 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 @@ -4593,8 +4547,6 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/evt/ext,eyt - common/pol0/psipol0,chipol0 - common/ipol/ipol common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi @@ -4878,6 +4830,9 @@ c subroutine ic_rt use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree, * ui=>im + use graydata_flags, only : ipol + use graydata_par, only : rwmax,psipol0,chipol0 + implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36) @@ -4890,7 +4845,6 @@ c complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) 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 @@ -4903,8 +4857,6 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/evt/ext,eyt - common/pol0/psipol0,chipol0 - common/ipol/ipol common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi @@ -5073,6 +5025,8 @@ c subroutine ic_rt2 use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree + use graydata_par, only : psipol0,chipol0 + implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36) @@ -5094,7 +5048,6 @@ c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/yyrfl/yyrfl common/evt/ext,eyt - common/pol0/psipol0,chipol0 common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi @@ -5189,12 +5142,13 @@ c c ray power weigth coefficient q(j) c subroutine pweigth + use graydata_par, only : rwmax + 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 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) @@ -5228,34 +5182,29 @@ c c subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, . bmxi,bmni,fci,intp) + use magsurf_data, only : rpstab,cvol,crri,crbav,cbmx,cbmn, + . carea,cfc,npsi + implicit real*8 (a-h,o-z) - parameter(nnintp=101) - dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4) - dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4) - dimension carea(nnintp,4),cfc(nnintp,4) c - common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc - common/pstab/rpstab - common/eqnn/nr,nz,npp,nintp -c - ip=int((nintp-1)*rpsi+1) + ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 -c if(ip.eq.nintp) ip=nintp-1 - ip=min(max(1,ip),nintp-1) +c if(ip.eq.npsi) ip=npsi-1 + ip=min(max(1,ip),npsi-1) c dps=rpsi-rpstab(ip) c - areai=spli(carea,nintp,ip,dps) - voli=spli(cvol,nintp,ip,dps) - dervoli=splid(cvol,nintp,ip,dps) - rrii=spli(crri,nintp,ip,dps) + areai=spli(carea,npsi,ip,dps) + voli=spli(cvol,npsi,ip,dps) + dervoli=splid(cvol,npsi,ip,dps) + rrii=spli(crri,npsi,ip,dps) c if(intp.eq.0) return c - rbavi=spli(crbav,nintp,ip,dps) - bmxi=spli(cbmx,nintp,ip,dps) - bmni=spli(cbmn,nintp,ip,dps) - fci=spli(cfc,nintp,ip,dps) + rbavi=spli(crbav,npsi,ip,dps) + bmxi=spli(cbmx,npsi,ip,dps) + bmni=spli(cbmn,npsi,ip,dps) + fci=spli(cfc,npsi,ip,dps) c return end @@ -5263,21 +5212,18 @@ c c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) + use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi + implicit real*8 (a-h,o-z) - parameter(nnintp=101) - dimension rpstab(nnintp) - dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4) - common/pstab/rpstab - common/eqnn/nr,nz,npp,nintp - common/cratj/cratja,cratjb,cratjpl - ip=int((nintp-1)*rpsi+1) + + ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 -c if(ip.eq.nintp) ip=nintp-1 - ip=min(max(1,ip),nintp-1) +c if(ip.eq.npsi) ip=npsi-1 + ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) - ratjai=spli(cratja,nintp,ip,dps) - ratjbi=spli(cratjb,nintp,ip,dps) - ratjpli=spli(cratjpl,nintp,ip,dps) + ratjai=spli(cratja,npsi,ip,dps) + ratjbi=spli(cratjb,npsi,ip,dps) + ratjpli=spli(cratjpl,npsi,ip,dps) return end c @@ -5285,6 +5231,10 @@ c c subroutine pabs_curr(i,j,k) use const_and_precisions, only : pi + use graydata_flags, only : iequil,dst + use graydata_par, only : sgnbphi + use graydata_anequil, only : b0,rr0m,rpam + implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) dimension psjki(jmx,kmx,nmx) @@ -5303,15 +5253,11 @@ c common/qw/q common/p0/p0mw c - common/dsds/dst common/dersdst/dersdst - common/iieq/iequil c - common/parban/b0,rr0m,zr0m,rpam common/absor/alpha,effjcd,akim,tau0 c common/psival/psinv - common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci c @@ -5359,14 +5305,13 @@ c subroutine ecrh_cd use const_and_precisions, only : mc2=>mc2_ use dispersion, only : harmnumber, warmdisp + use graydata_flags, only : iwarm,ilarm,ieccd,imx + implicit real*8 (a-h,o-z) parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) complex*16 ex,ey,ez c common/nharm/nhmin,nhmax - common/warm/iwarm,ilarm - common/imx/imx - common/ieccd/ieccd c common/xgxg/xg common/ygyg/yg @@ -5437,6 +5382,7 @@ c use green_func_p, only: SpitzFuncCoeff use const_and_precisions, only : pi,qesi=>e_,mesi=>me_,vcsi=>c_, * qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ + implicit none real*8 anum,denom,resp,resj,effjcd,coullog,dens real*8 yg,anpl,anpr,amu,anprre,anprim @@ -5950,25 +5896,22 @@ c end subroutine vlambda(alam,psi,fv,dfv) + use magsurf_data, only : ch,ch01,tjp,tlm,njpt,nlmt,npsi + implicit real*8 (a-h,o-z) - parameter (nnintp=101,nlam=41) - parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) - parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ - . njest*nnintp+nlest+54) - parameter(kwrk=nnintp+nlam+njest+nlest+3) - parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) + parameter (nlam=41) + parameter(ksp=3,nlest=nlam+ksp+1) external fpbisp dimension xxs(1),yys(1),ffs(1) - dimension ch01(lw01),ch((njest-4)*(nlest-4)) - dimension tjp(njest),tlm(nlest) - dimension iwrk(kwrk),wrk(lwrk) + integer, dimension(:), allocatable :: iwrk + real*8, dimension(:), allocatable :: wrk - common/coffvl/ch - common/coffdvl/ch01 - common/coffvlt/tjp,tlm - common/coffvln/njpt,nlmt + njest=npsi+ksp+1 + kwrk=npsi+nlam+njest+nlest+3 + lwrk=4*(npsi+nlam)+11*(njest+nlest)+njest*npsi+nlest+54 + allocate(iwrk(kwrk),wrk(lwrk)) njp=njpt nlm=nlmt @@ -6070,6 +6013,9 @@ c c subroutine pec(pabs,currt) use const_and_precisions, only : pi,one + use graydata_flags, only : ipec,ieccd,iequil,nnd,dst + use graydata_anequil, only : rr0m,rpam + implicit real*8(a-h,o-z) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) parameter(rtbc=one) @@ -6089,9 +6035,6 @@ c 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 @@ -6100,8 +6043,6 @@ c common/dpjjki/pdjki,currj c common/angles/alpha0,beta0 - common/iieq/iequil - common/parban/b0,rr0m,zr0m,rpam common/taumnx/taumn,taumx,pabstot,currtot c common/polcof/psipol,chipol @@ -6417,11 +6358,11 @@ c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) use const_and_precisions, only : emn1 + use graydata_flags, only : ipec,iequil + implicit real*8(a-h,o-z) dimension xx(nd),yy(nd) common/jmxmn/rhotp,rhotm,ypkp,ypkm - common/ipec/ipec,nnd - common/iieq/iequil c call vmaxmini(yy,nd,ymn,ymx,imn,imx) ypk=0.0d0 @@ -6656,14 +6597,14 @@ c ell=bb/aa logical function inside_plasma(rrm,zzm) - implicit none - real*8 rrm,zzm,psdbnd,psinv,zbmin,zbmax - integer iequil + use graydata_flags, only : iequil + use graydata_par, only : psdbnd + use interp_eqprof, only : zbmin,zbmax + + implicit none + real*8 rrm,zzm,psinv - common/densbnd/psdbnd - common/bound/zbmin,zbmax common/psival/psinv - common/iieq/iequil if(iequil.eq.1) then call equian(rrm,zzm) @@ -6688,6 +6629,7 @@ c ell=bb/aa . irfl) use reflections, only : inters_linewall,inside use const_and_precisions, only : ui=>im,pi + use interp_eqprof, only : rlim,zlim,nlim implicit none integer*4 irfl real*8 anv(3),anv0(3),xv(3),xvrfl(3) @@ -6695,11 +6637,6 @@ c ell=bb/aa real*8 smax,rrm,zzm complex*16 extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) - integer nbb,nlim - parameter(nbb=5000) - real*8 rlim(nbb),zlim(nbb) - - common/limiter/rlim,zlim,nlim anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) @@ -6767,19 +6704,18 @@ c wave vector and electric field after reflection in lab frame subroutine vacuum_rt(xvstart,anv,xvend,ivac) use reflections, only : inters_linewall,inside + use graydata_flags, only : dst + use interp_eqprof, only : rlim,zlim,nlim + implicit none integer*4 ivac - integer nbb,nlim,i - parameter(nbb=5000) - real*8 st,rrm,zzm,dst,dstvac,smax + integer i + real*8 st,rrm,zzm,dstvac,smax real*8 anv(3),xvstart(3),xvend(3),walln(3) real*8 xv0(3),anv0(3) - real*8 rlim(nbb),zlim(nbb) logical plfound,inside_plasma external inside_plasma - common/limiter/rlim,zlim,nlim - common/dsds/dst common/dstvac/dstvac c ivac=1 plasma hit before wall reflection @@ -6817,7 +6753,7 @@ c ivac=-1 vessel (and thus plasma) never crossed xvend=xv0+st*anv0 rrm=1.d-2*sqrt(xvend(1)**2+xvend(2)**2) zzm=1.d-2*xvend(3) - plfound=inside_plasma(rrm,zzm) + plfound=inside_plasma(rrm,zzm,0) if (st.ge.smax.or.plfound) exit i=i+1 end do