From 91a2e6cf07d6189eeaf5f8fc75c320ad1cb0b007 Mon Sep 17 00:00:00 2001
From: Michele Guerini Rocco <m.guerinirocco2@campus.unimib.it>
Date: Sun, 19 Dec 2021 16:39:19 +0100
Subject: [PATCH] src: implement toggling of output units

This actually implements the --units option
---
 src/gray_core.f90       | 511 ++++++++++++++++++++++------------------
 src/gray_jetto1beam.f90 |  13 -
 src/magsurf_data.f90    |  53 +++--
 src/main.f90            |   4 +
 src/units.f90           |  57 ++++-
 5 files changed, 374 insertions(+), 264 deletions(-)

diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 24ae752..a3d0880 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -1,34 +1,36 @@
+! This modules contains the core GRAY routines
 module gray_core
+
   use const_and_precisions, only : wp_
+
   implicit none
 
 contains
 
   subroutine gray_main(params, data, results, error, rhout)
     use const_and_precisions, only : zero, one, degree, comp_tiny
-    use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
+    use coreprofiles, only : temp, fzeff
     use dispersion,   only : expinit
-    use gray_params,  only : gray_parameters, gray_data, gray_results, print_parameters, &
+    use gray_params,  only : gray_parameters, gray_data, gray_results, &
+                             print_parameters, &
                              iwarm, ipec, istpr0, igrad, headw, headl, ipass
     use beams,        only : xgygcoeff, launchangles2n
     use beamdata,     only : pweight, rayi2jk
-    use equilibrium,  only : unset_eqspl, unset_rhospl, unset_q
     use errcodes,     only : check_err, print_errn, print_errhcd
     use magsurf_data, only : flux_average, dealloc_surfvec
     use beamdata,     only : init_btr, dealloc_beam, nray, nstep, dst
     use pec,          only : pec_init, spec, postproc_profiles, dealloc_pec, &
                              rhop_tab, rhot_tab
-    use limiter,      only : limiter_unset_globals=>unset_globals
     use utils,        only : vmaxmin
     use reflections,  only : inside
-    use multipass,    only : alloc_multipass, dealloc_multipass, initbeam,  &
-                             initmultipass, turnoffray, plasma_in, plasma_out, wall_out
-    use units,        only : ucenr
+    use multipass,    only : alloc_multipass, dealloc_multipass, initbeam,     &
+                             initmultipass, turnoffray, plasma_in, plasma_out, &
+                             wall_out
     use logger,       only : log_info, log_debug
 
     implicit none
 
-    ! Subroutine arguments
+    ! subroutine arguments
     type(gray_parameters), intent(in)  :: params
     type(gray_data),       intent(in)  :: data
     type(gray_results),    intent(out) :: results
@@ -554,8 +556,6 @@ contains
           end if
         end if
 
-        write(ucenr,*) ''
-
         call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
             pins_beam,ip)                                                                  !  *print power and current density profiles for current beam
             
@@ -609,7 +609,7 @@ contains
 
 
   subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
-    use const_and_precisions, only : wp_, zero
+    use const_and_precisions, only : zero
     implicit none
 ! arguments
     real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
@@ -642,7 +642,7 @@ contains
       ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt)
 !   beam tracing initial conditions  igrad=1
 !   !!!!!! check ray tracing initial conditions   igrad=0 !!!!!!
-    use const_and_precisions, only : wp_,zero,one,pi,half,two,degree,ui=>im
+    use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im
     use math, only : catand
     use gray_params, only : idst
     use beamdata, only : nray,nrayr,nrayth,rwmax
@@ -944,7 +944,6 @@ contains
 
   subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
 !     Runge-Kutta integrator
-    use const_and_precisions, only : wp_
 !    use gray_params, only : igrad
     use beamdata, only : h,hh,h6
     implicit none
@@ -980,7 +979,6 @@ contains
   subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad)
 ! Compute right-hand side terms of the ray equations (dery)
 ! used in R-K integrator
-    use const_and_precisions, only : wp_
     implicit none
 ! arguments
     real(wp_), dimension(6), intent(in) :: y
@@ -1049,7 +1047,7 @@ contains
 
 
   subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri)
-    use const_and_precisions, only : wp_,zero,half
+    use const_and_precisions, only : zero,half
     use beamdata, only : nray,nrayr,nrayth,twodr2
     implicit none
     real(wp_), intent(in) :: ak0
@@ -1199,7 +1197,6 @@ contains
 !     input vectors : dxv1, dxv2, dxv3, dff
 !     output vector : dgg
 !     dff=(1,0,0)
-    use const_and_precisions, only : wp_
     implicit none
 ! arguments
     real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
@@ -1220,7 +1217,6 @@ contains
 
   subroutine solg3(dxv1,dxv2,dxv3,dff,dgg)
 !      rhs "matrix" dff, result in dgg
-    use const_and_precisions, only : wp_
     implicit none
 ! arguments
     real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
@@ -1252,7 +1248,7 @@ contains
 
   subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, &
      xg,yg,derxg,deryg,ajphi)
-    use const_and_precisions, only : wp_,zero,ccj=>mu0inv
+    use const_and_precisions, only : zero,ccj=>mu0inv
     use gray_params, only : iequil
     use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi
     use coreprofiles, only : density
@@ -1406,7 +1402,7 @@ contains
 
   subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
      dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
-    use const_and_precisions, only : wp_,zero,one,half,two
+    use const_and_precisions, only : zero,one,half,two
     use gray_params, only : idst
     implicit none
 ! arguments
@@ -1533,7 +1529,7 @@ contains
 
   subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
      sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error)
-    use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_
+    use const_and_precisions, only : zero,pi,mc2=>mc2_
     use gray_params, only : iwarm,ilarm,ieccd,imx
     use coreprofiles, only : fzeff
     use equilibrium, only : sgnbphi
@@ -1623,7 +1619,7 @@ contains
 
 
   subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0)
-    use const_and_precisions, only : wp_,degree,zero,one,half,im
+    use const_and_precisions, only : degree,zero,one,half,im
     use beamdata, only : nray,nrayth
     use equilibrium, only : bfield
     use gray_params, only : ipol
@@ -1683,10 +1679,8 @@ contains
 
 
   subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
-    use const_and_precisions, only : wp_
 ! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
 !                   (based on an older code)
-    use const_and_precisions, only : wp_
     implicit none
 ! arguments
     integer, intent(in) :: nr,nz
@@ -1897,75 +1891,95 @@ bb: do
 
 
   subroutine print_headers(strheader)
-    use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm
+    ! Prints the headers of all gray_main output tables
+
+    use units, only : uprj0, uwbm, udisp, ucenr, uoutr, upec, usumm, &
+                      unit_active
+
     implicit none
+
     ! subroutine arguments
     character(len=*), dimension(:), intent(in) :: strheader
+
     ! local variables
-    integer :: i,l
-    
-    l=size(strheader)
-    do i=1,l
-      write(uprj0,'(1x,a)') strheader(i)
-      write(uprj0+1,'(1x,a)') strheader(i)
-      write(uwbm,'(1x,a)') strheader(i)
-      write(udisp,'(1x,a)') strheader(i)
-      write(ucenr,'(1x,a)') strheader(i)
-      write(uoutr,'(1x,a)') strheader(i)
-      write(upec,'(1x,a)') strheader(i)
-      write(usumm,'(1x,a)') strheader(i)
+    integer :: i, j
+    integer        :: main_units(8)
+    character(256) :: main_headers(8)
+
+    main_units = [uprj0, uprj0+1, uwbm, udisp, ucenr, uoutr, upec, usumm]
+
+    main_headers(1) = '#sst j k xt yt zt rt'
+    main_headers(2) = '#sst j k xt yt zt rt'
+    main_headers(3) = '#sst w1 w2'
+    main_headers(4) = '#sst Dr_Nr Di_Nr'
+    main_headers(5) = '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp ' &
+                      // 'Npl Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax ' &
+                      // 'iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
+    main_headers(6) = '#i k sst x y R z psin tau Npl alpha index_rt'
+    main_headers(7) = '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
+    main_headers(8) = '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav '  &
+                      // 'drhotjava drhotpav ratjamx ratjbmx stmx psipol ' &
+                      // 'chipol index_rt Jphimx dPdVmx drhotj drhotp P0 ' &
+                      // 'cplO cplX'
+
+    do i=1,size(main_units)
+      if (unit_active(main_units(i))) then
+        do j=1,size(strheader)
+          write (main_units(i), '(1x,a)') strheader(j)
+        end do
+        write (main_units(i), '(1x,a)') trim(main_headers(i))
+      end if
     end do
-    
-    write(uprj0,'(1x,a)')   '#sst j k xt yt zt rt'
-    write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt'
-    write(uwbm,'(1x,a)')    '#sst w1 w2'
-    write(udisp,'(1x,a)')   '#sst Dr_Nr Di_Nr'
-    write(ucenr,'(1x,a)')   '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp Npl '// &
-         'Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
-    write(uoutr,'(1x,a)')   '#i k sst x y R z psin tau Npl alpha index_rt'
-    write(upec,'(1x,a)')    '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
-    write(usumm,'(1x,a)')   '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
-        'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
-        'Jphimx dPdVmx drhotj drhotp P0 cplO cplX'
+
   end subroutine print_headers
 
 
-
   subroutine print_prof
-    use const_and_precisions, only : wp_
-    use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi
-    use coreprofiles, only : density, temp
-    use units, only : uprfin
+    ! Prints the (input) plasma kinetic profiles (unit 98)
+
+    use equilibrium,          only : psinr, nq, fq, frhotor, tor_curr_psi
+    use coreprofiles,         only : density, temp
+    use units,                only : uprfin, unit_active
+
     implicit none
-! local constants
-    real(wp_), parameter :: eps=1.e-4_wp_
-! local variables
+
+    ! local constants
+    real(wp_), parameter :: eps = 1.e-4_wp_
+
+    ! local variables
     integer :: i
-    real(wp_) :: psin,rhot,ajphi,dens,ddens
+    real(wp_) :: psin, rhot, ajphi, dens, ddens
 
-    write(uprfin,*)  ' #psi rhot ne Te q Jphi'
+    if (.not. unit_active(uprfin)) return
+
+    write (uprfin, *) '#psi rhot ne Te q Jphi'
     do i=1,nq
-      psin=psinr(i)
-      rhot=frhotor(sqrt(psin))
-      call density(psin,dens,ddens)
-      call tor_curr_psi(max(eps,psin),ajphi)
+      psin = psinr(i)
+      rhot = frhotor(sqrt(psin))
+      call density(psin, dens, ddens)
+      call tor_curr_psi(max(eps, psin), ajphi)
 
-      write(uprfin,"(12(1x,e12.5))") psin,rhot,dens,temp(psin),fq(psin),ajphi*1.e-6_wp_
+      write (uprfin, '(12(1x,e12.5))') &
+        psin, rhot, dens, temp(psin), fq(psin), ajphi*1.e-6_wp_
     end do
   end subroutine print_prof
 
 
-
   subroutine print_bres(bres)
-    use const_and_precisions, only : wp_
-    use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq
-    use units, only : ubres
+    ! Prints the EC resonance surface table (unit 70)
+
+    use equilibrium,          only : rmnm, rmxm, zmnm, zmxm, bfield, nq
+    use units,                only : ubres, unit_active
+
     implicit none
-! arguments
-    real(wp_) :: bres
-! local constants
-    integer, parameter :: icmx=2002
-! local variables
+
+    ! subroutine arguments
+    real(wp_), intent(in) :: bres
+
+    ! local constants
+    integer, parameter :: icmx = 2002
+
+    ! local variables
     integer :: j,k,n,nconts,nctot
     integer, dimension(10) :: ncpts
     real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
@@ -1973,98 +1987,109 @@ bb: do
     real(wp_) :: rv(nq), zv(nq)
     real(wp_), dimension(nq,nq) :: btotal
 
+    if (.not. unit_active(ubres)) return
     
-    dr = (rmxm-rmnm)/(nq-1)
-    dz = (zmxm-zmnm)/(nq-1)
+    dr = (rmxm - rmnm)/(nq - 1)
+    dz = (zmxm - zmnm)/(nq - 1)
     do j=1,nq
-      rv(j) = rmnm + dr*(j-1)
-      zv(j) = zmnm + dz*(j-1)
+      rv(j) = rmnm + dr*(j - 1)
+      zv(j) = zmnm + dz*(j - 1)
     end do  
 
-!   Btotal on psi grid
-    btmx=-1.0e30_wp_
-    btmn=1.0e30_wp_
+    ! Btotal on psi grid
+    btmx = -1.0e30_wp_
+    btmn = 1.0e30_wp_
     do k=1,nq
-      zzk=zv(k)
+      zzk = zv(k)
       do j=1,nq
-        rrj=rv(j)
-        call bfield(rrj,zzk,bbphi,bbr,bbz)
-        btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
-        if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
-        if(btotal(j,k).le.btmn) btmn=btotal(j,k)
-      enddo
-    enddo
+        rrj = rv(j)
+        call bfield(rrj, zzk, bbphi, bbr, bbz)
+        btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
+        if(btotal(j,k) >= btmx) btmx = btotal(j,k)
+        if(btotal(j,k) <= btmn) btmn = btotal(j,k)
+      end do
+    end do
 
-!   compute Btot=Bres/n with n=1,5
-    write(ubres,*)'#i Btot  R  z'
+    ! compute Btot=Bres/n with n=1,5
+    write (ubres, *) '#i Btot  R  z'
     do n=1,5
-      bbb=bres/dble(n)
-      if (bbb.ge.btmn.and.bbb.le.btmx) then
-        nconts=size(ncpts)
-        nctot=size(rrcb)
-        call cniteq(rv,zv,btotal,nq,nq,bbb,nconts,ncpts,nctot,rrcb,zzcb)
+      bbb = bres/dble(n)
+      if (bbb >= btmn .and. bbb <= btmx) then
+        nconts = size(ncpts)
+        nctot = size(rrcb)
+        call cniteq(rv, zv, btotal, nq, nq, bbb, &
+                    nconts, ncpts, nctot, rrcb, zzcb)
         do j=1,nctot
-          write(ubres,'(i6,12(1x,e12.5))') j,bbb,rrcb(j),zzcb(j)
+          write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j)
         end do
       end if
-      write(ubres,*)
+      write (ubres, *)
     end do
+
   end subroutine print_bres
 
 
+  subroutine print_maps(bres, xgcn, r0, anpl0)
+    ! Prints several input quantities on a regular grid (unit 72)
 
-  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 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
+    use units,        only : umaps, unit_active
+
     implicit none
-! arguments
+
+    ! subroutine arguments
     real(wp_), intent(in) :: bres,xgcn,r0,anpl0
-! local variables
+
+    ! local variables
     integer :: j,k
-    real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl
+    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)
+    if (.not. unit_active(umaps)) return
+
+    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)
+      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'
+    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
+      rj = r(j)
+      anpl = anpl0 * r0/rj
       do k=1,nq
-        zk=z(k)
+        zk = z(k)
         if (iequil < 2) then
-          call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br)
+          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)
+          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)
+        btot =  sqrt(br**2 + bphi**2 + bz**2)
         yg   =  btot/bres
         te   =  temp(psin)
-        call density(psin,ne,dne)
+        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
+        write (umaps,'(12(x,e12.5))') &
+          rj, zk, psin, br, bphi, bz, btot, ne, te, xg, yg, anpl
+      end do
+      write (umaps,*)
+    end do
 
   end subroutine print_maps
 
 
   subroutine print_surfq(qval)
+    ! Print constant ψ surfaces for a given `q` value
+
     use equilibrium,  only : psinr, nq, fq, frhotor, &
                              rmaxis, zmaxis, zbsup, zbinf
     use magsurf_data, only : contours_psi, npoints, print_contour
@@ -2094,16 +2119,16 @@ bb: do
     do i=1,size(qval)
       ! FIXME: check for non monotonous q profile
       call locate(abs(qpsi),nq,qval(i),i1)
-      if (i1>0.and.i1<nq) then
-        call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1), &
-          qval(i),psival)
-        rup=rmaxis
-        rlw=rmaxis
-        zup=(zbsup+zmaxis)/2.0_wp_
-        zlw=(zmaxis+zbinf)/2.0_wp_
-        call contours_psi(psival,rcn,zcn,rup,zup,rlw,zlw)
-        call print_contour(psival,rcn,zcn)
-        rhot=frhotor(sqrt(psival))
+      if (i1>0 .and. i1<nq) then
+        call intlin(abs(qpsi(i1)), psinr(i1), abs(qpsi(i1+1)), psinr(i1+1),  &
+                    qval(i),psival)
+        rup = rmaxis
+        rlw = rmaxis
+        zup = (zbsup + zmaxis)/2.0_wp_
+        zlw = (zmaxis + zbinf)/2.0_wp_
+        call contours_psi(psival, rcn, zcn, rup, zup, rlw, zlw)
+        call print_contour(psival, rcn, zcn)
+        rhot = frhotor(sqrt(psival))
         write (msg, '(4(x,a,"=",g0.3))') &
           'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot
         call log_info(msg, mod='gray_core', proc='print_surfq')
@@ -2113,36 +2138,44 @@ bb: do
   end subroutine print_surfq
 
 
-  subroutine print_projxyzt(stv,ywrk,iproj)
-    use const_and_precisions, only : wp_, comp_huge, zero, one
-    use beamdata, only : nray, nrayr, nrayth, rayi2jk
-    use units, only : uprj0,uwbm
+  subroutine print_projxyzt(stv, ywrk, iproj)
+    ! Prints the beam shape tables (unit 8, 12)
+
+    use const_and_precisions, only : comp_huge, zero, one
+    use beamdata,             only : nray, nrayr, nrayth, rayi2jk
+    use units,                only : uprj0, uwbm, unit_active
+
     implicit none
-! arguments
-    real(wp_), dimension(:), intent(in) :: stv
+
+    ! subroutine arguments
+    real(wp_), dimension(:),   intent(in) :: stv
     real(wp_), dimension(:,:), intent(in) :: ywrk
-    integer, intent(in) :: iproj
-! local variables
-    integer :: jk,jkz,uprj
+    integer,                   intent(in) :: iproj
+
+    ! local variables
+    integer :: jk, jkz, uprj
     integer, dimension(2) ::jkv
-    real(wp_), dimension(3) :: xv1,dir,dxv
-    real(wp_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti
-! common/external functions/variables
+    real(wp_), dimension(3) :: xv1, dir, dxv
+    real(wp_) :: rtimn, rtimx, csth1, snth1, &
+                 csps1, snps1, xti, yti, zti, rti
+
+    if (.not.(unit_active(uprj0) .or. unit_active(uwbm))) return
+
+    ! common/external functions/variables
 
     uprj = uprj0 + iproj
 
-    xv1  = ywrk(1:3,1)
-    dir  = ywrk(4:6,1)
-    dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2)
-    dir  = dir/dirm
+    xv1  = ywrk(1:3, 1)
+    dir  = ywrk(4:6, 1)
+    dir  = dir/norm2(dir)
     csth1 = dir(3)
     snth1 = sqrt(one - csth1**2)
     if(snth1 > zero) then
-      csps1=dir(2)/snth1
-      snps1=dir(1)/snth1
+      csps1 = dir(2)/snth1
+      snps1 = dir(1)/snth1
     else
-      csps1=one
-      snps1=zero
+      csps1 = one
+      snps1 = zero
     end if
 
     if(iproj==0) then
@@ -2158,116 +2191,148 @@ bb: do
       xti = dxv(1)*csps1 - dxv(2)*snps1
       yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1
       zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1
-      rti = sqrt(xti**2 + yti**2)
+      rti = hypot(xti, yti)
 
       jkv=rayi2jk(jk)
-      if(.not.(iproj==0 .and. jk==1))                  &
-         write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') stv(jk),jkv,xti,yti,zti,rti
-      if(iproj==1 .and. jkv(2)==nrayth) write(uprj,*)
+      if(.not.(iproj==0 .and. jk==1)) &
+        write(uprj, '(1x,e16.8e3,2i5,4(1x,e16.8e3))') &
+          stv(jk), jkv, xti, yti, zti, rti
+      if(iproj==1 .and. jkv(2)==nrayth) write (uprj, *)
 
       if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti
       if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti
     end do
-    write(uprj,*)
-    write(uwbm,'(3(1x,e16.8e3))') stv(1),rtimn,rtimx
+    write (uprj, *)
+    write (uwbm, '(3(1x,e16.8e3))') stv(1), rtimn, rtimx
+
   end subroutine print_projxyzt
 
 
+  subroutine print_output(i, jk, st, qj, xv, psinv, btot, bv, ak0, &
+                          anpl, anpr, anv, anprim, dens, tekev,    &
+                          alpha, tau, dids, nhm, nhf, iokhawa,     &
+                          index_rt, ddr, ddi, xg, yg, derxg)
+    ! Prints the ray variables and Hamiltonian error (units 4, 33, 17)
 
-  subroutine print_output(i,jk,st,qj,xv,psinv,btot,bv,ak0,anpl,anpr,anv, &
-           anprim,dens,tekev,alpha,tau,dids,nhm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg)
     use const_and_precisions, only : degree,zero,one
-    use equilibrium, only : frhotor
-    use gray_params, only : istpl0
-    use beamdata, only : nray,nrayth,jkray1
-    use units, only : ucenr,uoutr,udisp
+    use equilibrium,          only : frhotor
+    use gray_params,          only : istpl0
+    use beamdata,             only : nray,nrayth,jkray1
+    use units,                only : ucenr, uoutr, udisp, unit_active
+
     implicit none
-! arguments
-    integer, intent(in) :: i,jk,nhm,nhf,iokhawa,index_rt
-    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
-    real(wp_), intent(in) :: xg,yg
+
+    ! subroutine arguments
+    integer,  intent(in) :: i, jk, nhm, nhf, iokhawa, index_rt
+    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
+    real(wp_),               intent(in) :: xg, yg
     real(wp_), dimension(3), intent(in) :: derxg
 
-! local variables
-    real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn
+    ! local variables
+    real(wp_) :: stm, xxm, yym, zzm, rrm, phideg, rhot, akim, pt, didsn
     integer :: k
 
-    stm=st*1.0e-2_wp_
-    xxm=xv(1)*1.0e-2_wp_
-    yym=xv(2)*1.0e-2_wp_
-    zzm=xv(3)*1.0e-2_wp_
-    rrm=sqrt(xxm**2 + yym**2)
+    stm = st    * 1.0e-2_wp_
+    xxm = xv(1) * 1.0e-2_wp_
+    yym = xv(2) * 1.0e-2_wp_
+    zzm = xv(3) * 1.0e-2_wp_
+    rrm = hypot(xxm, yym)
 
-!     print central ray trajectory. dIds in A/m/W, ki in m^-1
-    if(jk.eq.1) then
-      phideg=atan2(yym,xxm)/degree
+    ! print central ray trajectory. dIds in A/m/W, ki in m⁻¹
+    if(jk == 1) then
+      phideg = atan2(yym, xxm)/degree
       if(psinv>=zero .and. psinv<=one) then
-        rhot=frhotor(sqrt(psinv))
+        rhot = frhotor(sqrt(psinv))
       else
-        rhot=1.0_wp_
+        rhot = 1.0_wp_
       end if
-      akim=anprim*ak0*1.0e2_wp_
-      pt=qj*exp(-tau)
-      didsn=dids*1.0e2_wp_/qj
+      akim  = anprim * ak0 * 1.0e2_wp_
+      pt    = qj*exp(-tau)
+      didsn = dids * 1.0e2_wp_/qj
 
-      write(ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') stm,rrm,zzm,phideg,       &
-         psinv,rhot,dens,tekev,btot,bv,anpr,anpl,anv,akim,alpha,tau,pt,didsn,  &
-         nhm,nhf,iokhawa,index_rt,ddr,xg,yg,derxg
+      ! print central ray data
+      if (unit_active(ucenr)) then
+        write (ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))')  &
+           stm, rrm, zzm, phideg, psinv, rhot, dens, tekev, &
+           btot, bv, anpr, anpl, anv, akim, alpha, tau, pt, &
+           didsn, nhm, nhf, iokhawa, index_rt, ddr, xg, yg, &
+           derxg
+      end if
     end if
 
-!     print conservation of dispersion relation
-    if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
+    ! print conservation of dispersion relation (Hamiltonian error)
+    if(unit_active(udisp) .and. jk==nray) then
+      write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
+    end if
 
-!     print outer trajectories
-    if(mod(i,istpl0)==0)  then
+    ! print outer ray trajectories
+    if(mod(i, istpl0) == 0)  then
       k = jk - jkray1 + 1
-      if(k>0 .and. k<=nrayth) then
-        write(uoutr,'(2i5,9(1x,e16.8e3),i5)') i,k,stm,xxm,yym,rrm,zzm, &
-           psinv,tau,anpl,alpha,index_rt
+      if(unit_active(uoutr) .and. k>0 .and. k<=nrayth) then
+        write(uoutr, '(2i5,9(1x,e16.8e3),i5)') &
+          i, k, stm, xxm, yym, rrm, zzm,  &
+          psinv, tau, anpl, alpha, index_rt
       end if
     end if
   end subroutine print_output
 
 
+  subroutine print_pec(rhop_tab, rhot_tab, jphi, jcd, &
+                       dpdv, currins, pins, index_rt)
+    ! Prints the ECRH and CD profiles (unit 48)
+
+    use units,                only : upec, unit_active
 
-  subroutine print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
-    use const_and_precisions, only : wp_
-    use units, only : upec
     implicit none
-! arguments
-    real(wp_), dimension(:), intent(in) :: rhop_tab,rhot_tab,jphi,jcd,dpdv, &
-                                           currins,pins
-    integer, intent(in) :: index_rt
-! local variables
+
+    ! subroutine arguments
+    real(wp_), dimension(:), intent(in) :: rhop_tab, rhot_tab, jphi, jcd, &
+                                           dpdv, currins, pins
+    integer,                 intent(in) :: index_rt
+
+    ! local variables
     integer :: i
-    
-    do i=1,size(rhop_tab)
-      write(upec,'(7(1x,e16.8e3),i5)') rhop_tab(i),rhot_tab(i), &
-          jphi(i),jcd(i),dpdv(i),currins(i),pins(i),index_rt
+
+    if (.not. unit_active(upec)) return
+
+    do i=1, size(rhop_tab)
+      write(upec, '(7(1x, e16.8e3), i5)') rhop_tab(i), rhot_tab(i),  &
+          jphi(i), jcd(i), dpdv(i), currins(i), pins(i), index_rt
     end do  
-    write(upec,*) ''
+    write(upec, *) ''
   end subroutine print_pec
 
 
+  subroutine print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, &
+                          drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp,      &
+                          drhotj, ratjamx, ratjbmx,  stmx, psipol, chipol,      &
+                          index_rt, p0, cpl1, cpl2)
+    ! Prints global data and final results (unit 7)
+
+    use units, only : ucenr, usumm, unit_active
 
-  subroutine print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav,  &
-        drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
-        stmx,psipol,chipol,index_rt,p0,cpl1,cpl2)
-    use const_and_precisions, only : wp_
-    use units, only : usumm
     implicit none
-    real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav,  &
-        drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx,    &
-        stmx,psipol,chipol,p0,cpl1,cpl2
+
+    ! subroutine arguments
+    real(wp_),  intent(in) :: pabs, icd, dpdvp, jphip, rhotpav, rhotjava, &
+                              drhotpav, drhotjava, dpdvmx, jphimx, rhotp, &
+                              rhotj, drhotp, drhotj, ratjamx, ratjbmx,    &
+                              stmx, psipol, chipol, p0, cpl1, cpl2
     integer, intent(in) :: index_rt
-    
-    write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') icd,pabs,jphip,dpdvp,    &
-         rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
-         stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp,p0,      &
-         cpl1,cpl2
-!    write(usumm,*) ''
+
+    if (unit_active(ucenr)) then
+      write(ucenr, *) ''
+    end if
+
+    if (unit_active(usumm)) then
+      write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') &
+           icd, pabs, jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
+           drhotjava, drhotpav, ratjamx, ratjbmx, stmx, psipol,      &
+           chipol, index_rt, jphimx, dpdvmx, drhotj, drhotp,         &
+           p0, cpl1, cpl2
+    end if
   end subroutine print_finals
 
 end module gray_core
diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90
index 8e6f3d4..4d5475a 100644
--- a/src/gray_jetto1beam.f90
+++ b/src/gray_jetto1beam.f90
@@ -127,19 +127,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
     call unset_prfspl
   end block free_memory
 
-  ! Close output units used for debugging
-  close_debug: block
-    use units, only : ucenr, usumm, uprj0, uwbm, udisp, ubres, &
-                      ucnt, uoutr, ueq, uprfin, uflx, upec
-    integer :: i, debug_units(13)
-
-    debug_units = [ucenr, usumm, uprj0, uprj0+1, uwbm, udisp, &
-                   ubres, ucnt, uoutr, ueq, uprfin, uflx, upec]
-    do i=1,size(debug_units)
-      close(debug_units(i))
-    end do
-  end block close_debug
-
   ! Copy over the results
   pabs = res%pabs
   icd  = res%icd
diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90
index 1bbf45d..f6febad 100644
--- a/src/magsurf_data.f90
+++ b/src/magsurf_data.f90
@@ -535,44 +535,57 @@ contains
   end subroutine contours_psi
 
 
+  subroutine print_contour(psin, rc, zc)
+    ! Prints the flux surface contours (unit 71)
 
-  subroutine print_contour(psin,rc,zc)
     use const_and_precisions, only : wp_, comp_tiny
-    use units, only : ucnt
+    use units,                only : ucnt, unit_active
+
     implicit none
-! arguments
+
+    ! subroutine arguments
     real(wp_), intent(in) :: psin
     real(wp_), dimension(:), intent(in) :: rc, zc
-! local variables
-    integer :: npoints,ic
+
+    ! local variables
+    integer :: npoints, ic
+
+    if (.not. unit_active(ucnt)) return
 
     if (psin < comp_tiny) then
-      write(ucnt,*)' #i psin R  z'
+      write(ucnt, *)' #i psin R  z'
     else
-      npoints=size(rc)
+      npoints = size(rc)
       do ic=1,npoints
-        write(ucnt,'(i6,12(1x,e12.5))') ic,psin,rc(ic),zc(ic)
+        write(ucnt, '(i6,12(1x,e12.5))') ic, psin, rc(ic), zc(ic)
       end do
-      write(ucnt,*)
-      write(ucnt,*)
+      write(ucnt, *)
+      write(ucnt, *)
     end if
   end subroutine print_contour
 
 
+  subroutine print_fluxav(psin, rhot, bav, bmx, bmn, area, vol, &
+                          currp, ajphiav, ffc, ratja, ratjb, qq)
+    ! Prints the flux-averaged quantities table (unit 56)
 
-  subroutine print_fluxav(psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb,qq)
     use const_and_precisions, only : wp_, comp_tiny
-    use units, only : uflx
+    use units,                only : uflx, unit_active
+
     implicit none
-! arguments
-    real(wp_), intent(in) :: psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb,qq
-    
+
+    ! subroutine arguments
+    real(wp_), intent(in) :: psin, rhot, bav, bmx, bmn, area, vol, &
+                             currp, ajphiav, ffc, ratja, ratjb, qq
+
+    if (.not. unit_active(uflx)) return
+
     if (psin < comp_tiny) &
-      write(uflx,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb qq'
-    write(uflx,'(20(1x,e12.5))') psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb,qq
+      write(uflx, *)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol' // &
+                    ' |I_pl| <J_phi> fc ratJa ratJb qq'
+    write(uflx, '(20(1x,e12.5))') &
+      psin, rhot, bav, bmx, bmn, area, vol, currp, ajphiav,  &
+      ffc, ratja, ratjb, qq
   end subroutine print_fluxav
 
 end module magsurf_data
diff --git a/src/main.f90 b/src/main.f90
index 26834c5..2355770 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -1,6 +1,7 @@
 program main
   use const_and_precisions, only : wp_, one, zero
   use logger,               only : INFO, ERROR, set_log_level, log_message
+  use units,                only : set_active_units
   use gray_cli,             only : cli_options, parse_cli_options
   use gray_core,            only : gray_main
   use gray_params,          only : gray_parameters, gray_data, gray_results, &
@@ -23,6 +24,9 @@ program main
   if (opts%quiet) opts%verbose = ERROR
   call set_log_level(opts%verbose)
 
+  ! Activate the given output units
+  call set_active_units(opts%units)
+
   ! Load the parameters and also copy them into
   ! global variables exported by the gray_params
   call read_parameters(opts%params_file, params)
diff --git a/src/units.f90 b/src/units.f90
index ce30d71..f4dffca 100644
--- a/src/units.f90
+++ b/src/units.f90
@@ -1,15 +1,56 @@
+! This modules contains the definitions of the unit numbers used for
+! all the GRAY output ifles and a mechanism to toggle them
 module units
-! STANDARD
+
+  ! Unit numbers
+#ifdef JINTRAC
+  ! JINTRAC
+  integer, parameter :: uprm   =602, ubeam =603, uprf  =644, ueq    =644
+  integer, parameter :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
+  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
+  integer, parameter :: udisp  =617, upec  =648, usumm =607, umaps  =632
+  integer, parameter :: uhead  =638
+#else
+  ! STANDARD
   integer, parameter :: uprm   =  2, ubeam = 97, uprf  = 98, ueq    = 99
   integer, parameter :: uprfin = 55, uflx  = 56, ubres = 70, ucnt   = 71
   integer, parameter :: uprj0  =  8, uwbm  = 12, ucenr =  4, uoutr  = 33
   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 :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
-!  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
-!  integer, parameter :: udisp  =617, upec  =648, usumm =607, umaps  =632
-!  integer, parameter :: uhead  =638
+#endif
 
-end module units
\ No newline at end of file
+  ! List of active units
+  integer, allocatable, save :: active_units(:)
+
+contains
+
+  subroutine set_active_units(units)
+    ! Activate the given units
+    !
+    ! All units are inactive by default and no output
+    ! will be directed to them.
+
+    implicit none
+
+    ! subroutine arguments
+    integer, intent(in) :: units(:)
+
+    active_units = units
+
+  end subroutine set_active_units
+
+
+  function unit_active(unit) result(on)
+    ! Checks whether the given `unit` is active
+
+    implicit none
+
+    ! function arguments
+    integer, intent(in)  :: unit
+    logical              :: on
+
+    on = any(active_units == unit)
+
+  end function unit_active
+
+end module units