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

@@ -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     ipec=0/1/-1 :pec profiles grid equispaced in psi/rhop/as input
-c     nnd       :number of grid steps for pec profiles +1
+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)
+        if (ipec.ge.0) then
-        rt=dble(it-1)*drt
+          drt=1.0d0/dble(nd-1)
-        if(it.lt.nd) then
+          rt=dble(it-1)*drt
-          rt1=rt+drt/2.0d0
+          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
+            do ind=ind1,ind2,iind     !!!!!!!!!! is it safe?
-              iind=ind
+              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)
+      h=1.0d0/dble(nd-1)    !!!!!!!!!! equi-spaced grid expected in
-      rhotpav=0.0d0
+      rhotpav=0.0d0         !!!!!!!!!! simpson subroutine
-      drhotpav=0.0d0
+      drhotpav=0.0d0        !!!!!!!!!! wrong integrals until fixed
-      rhotjav=0.0d0
+      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))

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=dpdv*1.e6_r8
-  dpdv=0._r8
+  jcd=jcd*1.e6_r8
-  jcd=0._r8
-! =================
   ierr=iercom
 end subroutine gray_main