From b02c254a3585ce6987a95c10fa91ea9925b4592e Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 12 Jun 2014 12:51:03 +0000 Subject: [PATCH] added output of dP/dV and Jcd on arbitrary input grid. ijetto=2 (null psi outside plasma boundary) not handled yet --- src/gray-externals.f | 64 +++++++++++++++++++++++++++++++------------- src/gray_main.f90 | 20 +++++++++----- 2 files changed, 59 insertions(+), 25 deletions(-) diff --git a/src/gray-externals.f b/src/gray-externals.f index 8aa8a92..8f2347d 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -37,10 +37,12 @@ c ray integration: end end - subroutine after_gray_integration + subroutine after_gray_integration(rhopin,nrho,dpdvout,ajcdout) implicit real*8 (a-h,o-z) parameter(zero=0.0d0) character*255 filenmeqq,filenmprf,filenmbm + + dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) c common/ss/st common/ibeam/ibeam @@ -92,7 +94,7 @@ c compute power and current density profiles for all rays c pabs=pabstot currt=currtot - call pec(pabs,currt) + call pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) c return 99 format(20(1x,e16.8e3)) @@ -603,10 +605,12 @@ c read(602,*) ieccd if (ieccd.ge.10) call Setup_SpitzFunc c -c ipec=0/1 :pec profiles grid in psi/rhop -c nnd :number of grid steps for pec profiles +1 +c ipec=0/1/-1 :pec profiles grid equispaced in psi/rhop/as input +c nnd :number of grid steps for pec profiles +1 c read(602,*) ipec,nnd +c input radial grid forced on output. nnd is ignored in this case + ipec=-1 c c dst integration step c nstep maximum number of integration steps @@ -5921,11 +5925,12 @@ c c c c - subroutine pec(pabs,currt) + subroutine pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) use itm_constants, only : pi=>itm_pi 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) @@ -5962,7 +5967,12 @@ c common/polcof/psipol,chipol c stf=istep*dst - nd=nnd + if (ipec.ge.0) then + nd=nnd + else + ! ipec=-1 uses rho_pol input grid + nd=nrho + end if if(pabstot.gt.0.0d0) then do ll=1,llmx @@ -5974,12 +5984,23 @@ c areai0=0.0d0 rtab1(0)=0.0d0 do it=1,nd - drt=1.0d0/dble(nd-1) - rt=dble(it-1)*drt - if(it.lt.nd) then - rt1=rt+drt/2.0d0 + if (ipec.ge.0) then + drt=1.0d0/dble(nd-1) + rt=dble(it-1)*drt + if(it.lt.nd) then + rt1=rt+drt/2.0d0 + else + rt1=rt + end if else - rt1=rt +c radial coordinate of i-(i+1) interval left point + rt=rhopin(it) + if(it.lt.nd) then +c radial coordinate of i-(i+1) interval mid point + rt1=(rt+rhopin(it+1))/2.0d0 + else + rt1=rt + end if end if rtab(it)=rt rtab1(it)=rt1 @@ -6050,7 +6071,7 @@ c else if(xxi(i).gt.rtbc) exit if(xxi(i).lt.xxi(i-1)) then - isev(is)=i + isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 is=is+1 idecr=-1 end if @@ -6076,8 +6097,8 @@ c ind2=nd iind=1 end if - do ind=ind1,ind2,iind - iind=ind + do ind=ind1,ind2,iind !!!!!!!!!! is it safe? + iind=ind !!!!!!!!!! iind reused in the loop! if (idecr.eq.-1) iind=ind-1 rt1=rtab1(iind) call locatex(xxi,iise,iise0,iise,rt1,itb1) @@ -6098,10 +6119,10 @@ c end do end do - h=1.0d0/dble(nd-1) - rhotpav=0.0d0 - drhotpav=0.0d0 - rhotjav=0.0d0 + h=1.0d0/dble(nd-1) !!!!!!!!!! equi-spaced grid expected in + rhotpav=0.0d0 !!!!!!!!!! simpson subroutine + drhotpav=0.0d0 !!!!!!!!!! wrong integrals until fixed + rhotjav=0.0d0 !!!!!!!!!! (d)rhotpav,rhotjav,(d)rhotjava affected rhotjava=0.0d0 rhot2java=0.0d0 @@ -6265,6 +6286,13 @@ c dPdV [MW/m^3], Jcd [MA/m^2] write(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), . currins(i),pins(i),pinsr,real(index_rt) end do +c fill output arguments only if ipec=-1 + if (ipec.eq.-1) then + do i=1,nd + dpdvout(i)=dpdv(i) + ajcdout(i)=ajcdbv(i) + end do + end if return 49 format(i5,20(1x,e12.5)) diff --git a/src/gray_main.f90 b/src/gray_main.f90 index 79e7fe2..eae1e39 100644 --- a/src/gray_main.f90 +++ b/src/gray_main.f90 @@ -18,6 +18,7 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & integer :: istop, iercom, igrad, iopmin, iowmin, index_rt, ipass real(r8) :: sox, p0mw, powrfl, taumn, taumx, pabstot, currtot real(r8) :: p0mw1, powtr, pabstott, currtott + real(r8), dimension(nrho) :: rhopol,dpdv1pass, jcd1pass common/istop/istop common/ierr/iercom @@ -49,10 +50,13 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & call gray_integration ! postprocessing - call after_gray_integration + rhopol=sqrt(psijet) + call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass) pabstott=pabstot currtott=currtot + dpdv=dpdv1pass + jcd=jcd1pass powtr=p0mw-pabstot if (iowmin==2 .and. ipass>1) then @@ -67,9 +71,11 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & call paraminit call ic_rt2 call gray_integration - call after_gray_integration + call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass) pabstott=pabstott+pabstot currtott=currtott+currtot + dpdv=dpdv+dpdv1pass + jcd=jcd+jcd1pass index_rt=3 sox=-sox @@ -79,17 +85,17 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & call paraminit call ic_rt2 call gray_integration - call after_gray_integration + call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass) pabstott=pabstott+pabstot currtott=currtott+currtot + dpdv=dpdv+dpdv1pass + jcd=jcd+jcd1pass end if write(*,*) write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8 pec=pabstott*1.e6_r8 icd=currtott*1.e6_r8 -! NOT YET COMPLETE! - dpdv=0._r8 - jcd=0._r8 -! ================= + dpdv=dpdv*1.e6_r8 + jcd=jcd*1.e6_r8 ierr=iercom end subroutine gray_main