From e92ff7cee1c115fa0b3526e02cfb3bc54a1e9570 Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Fri, 19 Jun 2015 13:07:49 +0000 Subject: [PATCH] pec subroutine(s) updated --- Makefile | 4 +- src/gray.f | 702 +++++++++++++++++++++++++---------------------------- 2 files changed, 327 insertions(+), 379 deletions(-) diff --git a/Makefile b/Makefile index 077048a..520296b 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ MAINOBJ=gray.o OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \ graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \ interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ - reflections.o simplespline.o utils.o beamdata.o + reflections.o simplespline.o utils.o beamdata.o # Alternative search paths vpath %.f90 src @@ -28,7 +28,7 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \ graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \ interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ - reflections.o simplespline.o utils.o beamdata.o + reflections.o simplespline.o utils.o beamdata.o conical.o: const_and_precisions.o dierckx.o: const_and_precisions.o dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o diff --git a/src/gray.f b/src/gray.f index 6983c6f..c106093 100644 --- a/src/gray.f +++ b/src/gray.f @@ -118,16 +118,26 @@ c subroutine after_gray_integration use const_and_precisions, only : wp_,zero use graydata_flags, only : ibeam,iwarm,iequil,iprof, - . filenmeqq,filenmprf,filenmbm + . nnd,ipec,ieccd,filenmeqq,filenmprf,filenmbm use graydata_anequil, only : dens0,te0 use beamdata, only : nrayr implicit none c local variables - integer :: iproj,nfilp + integer :: iproj,nfilp,i real(wp_) :: currtka,pabs,currt + + real(wp_), dimension(nnd) :: dpdv,ajphiv,pins,currins + real(wp_), dimension(nnd) :: dvol,darea,ratjav,ratjbv,ratjplv + real(wp_), dimension(nnd) :: ajcdav,ajcdbv,ajplv + real(wp_), dimension(nnd) :: rt_in,rhop_tab,rhot_tab + real(wp_), dimension(0:nnd) :: rtabpsi1 + real(wp_) :: rhotpav,drhotpav + real(wp_) :: rhotjav,rhotjava,drhotjav,drhotjava c common/external functions/variables integer :: index_rt real(wp_) :: st,taumn,taumx,pabstot,currtot + +! arguments c common/ss/st common/index_rt/index_rt @@ -164,13 +174,30 @@ c end if c -c compute power and current density profiles for all rays +c compute power and current density profiles for all rays c pabs=pabstot currt=currtot - call pec(pabs,currt) -c + + call pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, + . dvol,darea,ratjav,ratjbv,ratjplv) + + call spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,pins,currins) + + call postproc_profiles(pabs,currt,rhot_tab,dvol,darea,dpdv,ajphiv, + . rhotpav,drhotpav,rhotjava,drhotjava) + + write(*,*) 'rhotpav,drhotpav,rhotjava,drhotjava' + write(*,99) rhotpav,drhotpav,rhotjava,drhotjava + + + do i=1,nnd + write(48,99) rhop_tab(i),rhot_tab(i),ajphiv(i), + . ajphiv(i)*ratjbv(i),dpdv(i),currins(i),pins(i) + end do + return +99 format(16(1x,e16.8e3)) end c c @@ -605,7 +632,7 @@ c . 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// . 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '// . 'drhotp' - write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' + write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins' write(66,*) "# psipol0 chipol0 powrfl" c else @@ -6374,376 +6401,335 @@ c 99 format(i5,12(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end + + + + subroutine pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, + . dvol,darea,ratjav,ratjbv,ratjplv) + use const_and_precisions, only : wp_,zero,one,izero + use graydata_flags, only : nnd + implicit none + integer :: it,ipec + real(wp_) :: drt,rt,rt1,rhop1 + real(wp_) :: ratjai,ratjbi,ratjpli + real(wp_) :: voli0,voli1,dervoli,areai0,areai1 + real(wp_) :: rrii,rbavi,bmxi,bmni,fci + real(wp_), dimension(nnd), intent(in) :: rt_in + real(wp_), dimension(nnd), intent(out) :: rhop_tab,rhot_tab + real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1 + real(wp_), dimension(nnd), intent(out) :: dvol,darea + real(wp_), dimension(nnd), intent(out) :: ratjav,ratjbv,ratjplv + real(wp_), external :: frhotor,frhopol + +c ipec positive build equidistant grid dimension nnd +c ipec negative read input grid c -c -c - subroutine pec(pabs,currt) - use const_and_precisions, only : wp_,pi,one - use numint, only : simpson - use graydata_flags, only : ipec,ieccd,nnd,dst - use utils, only : locatex, locate, intlin - use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv,pdjki +c ipec=+/-1 rho_pol grid +c ipec=+/-2 rho_tor grid + + voli0=zero + areai0=zero + rtabpsi1(0)=zero + + do it=1,nnd + if(ipec > 0) then +c build equidistant radial grid + drt=one/dble(nnd-1) + rt=dble(it-1)*drt + else +c read radial grid from input + rt=rt_in(it) + drt=(rt_in(it+1)-rt)/2.0_wp_ + end if +c radial coordinate of i-(i+1) interval mid point + if(it < nnd) then + rt1=rt+drt/2.0_wp_ + else + rt1=one + end if + if (abs(ipec) == 1) then + rhop_tab(it)=rt + rhot_tab(it)=frhotor(rt*rt) + rhop1=rt1 + else + rhot_tab(it)=rt + rhop_tab(it)=frhopol(rt) + rhop1=frhopol(rt1) + end if +c psi grid at mid points, dimension nnd+1, for use in pec_tab + rtabpsi1(it)=rhop1**2 + + call valpsispl(rhop1,voli1,dervoli,areai1,rrii, + . rbavi,bmxi,bmni,fci,izero) + dvol(it)=abs(voli1-voli0) + darea(it)=abs(areai1-areai0) + voli0=voli1 + areai0=areai1 + + call ratioj(rhop_tab(it),ratjai,ratjbi,ratjpli) + ratjav(it)=ratjai + ratjbv(it)=ratjbi + ratjplv(it)=ratjpli + end do + end subroutine pec_init + + + + subroutine spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv, + . pins,currins) + use const_and_precisions, only : wp_,one,zero,izero + use graydata_flags, only : nnd,ipec,ieccd + use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv implicit none c local constants - integer, parameter :: nndmx=5001,llmx=21 real(wp_), parameter :: rtbc=one c arguments - real(wp_) :: pabs,currt + real(wp_), intent(in) :: pabs,currt + real(wp_), dimension(nstep):: xxi,ypt,yamp + real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv + real(wp_), dimension(nnd), intent(out) :: pins,currins + real(wp_), dimension(nnd), intent(in) :: dvol,darea + real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 c local variables - integer :: ll,it,nd,intp,kkk,i,j,k,ii,is,idecr,ise0,iis,iis1, - . iise0,iise,ind1,ind2,iind,indi,iif1,istmx,i1,ind,itb1 - integer, dimension(llmx) :: isev - real(wp_) :: stf,voli0,areai0,drt,rt,rt1,psit,psit1,voli1, - . dervoli,areai1,rrii,rbavi,bmxi,bmni,fci,ppa1,cci1,ppa2,cci2, - . dppa,didst,rhotpav,rhot2pav,rhotjav,rhotjava,rhot2java,spds, - . sccs,sccsa,drhot2pav,drhotpav,drhot2java,drhotjava,facpds, - . facjs,rhpp,rhpj,dpdvp,ajphip,ratjamx,ratjbmx,ratjplmx,rhopp, - . drhopp,rhotjfi,rhopfi,ajmxfi,drhotjfi,drhopfi,xps,dpdvmx, - . rhotp,drhotp,stmx,pins_02,pins_05,pins_085,xrhot,currtka,rhop, - . pinsr,ratjpli,ratjai,ratjbi,frhotor,frhopol,fdadrhot, - . fdvdrhot,h,psin - real(wp_), dimension(nstep) :: xxi,ypt,yamp - real(wp_), dimension(nndmx) :: ajphiv,dpdv,dvolt,darea,ratjav, - . ratjbv,ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi,rtab,rhotv - real(wp_), dimension(0:nndmx) :: rtab1 -c common/external functions/variables - integer :: istep,index_rt - real(wp_) :: alpha0,beta0,taumn,taumx,pabstot,currtot,psipol, - . chipol -c - common/istep/istep - common/index_rt/index_rt - common/angles/alpha0,beta0 - common/taumnx/taumn,taumx,pabstot,currtot - common/polcof/psipol,chipol -c - stf=istep*dst - nd=nnd - - if(pabstot.gt.0.0_wp_) then - do ll=1,llmx - isev(ll)=0 - end do - - intp=0 - voli0=0.0_wp_ - areai0=0.0_wp_ - rtab1(0)=0.0_wp_ - do it=1,nd - drt=1.0_wp_/dble(nd-1) - rt=dble(it-1)*drt - if(it.lt.nd) then - rt1=rt+drt/2.0_wp_ - else - rt1=rt - end if - rtab(it)=rt - rtab1(it)=rt1 - dpdv(it)=0.0_wp_ - ajphiv(it)=0.0_wp_ - if (ipec.eq.0) then - psit=rt - psit1=rt1 - else - psit=rt**2 - psit1=rt1**2 - end if - rhotv(it)=frhotor(psit) - call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii, - . rbavi,bmxi,bmni,fci,intp) - dvolt(it)=abs(voli1-voli0) - darea(it)=abs(areai1-areai0) - voli0=voli1 - areai0=areai1 - call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli) - ratjav(it)=ratjai - ratjbv(it)=ratjbi - ratjplv(it)=ratjpli - end do - + integer :: i,ii,j,k,kkk + real(wp_) :: spds,sccs,facpds,facjs + real(wp_), dimension(nnd) :: wdpdv,wajphiv + +c calculation of dP and dI over radial grid + dpdv=zero + ajphiv=zero kkk=1 do j=1,nrayr - if(j.gt.1) kkk=nrayth + if(j > 1) kkk=nrayth do k=1,kkk - ise0=0 ii=iiv(j,k) - if (ii.lt.nstep) then - if(psjki(j,k,ii+1).ne.0.0_wp_) ii=ii+1 - end if - idecr=-1 - is=1 + if (ii < nstep .and. psjki(j,k,ii+1) /= zero) ii=ii+1 + xxi=zero + ypt=zero + yamp=zero do i=1,ii - if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i)) - if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i))) - if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then + xxi(i)=abs(psjki(j,k,i)) + if(xxi(i) <= one) then ypt(i)=ppabs(j,k,i) yamp(i)=ccci(j,k,i) - else - ypt(i)=0.0_wp_ - yamp(i)=0.0_wp_ - end if - if(ise0.eq.0) then - if(xxi(i).lt.rtbc) then - ise0=i - isev(is)=i-1 - is=is+1 - end if - else - if (idecr.eq.-1) then - if(xxi(i).gt.xxi(i-1)) then - isev(is)=i-1 - is=is+1 - idecr=1 - end if - else - if(xxi(i).gt.rtbc) exit - if(xxi(i).lt.xxi(i-1)) then -c isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 - isev(is)=i-1 - is=is+1 - idecr=-1 - end if - end if end if end do -c - isev(is)=i-1 - ppa1=0.0_wp_ - cci1=0.0_wp_ - do iis=1,is-1 - iis1=iis+1 - iise0=isev(iis) - iise=isev(iis1) - if (mod(iis,2).ne.0) then - idecr=-1 - ind1=nd - ind2=2 - iind=-1 - else - idecr=1 - ind1=1 - ind2=nd - iind=1 - end if - do ind=ind1,ind2,iind !!!!!!!!!! is it safe? -c iind=ind !!!!!!!!!! iind reused in the loop! - indi=ind !!!!!!!!!! iind reused in the loop! - if (idecr.eq.-1) indi=ind-1 - rt1=rtab1(indi) - call locatex(xxi,iise,iise0,iise,rt1,itb1) - if(itb1.ge.iise0.and.itb1.lt.iise) then - call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1), - . ypt(itb1+1),rt1,ppa2) - call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1), - . yamp(itb1+1),rt1,cci2) - dppa=ppa2-ppa1 - dpdv(ind)=dpdv(ind)+dppa - didst=(cci2-cci1) - ajphiv(ind)=ajphiv(ind)+didst - ppa1=ppa2 - cci1=cci2 - end if - end do - end do + call pec_tab(xxi,ypt,yamp,ii,rtabpsi1, + . wdpdv,wajphiv) + dpdv=dpdv+wdpdv + ajphiv=ajphiv+wajphiv end do end do - - h=1.0_wp_/dble(nd-1) - rhotpav=0.0_wp_ - drhotpav=0.0_wp_ - rhotjav=0.0_wp_ - rhotjava=0.0_wp_ - rhot2java=0.0_wp_ - - fi=dpdv/h - call simpson (nd,h,fi,spds) - fi=rhotv*dpdv/h - call simpson (nd,h,fi,rhotpav) - fi=rhotv*rhotv*dpdv/h - call simpson (nd,h,fi,rhot2pav) - rhotpav=rhotpav/spds - rhot2pav=rhot2pav/spds - - if (ieccd.ne.0) then - fi=ajphiv/h - call simpson (nd,h,fi,sccs) - fi=rhotv*ajphiv/h - call simpson (nd,h,fi,rhotjav) - fi=abs(ajphiv)/h - call simpson (nd,h,fi,sccsa) - fi=rhotv*abs(ajphiv)/h - call simpson (nd,h,fi,rhotjava) - fi=rhotv*rhotv*abs(ajphiv)/h - call simpson (nd,h,fi,rhot2java) - rhotjav=rhotjav/sccs - rhotjava=rhotjava/sccsa - rhot2java=rhot2java/sccsa - end if - -c factor sqrt(8)=2 sqrt(2) to match with full width -c of gaussian profile - drhot2pav=rhot2pav-rhotpav**2 - drhotpav=sqrt(8.0_wp_*drhot2pav) - drhot2java=rhot2java-rhotjava**2 - drhotjava=sqrt(8.0_wp_*drhot2java) - - spds=0.0_wp_ - sccs=0.0_wp_ - do i=1,nd + +c here dpdv is still dP (not divided yet by dV) +c here ajphiv is still dI (not divided yet by dA) + + spds=zero + sccs=zero + do i=1,nnd spds=spds+dpdv(i) sccs=sccs+ajphiv(i) pins(i)=spds currins(i)=sccs - dpdv(i)=dpdv(i)/dvolt(i) + dpdv(i)=dpdv(i)/dvol(i) ajphiv(i)=ajphiv(i)/darea(i) end do - - facpds=1.0_wp_ - facjs=1.0_wp_ - if(spds.gt.0.0_wp_) facpds=pabs/spds - if(sccs.ne.0.0_wp_) facjs=currt/sccs - - do i=1,nd + + facpds=one + facjs=one + if(spds > zero) facpds=pabs/spds + if(sccs /= zero) facjs=currt/sccs + + do i=1,nnd dpdv(i)=facpds*dpdv(i) ajphiv(i)=facjs*ajphiv(i) - ajcdav(i)=ajphiv(i)*ratjav(i) - ajcdbv(i)=ajphiv(i)*ratjbv(i) - ajplv(i)=ajphiv(i)*ratjplv(i) - end do - - rhpp=frhopol(rhotpav) - rhpj=frhopol(rhotjava) - dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) - ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) - call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) + end do - call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, - . drhotp,drhopp) - if(ieccd.ne.0) then - call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, - . drhotjfi,drhopfi) - xps=rhopfi - else - rhotjfi=0.0_wp_ - rhopfi=0.0_wp_ - ajmxfi=0.0_wp_ - ajphip=0.0_wp_ - drhotjfi=0.0_wp_ - drhopfi=0.0_wp_ - xps=rhopp - ratjamx=1.0_wp_ - ratjbmx=1.0_wp_ - ratjplmx=1.0_wp_ - end if - - iif1=iiv(1,1) - istmx=1 - do i=2,iif1 - if(psjki(1,1,i).ge.0.0_wp_) then - if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i - end if - end do - stmx=istmx*dst - - pins_02=0.0_wp_ - pins_05=0.0_wp_ - pins_085=0.0_wp_ - - xrhot=0.2_wp_ - call locate(rhotv,nd,xrhot,i1) - call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), - . xrhot,pins_02) - xrhot=0.5_wp_ - call locate(rhotv,nd,xrhot,i1) - call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), - . xrhot,pins_05) - xrhot=0.85_wp_ - call locate(rhotv,nd,xrhot,i1) - call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), - . xrhot,pins_085) - - else - ajmxfi=0.0_wp_ - ajphip=0.0_wp_ - dpdvp=0.0_wp_ - dpdvmx=0.0_wp_ - rhotjfi=1.0_wp_ - rhotjav=1.0_wp_ - rhotjava=1.0_wp_ - rhotp=1.0_wp_ - rhotpav=1.0_wp_ - drhotjfi=0.0_wp_ - drhotjava=0.0_wp_ - drhotp=0.0_wp_ - drhotpav=0.0_wp_ - ratjamx=1.0_wp_ - ratjbmx=1.0_wp_ - taumn=0 - taumx=0 - stmx=stf - pins_02=0.0_wp_ - pins_05=0.0_wp_ - pins_085=0.0_wp_ -c end of pabstot > 0 - end if - -c dPdV [MW/m^3], Jcd [MA/m^2] - - if(ieccd.eq.0) currt=0.0_wp_ - currtka=currt*1.0e3_wp_ - - write(6,*)' ' - write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '// - .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// - .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' - write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp, - . rhotjfi,rhotjava,rhotp,rhotpav, - . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, - . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp - - write(7,99) currtka,pabstot,ajphip,dpdvp, - . rhotjfi,rhotjava,rhotp,rhotpav, - . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, - . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp - - do i=1,nd - if (ipec.eq.0) then - psin=rtab(i) - rhop=sqrt(rtab(i)) +c now dpdv is dP/dV [MW/m^3] +c now ajphiv is J_phi=dI/dA [MA/m^2] + + end subroutine spec + + + + subroutine pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv) +c Power and current projected on psi grid - mid points + use const_and_precisions, only : wp_,one,zero + use graydata_flags, only : nnd + use beamdata, only : nstep + use utils, only : locatex,intlin +c arguments + integer, intent(in) :: ii + integer, parameter :: llmx = 21 + integer, dimension(llmx) ::isev + real(wp_), dimension(nstep), intent(in) :: xxi,ypt,yamp + real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 + real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv +c local variables + real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1 + integer :: i,is,ise0,idecr,iise0,iise,iis,iis1 + integer :: ind1,ind2,iind,ind,indi,itb1 + + isev=0 + ise0=0 + idecr=-1 + is=1 + wdpdv=zero + wajphiv=zero + do i=1,ii + if(ise0 == 0) then + if(xxi(i) < one) then + ise0=i + isev(is)=i-1 + is=is+1 + end if else - psin=rtab(i)**2 - rhop=rtab(i) + if (idecr == -1) then + if(xxi(i) > xxi(i-1)) then + isev(is)=i-1 + is=is+1 + idecr=1 + end if + else + if(xxi(i) > one) exit + if(xxi(i) < xxi(i-1)) then + isev(is)=i-1 + is=is+1 + idecr=-1 + end if + end if end if - pinsr=0.0_wp_ - if(pabstot.gt.0) pinsr=pins(i)/pabstot - write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), - . currins(i),pins(i),pinsr,real(index_rt) end do + + isev(is)=i-1 + ppa1=zero + cci1=zero + + do iis=1,is-1 + iis1=iis+1 + iise0=isev(iis) + iise=isev(iis1) + if (mod(iis,2) /= 0) then + idecr=-1 + ind1=nnd + ind2=2 + iind=-1 + else + idecr=1 + ind1=1 + ind2=nnd + iind=1 + end if + do ind=ind1,ind2,iind + indi=ind + if (idecr == -1) indi=ind-1 + rt1=rtabpsi1(indi) + call locatex(xxi,iise,iise0,iise,rt1,itb1) + if(itb1 >= iise0 .and. itb1 < iise) then + call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),ypt(itb1+1), + . rt1,ppa2) + call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1), + . rt1,cci2) + dppa=ppa2-ppa1 + wdpdv(ind)=wdpdv(ind)+dppa + didst=cci2-cci1 + wajphiv(ind)=wajphiv(ind)+didst + ppa1=ppa2 + cci1=cci2 + end if + end do + end do + + end subroutine pec_tab + + + + subroutine postproc_profiles(pabs,currt,rhot_tab,dvol,darea, + . dpdv,ajphiv,rhotpav,drhotpav,rhotjava,drhotjava) +c radial average values over power and current density profile + use const_and_precisions, only : wp_,one,zero,izero,pi + use graydata_flags, only : nnd,ieccd + implicit none + real(wp_), intent(in) :: pabs,currt + real(wp_), intent(out) :: rhotpav,rhotjava + real(wp_), intent(out) :: drhotpav,drhotjava + real(wp_) :: rhopjava,rhoppav + real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp + real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi + real(wp_) :: ratjamx,ratjbmx,ratjplmx + real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea + real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv + real(wp_), external :: frhopol,fdadrhot,fdvdrhot + + integer :: i + real(wp_) :: sccsa + real(wp_) :: rhotjav,rhot2pav,rhot2java + + rhotpav=zero + rhot2pav=zero + rhotjav=zero + rhotjava=zero + rhot2java=zero + + if (pabs > zero) then + rhotpav=sum(rhot_tab(1:nnd)*dpdv(1:nnd)*dvol(1:nnd))/pabs + rhot2pav=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* + . dpdv(1:nnd)*dvol(1:nnd))/pabs + end if + + if (ieccd /= 0) then + if (abs(currt) > zero) then + rhotjav=sum(rhot_tab(1:nnd)*ajphiv(1:nnd)*darea(1:nnd))/currt + end if + sccsa=sum(abs(ajphiv(1:nnd))*darea(1:nnd)) + if (sccsa > zero) then + rhotjava=sum(rhot_tab(1:nnd)*abs(ajphiv(1:nnd)) + . *darea(1:nnd))/sccsa + rhot2java=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* + . abs(ajphiv(1:nnd))*darea(1:nnd))/sccsa + end if + end if + +c factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile + drhotpav=sqrt(8._wp_*(rhot2pav-rhotpav**2)) + drhotjava=sqrt(8._wp_*(rhot2java-rhotjava**2)) - return - 99 format(30(1x,e12.5)) - end -c -c -c - subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) + rhoppav=frhopol(rhotpav) + rhopjava=frhopol(rhotjava) + dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhoppav)) + ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhopjava)) + call ratioj(rhopjava,ratjamx,ratjbmx,ratjplmx) + + call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp) + if(ieccd.ne.0) then + call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi) + else + rhotjfi=0.0d0 + ajmxfi=0.0d0 + ajphip=0.0d0 + drhotjfi=0.0d0 + ratjamx=1.0d0 + ratjbmx=1.0d0 + ratjplmx=1.0d0 + end if + + end subroutine postproc_profiles + + + + subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe) use const_and_precisions, only : wp_,emn1 use graydata_flags, only : ipec use utils, only : locatex, locate, intlin, vmaxmini implicit none c arguments integer :: nd - real(wp_) :: rhotmx,rhopmx,ypk,drhot,drhop real(wp_), dimension(nd) :: xx,yy + real(wp_), intent(out) :: xpk,ypk,dxxe c local variables integer :: imn,imx,ipk,ie1,ie2 - real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,xpk,yye,rte1,rte2, - . psie1,psie2,rhopmn,rhote2,rhote1,rhotmn -c common/external functions/variables - real(wp_) :: rhotp,rhotm,ypkp,ypkm - real(wp_) :: frhotor -c - common/jmxmn/rhotp,rhotm,ypkp,ypkm + real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2 + real(wp_) :: ypkp,ypkm c call vmaxmini(yy,nd,ymn,ymx,imn,imx) ypk=0.0_wp_ @@ -6789,47 +6775,9 @@ c call locate(yy,nd,yye,ie2) call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2) end if -c - ipk=1 - if(ymx.ne.0.and.ymn.ne.0) ipk=2 -c - drhop=0.0_wp_ - drhot=0.0_wp_ - psie1=0.0_wp_ - psie2=1.0_wp_ - rhopmx=1.0_wp_ - rhopmn=0.0_wp_ - if (ie1.gt.0.or.ie2.gt.0) then - if(ipec.eq.0) then - rhopmx=sqrt(xpkp) - rhopmn=sqrt(xpkm) - psie2=rte2 - psie1=rte1 - drhop=sqrt(rte2)-sqrt(rte1) - else - rhopmx=xpkp - rhopmn=xpkm - drhop=rte2-rte1 - psie2=rte2**2 - psie1=rte1**2 - end if - end if -c - rhotmx=frhotor(rhopmx**2) - rhotmn=frhotor(rhopmn**2) - rhote2=frhotor(psie2) - rhote1=frhotor(psie1) - drhot=rhote2-rhote1 -c - if(ipk.eq.2) then - drhop=-drhop - drhot=-drhot - end if + dxxe=rte2-rte1 + if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe - ypk=ypkp - rhotp=rhotmx - rhotm=rhotmn -c return end c