diff --git a/src/gray.f b/src/gray.f index 58d7bf3..bd035df 100644 --- a/src/gray.f +++ b/src/gray.f @@ -168,6 +168,7 @@ c if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES' if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm + call print_tox end if c c compute power and current density profiles for all rays @@ -282,6 +283,7 @@ c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma if(index_rt.eq.1) then if(j.eq.1) then psinv11=psinv + call calc_tox(i) if(ipolc.eq.0.and.iov(j,k).eq.1) then call pol_limit(qqin,uuin,vvin) ipolc=1 @@ -370,6 +372,77 @@ c end c c +c + subroutine calc_tox(i) + implicit real*8 (a-h,o-z) +c computation of O-X transmission function +c calculated only for xg=(wpl/w)**2>0.5 +c values stored at maximum xg and at maximum transmission + parameter(pi=3.14159265358979d0) + dimension derxg(3),deryg(3),xv(3),anv(3) +c + common/xgxg/xg + common/ygyg/yg + common/xv/xv + common/anv/anv + common/nplr/anpl,anpr + common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/nprw/anprr,anpri + common/parwv/ak0,akinv,fhz + common/dxgyg/derxg,deryg + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx + + anplopt=sqrt(yg/(yg+1.0d0)) + thkb=1.8d2/pi*acos(anpl/sqrt(an2)) + dlen=sqrt(derxg(1)**2+derxg(2)**2+derxg(3)**2) + if (dlen.gt.0.0d0) then + dlen=xg/dlen + tox=exp(-pi*ak0*dlen*sqrt(0.5d0*yg)* + . (2*(1.0d0+yg)*(anplopt-abs(anpl))**2+anpr**2)) + else + tox=0.0d0 + end if + if (xg.gt.0.5d0.and.xg.le.1.0d0) then + if (xg.gt.xgmax) then + ixmx=i + toxxmx=tox + xgmax=xg + ygxmx=yg + anplxmx=anpl + anprxmx=anpr + thkbxmx=thkb + end if +! if (tox.gt.toxmax) then +! itmx=i +! toxmax=tox +! xgtmx=xg +! ygtmx=yg +! anpltmx=anpl +! anprtmx=anpr +! thkbtmx=thkb +! end if +! write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg + end if +c + return + 98 format(1x,i8,30(1x,e16.8e3)) + end +c +c +c + subroutine print_tox + implicit real*8 (a-h,o-z) + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx + write(44,98) ixmx,toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx +! write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx +c + return + 98 format(1x,i8,30(1x,e16.8e3)) + end +c +c c subroutine print_output(i,j,k) implicit real*8 (a-h,o-z) @@ -518,8 +591,8 @@ c write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' - write(8,*) ' #istep j k xt yt zt rt psin' - write(9,*) ' #istep j k xt yt zt rt psin' + write(8,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn' + write(9,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn' write(82,*) ' #istep j k xt yt zwspl zwparab srspl,etre,etim' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' @@ -528,6 +601,7 @@ c write(7,*)'#Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' + write(44,*) '#i Tox X Y Npl Nperp thetakB' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else @@ -852,7 +926,11 @@ c versus psi, rhop, rhot end if if (iequil.eq.1) call surf_anal - +c +c print psi,X,Y on the R,z grid +c + call print_2d +c if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 @@ -1584,6 +1662,34 @@ c end c c +c + subroutine print_2d + implicit real*8 (a-h,o-z) + parameter(nnw=501,nnh=501) + dimension rv(nnw),zv(nnh),psin(nnw,nnh) + dimension btotal(nnw,nnh) +c + common/cpsin/rv,zv,psin + common/psival/psinv + common/eqnn/nr,nz,npp,nintp + common/parpl/brr,bphi,bzz,ajphi + common/xgxg/xg + common/ygyg/yg +c + do j=1,nr + rjcm=rv(j)*1.d2 + do k=1,nz + zkcm=zv(k)*1.d2 + call plas_deriv(rjcm,0.d0,zkcm) + write(91,111) rv(j),zv(k),psinv,brr,bphi,bzz,xg,yg + enddo + write(91,*) ' ' + enddo + return + 111 format(12(1x,e12.5)) + end +c +c c subroutine print_prof implicit real*8 (a-h,o-z) @@ -2607,6 +2713,9 @@ c common/ierr/ierr common/istop/istop common/ipol/ipolc +c used only for computation of O-X transmission function + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx c istpr=0 istpl=1 @@ -2615,6 +2724,11 @@ c s11=0 istop=0 ipolc=0 +c + ixmx=0 +! itmx=0 + xgmax=0.0d0 +! toxmax=0.0d0 c return end @@ -2878,6 +2992,7 @@ 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) + dimension vgv(3),vgv11(3) c common/nray/nrayr,nrayth common/dsds/dst @@ -2890,6 +3005,9 @@ c common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/igrad/igrad + + + common/vgv11/vgv,vgv11 c h=dst hh=h*0.5d0 @@ -2930,6 +3048,7 @@ c . +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq)) end do end do + if(j.eq.1) vgv11=vgv end do c call updatepos @@ -2994,11 +3113,12 @@ c implicit real*8 (a-h,o-z) parameter(ndim=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),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) dimension yy11(6,3),derdxv11(3,3) dimension dery0(3) + dimension vgv(3),vgv11(3) c save yy11,derdxv11 @@ -3019,6 +3139,7 @@ c common/anv/anv common/xv/xv common/idst/idst + common/vgv11/vgv,vgv11 c common/iplane/iplane c @@ -5760,6 +5881,7 @@ c parameter(nxmax=2*jmx-1) dimension txs(nxest),tys(nxest) parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax) dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) + dimension vgv(3),vgv11(3) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk @@ -5778,6 +5900,7 @@ c common/dery0/dery0,dery0mod common/iplane/iplane common/gradjk/gri + common/vgv11/vgv,vgv11 c x4m=0.0d0 x3ym=0.0d0 @@ -5947,10 +6070,16 @@ c dirx=dery0(1) c diry=dery0(2) c dirz=dery0(3) c - dirx=ywrk(4,j,k) - diry=ywrk(5,j,k) - dirz=ywrk(6,j,k) + dirxn=ywrk(4,j,k) + diryn=ywrk(5,j,k) + dirzn=ywrk(6,j,k) + dirn=sqrt(dirxn*dirxn+diryn*diryn+dirzn*dirzn) + dirx=vgv11(1) + diry=vgv11(2) + dirz=vgv11(3) dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) +c + vgdotn=dirxn*dirx+diryn*diry+dirzn*dirz c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir @@ -6024,7 +6153,7 @@ c . snth1*snps1*dx/dr+snth1*csps1*dy/dr+csth1*dz/dr if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr)) . write(nfile,111) istep,j,k,xti,yti,zti,rti, - . xdpsi,ydpsi,zdpsi + . xdpsi,ydpsi,zdpsi,psinv11,dir,dirn,vgdotn c x4m=x4m+xti**4 x3ym=x3ym+xti**3*yti