some additional output prints (2D maps and Bcart, Ncart along central ray)

src/graycore.f90

@@ -8,7 +8,7 @@ contains
               eqp,psrad,terad,derad,zfc,prfp,  rlim,zlim,                      &
               p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0,           &
               psipol0,chipol0,  dpdv,jcd,pabs,icd,  outp,rtrp,hcdp,ierr, rhout)
-    use const_and_precisions, only : zero, one
+    use const_and_precisions, only : zero, one, degree
     use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
     use dispersion, only : expinit
     use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
@@ -59,7 +59,7 @@ contains
     real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip
     real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
 
-    real(wp_), dimension(3) :: xv,anv0,anv
+    real(wp_), dimension(3) :: xv,anv0,anv,bv
     real(wp_), dimension(:,:), allocatable :: yw,ypw,gri
     real(wp_), dimension(:,:,:), allocatable :: xc,du1,ggri
     integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1
@@ -115,6 +115,8 @@ contains
 ! print ne, Te, q, Jphi versus psi, rhop, rhot
     call print_bres(bres)
     call print_prof
+    call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), &
+      sin(beta0*degree))
 ! ======= pre-proc prints  END  ======
 
 ! ======= main loop BEGIN ======
@@ -152,7 +154,7 @@ contains
         xv = yw(1:3,jk)
         anv = yw(4:6,jk)
         call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
-           psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn)
+           psinv,dens,btot,bv,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn)
 ! update global error code and print message
         if (ierrn/=0) then
           ierr = ior(ierr,ierrn)
@@ -201,8 +203,8 @@ contains
         dids0(jk)=dids
         ccci0(jk)=ccci(jk,i)
 
-        call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr,     &
+        call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,bv,ak0,anpl,anpr,  &
-           anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
+           anv,anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
 
       end do
 
@@ -596,8 +598,9 @@ contains
 
       ddr = anx**2 + any**2 + anz**2 - an20
       ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
-      call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,ak0,zero,zero,zero, &
+      call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,(/zero,zero,zero/), &
-           zero,zero,zero,zero,zero,0,0,1,ddr,ddi) ! st=0, index_rt=1, Btot=0, psin=-1
+           ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, &
+           0,0,1,ddr,ddi) ! st=0, index_rt=1, B=0, N=0, psin=-1
     end do
   end subroutine ic_gb
 
@@ -666,7 +669,7 @@ contains
 
 
   subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
-     xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr)
+     bv,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr)
 ! Compute right-hand side terms of the ray equations (dery)
 ! used after full R-K step and grad(S_I) update
     use errcodes, only : pnpl
@@ -679,10 +682,11 @@ contains
     real(wp_), dimension(6), intent(out) :: dery
     real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr
     real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm
+    real(wp_), dimension(3), intent(out) :: bv
     integer, intent(out) :: ierr
 ! local variables
     real(wp_) :: gr2,ajphi
-    real(wp_), dimension(3) :: dgr2,bv,derxg,deryg
+    real(wp_), dimension(3) :: dgr2,derxg,deryg
     real(wp_), dimension(3,3) :: derbv
     real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
 
@@ -1648,8 +1652,8 @@ bb: do
     write(uprj0+1,*) ' #sst j k xt yt zt rt'
     write(uwbm,*)    ' #sst w1 w2'
     write(udisp,*)   ' #sst Dr_Nr Di_Nr'
-    write(ucenr,*)   ' #sst R z phi psin rhot ne Te Btot Nperp Npl ki '// &
+    write(ucenr,*)   ' #sst R z phi psin rhot ne Te Btot Bx By Bx Nperp Npl '// &
-         'alpha tau Pt dIds nhmax iohkw index_rt ddr'
+         'Nx Ny Nz ki alpha tau Pt dIds nhmax iohkw index_rt ddr'
     write(uoutr,*)   ' #i k sst x y R z psin tau Npl alpha index_rt'
     write(upec,*)    ' #rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
     write(usumm,*)   ' #Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
@@ -1742,6 +1746,59 @@ bb: do
 
 
 
+  subroutine print_maps(bres,xgcn,r0,anpl0)
+    use const_and_precisions, only : wp_
+    use gray_params, only : iequil
+    use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, &
+                            equinum_fpol, nq
+    use coreprofiles, only : density, temp
+    use units, only : umaps
+    implicit none
+! arguments
+    real(wp_), intent(in) :: bres,xgcn,r0,anpl0
+! local variables
+    integer :: j,k
+    real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl
+    real(wp_), dimension(nq) :: r, z
+
+    
+    dr = (rmxm-rmnm)/(nq-1)
+    dz = (zmxm-zmnm)/(nq-1)
+    do j=1,nq
+      r(j) = rmnm + dr*(j-1)
+      z(j) = zmnm + dz*(j-1)
+    end do
+    
+    write(umaps,*)'#R z psin Br Bphi Bz Btot ne Te X Y Npl'
+    do j=1,nq
+      rj=r(j)
+      anpl=anpl0*r0/rj
+      do k=1,nq
+        zk=z(k)
+        if (iequil < 2) then
+          call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br)
+        else
+          call equinum_psi(rj,zk,psinv=psin,dpsidr=bz,dpsidz=br)
+          call equinum_fpol(psin,fpolv=bphi)
+        end if
+        br   = -br/rj
+        bphi =  bphi/rj
+        bz   =  bz/rj
+        btot =  sqrt(br**2+bphi**2+bz**2)
+        yg   =  btot/bres
+        te   =  temp(psin)
+        call density(psin,ne,dne)
+        xg   =  xgcn*ne
+        write(umaps,'(12(x,e12.5))') rj,zk,psin,br,bphi,bz,btot,ne,te,xg,yg,anpl
+      enddo
+      write(umaps,*)
+    enddo
+
+  end subroutine print_maps
+
+
+
+
   subroutine print_surfq(qval)
     use const_and_precisions, only : wp_, one
     use equilibrium, only : psinr,nq,fq,frhotor,rmaxis,zmaxis, &
@@ -1844,8 +1901,8 @@ bb: do
 
 
 
-  subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, &
+  subroutine print_output(i,jk,st,qj,xv,psinv,btot,bv,ak0,anpl,anpr,anv, &
-           dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
+           anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
     use const_and_precisions, only : degree,zero,one
     use equilibrium, only : frhotor
     use gray_params, only : istpl0
@@ -1854,7 +1911,7 @@ bb: do
     implicit none
 ! arguments
     integer, intent(in) :: i,jk,nhf,iokhawa,index_rt
-    real(wp_), dimension(3), intent(in) :: xv
+    real(wp_), dimension(3), intent(in) :: xv,bv,anv
     real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim
     real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi
 ! local variables
@@ -1879,8 +1936,8 @@ bb: do
       pt=exp(-tau)
       didsn=dids*1.0e2_wp_/qj
 
-      write(ucenr,'(16(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, &
+      write(ucenr,'(22(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg,       &
-         psinv,rhot,dens,tekev,btot,anpr,anpl,akim,alpha,tau,pt,didsn,        &
+         psinv,rhot,dens,tekev,btot,bv,anpr,anpl,anv,akim,alpha,tau,pt,didsn,  &
          nhf,iokhawa,index_rt,ddr
     end if
 

src/units.f90

@@ -1,13 +1,15 @@
 module units
 ! STANDARD
-!  integer, parameter :: uprm   =  2, ubeam = 97, uprf  = 98, ueq    = 99
+  integer, parameter :: uprm   =  2, ubeam = 97, uprf  = 98, ueq    = 99
-!  integer, parameter :: uprfin = 55, uflx  = 56, ubres = 70, ucnt   = 71
+  integer, parameter :: uprfin = 55, uflx  = 56, ubres = 70, ucnt   = 71
-!  integer, parameter :: uprj0  =  8, uwbm  = 12, ucenr =  4, uoutr  = 33
+  integer, parameter :: uprj0  =  8, uwbm  = 12, ucenr =  4, uoutr  = 33
-!  integer, parameter :: udisp  = 17, upec  = 48, usumm =  7
+  integer, parameter :: udisp  = 17, upec  = 48, usumm =  7, umaps  = 72
+  integer, parameter :: uhead  = 78
 ! JINTRAC
-  integer, parameter :: uprm   =602, ubeam =603, uprf  =644, ueq    =644
+!  integer, parameter :: uprm   =602, ubeam =603, uprf  =644, ueq    =644
-  integer, parameter :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
+!  integer, parameter :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
-  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
+!  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
-  integer, parameter :: udisp  =617, upec  =648, usumm =607, uhead  =638
+!  integer, parameter :: udisp  =617, upec  =648, usumm =607, umaps  =632
+!  integer, parameter :: uhead  =638
 
 end module units