fix memory errors and use automatic deallocation

src/gray_cli.f90

@@ -24,8 +24,7 @@ module gray_cli
 
   private
   public :: cli_options, print_cli_options, parse_cli_options, &
-            deinit_cli_options, parse_param_overrides, &
+            parse_param_overrides, print_version, get_next_command
-            print_version, get_next_command
 
 contains
 
@@ -285,18 +284,4 @@ contains
     i = i + 1                                ! increment counter
   end subroutine
 
-
-  subroutine deinit_cli_options(opts)
-    ! Frees all memory allocated by the parse_cli_options subroutine
-
-    ! subroutine arguments
-    type(cli_options), intent(inout) :: opts
-
-    if (allocated(opts%output_dir))   deallocate(opts%output_dir)
-    if (allocated(opts%params_file))  deallocate(opts%params_file)
-    if (allocated(opts%config_file))  deallocate(opts%config_file)
-    if (allocated(opts%sum_filelist)) deallocate(opts%sum_filelist)
-    if (allocated(opts%tables))       deallocate(opts%tables)
-  end subroutine deinit_cli_options
-
 end module gray_cli

src/gray_params.f90

@@ -6,7 +6,6 @@ module gray_params
   implicit none
 
   integer, parameter :: lenfnm = 256
-  integer, parameter :: headw = 132, headl = 21
 
   ! Polarisation mode
   enum, bind(c)
@@ -266,106 +265,149 @@ module gray_params
 
 contains
 
-  subroutine write_gray_header(params, file)
+  function make_gray_header(params) result(header)
     ! Writes the Gray output header to a file
 
-    ! subroutine arguments
+    ! function arguments
     type(gray_parameters), intent(in) :: params
-    integer :: file
 
     ! local variables
-    character(len=8)  :: date
+    character(len=200)          :: line
-    character(len=10) :: time
+    character(len=21*len(line)) :: header
+    character(len=8)            :: date
+    character(len=10)           :: time
 #ifndef REVISION
     character(len=*), parameter :: REVISION="unknown"
 #endif
 
+    header = ""
+
     ! Date and time
     call date_and_time(date, time)
-    write(file, '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
+    write(line, '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
                       date(1:4), date(5:6), date(7:8), &
                       time(1:2), time(3:4), time(5:10)
+    call append(header, line)
 
     ! Git revision
-    write(file, '("# GRAY Git revision: ",a)') REVISION
+    write(line, '("# GRAY Git revision: ",a)') REVISION
+    call append(header, line)
 
     ! MHD equilibrium parameters
     if (params%equilibrium%iequil > 0) then
-      write(file, '("# EQL input: ",a)') trim(params%equilibrium%filenm)
+      write(line, '("# EQL input: ",a)') trim(params%equilibrium%filenm)
+      call append(header, line)
       ! TODO add missing values
-      write(file, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
+      write(line, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
         0, 0, 0, 0, 0
+      call append(header, line)
     else
-      write(file, '("# EQL input: N/A (vacuum)")')
+      write(line, '("# EQL input: N/A (vacuum)")')
-      write(file, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
+      call append(header, line)
+      write(line, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
+      call append(header, line)
     end if
 
-    write(file, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
+    write(line, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
       params%equilibrium%iequil, params%equilibrium%sgnb, &
       params%equilibrium%sgni, params%equilibrium%factb
+    call append(header, line)
     if (params%equilibrium%iequil > 1) then
-      write(file, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
+      write(line, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
         params%equilibrium%icocos, params%equilibrium%ipsinorm, &
         params%equilibrium%idesc, params%equilibrium%ifreefmt
-      write(file, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
+      call append(header, line)
+      write(line, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
         params%equilibrium%ssplps, params%equilibrium%ssplf, &
         params%equilibrium%ixp
+      call append(header, line)
     else
-      write(file, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
+      write(line, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
-      write(file, '("# EQL ssplps ssplf ixp: N/A (analytical)")')
+      call append(header, line)
+      write(line, '("# EQL ssplps ssplf ixp: N/A (analytical)")')
+      call append(header, line)
     end if
 
     ! Profiles parameters
     if (params%equilibrium%iequil > 0) then
-      write(file, '("# PRF input: ",a)') trim(params%profiles%filenm)
+      write(line, '("# PRF input: ",a)') trim(params%profiles%filenm)
-      write(file, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
+      call append(header, line)
+      write(line, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
         params%profiles%iprof, params%profiles%iscal, &
         params%profiles%factne, params%profiles%factte
+      call append(header, line)
       if (params%profiles%iprof > 0) then
-        write(file, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
+        write(line, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
           params%profiles%irho, params%profiles%psnbnd, params%profiles%sspld
+        call append(header, line)
       else
-        write(file, '("# PRF irho psnbnd sspld: N/A (analytical)")')
+        write(line, '("# PRF irho psnbnd sspld: N/A (analytical)")')
+        call append(header, line)
       end if
       ! TODO: add missing values
-      write(file, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0
+      write(line, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0
+      call append(header, line)
     else
-      write(file, '("# PRF input: N/A (vacuum)")')
+      write(line, '("# PRF input: N/A (vacuum)")')
-      write(file, '("# PRF iprof iscal factne factte: N/A (vacuum)")')
+      call append(header, line)
-      write(file, '("# PRF irho psnbnd sspld: N/A (vacuum)")')
+      write(line, '("# PRF iprof iscal factne factte: N/A (vacuum)")')
-      write(file, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
+      call append(header, line)
+      write(line, '("# PRF irho psnbnd sspld: N/A (vacuum)")')
+      call append(header, line)
+      write(line, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
+      call append(header, line)
     end if
 
     ! Antenna parameters
-    write(file, '("# ANT input: ",a)') trim(params%antenna%filenm)
+    write(line, '("# ANT input: ",a)') trim(params%antenna%filenm)
-    write(file, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
+    call append(header, line)
+    write(line, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
       params%antenna%ibeam, params%antenna%iox, &
       params%antenna%psi, params%antenna%chi
-    write(file, '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
+    call append(header, line)
+    write(line, '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
       params%antenna%alpha, params%antenna%beta, &
       params%antenna%power, params%antenna%fghz
-    write(file, '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
+    call append(header, line)
-    write(file, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
+    write(line, '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
+    call append(header, line)
+    write(line, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
       params%antenna%w, params%antenna%ri, params%antenna%phi
+    call append(header, line)
 
     ! Other parameters
-    write(file, '("# RFL rwall:",1x,e12.5)') params%misc%rwall
+    write(line, '("# RFL rwall:",1x,e12.5)') params%misc%rwall
+    call append(header, line)
 
     ! code parameters
-    write(file, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
+    write(line, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
       params%raytracing%igrad, params%raytracing%idst, &
       params%raytracing%ipass, params%raytracing%ipol
-    write(file, '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
+    call append(header, line)
+    write(line, '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
       params%raytracing%nrayr, params%raytracing%nrayth, &
       params%raytracing%nstep, params%raytracing%rwmax, &
       params%raytracing%dst
-    write(file, '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
+    call append(header, line)
+    write(line, '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
       params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, &
       params%ecrh_cd%imx, params%ecrh_cd%ieccd
-    write(file, '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
+    call append(header, line)
+    write(line, '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
       params%output%ipec, params%output%nrho, &
       params%output%istpr, params%output%istpl
-  end subroutine write_gray_header
+    header = trim(header) // trim(line)
+
+  contains
+
+    subroutine append(header, line)
+      ! Appends a line to the header
+      character(len=*), intent(inout) :: header
+      character(len=*), intent(in)    :: line
+      header = trim(header) // trim(line) // new_line('a')
+    end subroutine
+
+  end function make_gray_header
 
 
   function update_parameter(params, name, value) result(err)

src/main.f90

@@ -3,13 +3,16 @@ program main
   use logger,               only : INFO, ERROR, set_log_level, log_message
   use utils,                only : dirname
   use gray_cli,             only : cli_options, parse_cli_options, &
-                                   deinit_cli_options, parse_param_overrides
+                                   parse_param_overrides
   use gray_core,            only : gray_main
   use gray_params,          only : gray_parameters, gray_data, gray_results, &
                                    read_gray_params, read_gray_config,       &
+                                   make_gray_header,                         &
                                    params_set_globals => set_globals
   implicit none
 
+  no_save: block
+
   ! CLI options
   type(cli_options) :: opts
 
@@ -75,25 +78,17 @@ program main
 
   call init_misc(params, data)
 
-  ! Activate the given output tables
-  params%misc%active_tables = opts%tables
-
   ! Get back to the initial directory
   err = chdir(cwd)
 
-  ! Change the current directory to output files there
-  if (allocated(opts%output_dir)) then
-    if (chdir(opts%output_dir) /= 0) then
-      call log_message(level=ERROR, mod='main', &
-        msg='chdir to output_dir ('//opts%output_dir//') failed!')
-      call exit(1)
-    end if
-  end if
-
   if (allocated(opts%sum_filelist)) then
     call log_message(level=INFO, mod='main', msg='summing profiles')
-    call sum_profiles(params, results)
+    call sum_profiles(params, opts%sum_filelist, opts%output_dir, results)
   else
+    ! Activate the given output tables
+    params%misc%active_tables = opts%tables
+
+    ! Finally run the simulation
     call gray_main(params, data, results, err)
   end if
 
@@ -105,15 +100,23 @@ program main
     call log_message(msg, level=INFO, mod='main')
   end block print_res
 
+  ! Change the current directory to output files there
+  if (allocated(opts%output_dir)) then
+    if (chdir(opts%output_dir) /= 0) then
+      call log_message(level=ERROR, mod='main', &
+        msg='chdir to output_dir ('//opts%output_dir//') failed!')
+      call exit(1)
+    end if
+  end if
+
+  ! Write results to file
   write_res: block
     integer :: i
-
-    ! Write results to file
     do i = 1, size(results%tables%iter)
       if (.not. associated(results%tables%iter(i)%ptr)) cycle
 
       associate (tbl => results%tables%iter(i)%ptr)
-      call tbl%save(err, make_header=make_header)
+      call tbl%save(err, header=make_gray_header(params))
       if (err /= 0) then
         call log_message('failed to save '//trim(tbl%title)//' table', &
                          level=ERROR, mod='main')
@@ -123,22 +126,20 @@ program main
   end block write_res
 
   ! Free memory
-  call deinit_equilibrium(data%equilibrium)
+  free_mem: block
-  call deinit_profiles(data%profiles)
+    use coreprofiles, only : unset_profiles_spline
-  call deinit_misc
+    use equilibrium,  only : unset_equil_spline
-  call deinit_cli_options(opts)
+    use limiter,      only : limiter_unset_globals=>unset_globals
-  deallocate(results%dpdv, results%jcd)
+
+    call unset_profiles_spline
+    call unset_equil_spline
+    call limiter_unset_globals
+  end block free_mem
+
+  end block no_save
 
 contains
 
-  subroutine make_header(file)
-    ! Writes the gray output files header
-    use gray_params, only : write_gray_header
-    integer, intent(in) :: file
-    call write_gray_header(params, file)
-  end subroutine
-
-
   subroutine init_equilibrium(params, data, err)
     ! Reads the MHD equilibrium file (either in the G-EQDSK format
     ! or an analytical description) and initialises the respective
@@ -191,25 +192,6 @@ contains
   end subroutine init_equilibrium
 
 
-  subroutine deinit_equilibrium(data)
-    ! Free all memory allocated by the init_equilibrium subroutine.
-    use gray_params, only : equilibrium_data
-    use equilibrium, only : unset_equil_spline
-
-    ! subroutine arguments
-    type(equilibrium_data), intent(inout) :: data
-
-    ! Free the MHD equilibrium arrays
-    if (allocated(data%rv))   deallocate(data%rv, data%zv, data%fpol, data%qpsi)
-    if (allocated(data%psin)) deallocate(data%psin, data%psinr)
-    if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd)
-    if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
-
-    ! Unset global variables of the `equilibrium` module
-    call unset_equil_spline
-  end subroutine deinit_equilibrium
-
-
   subroutine init_profiles(params, factb, launch_pos, data, err)
     ! Reads the plasma kinetic profiles file (containing the elecron
     ! temperature, density and plasma effective charge) and initialises
@@ -272,23 +254,6 @@ contains
   end subroutine init_profiles
 
 
-  subroutine deinit_profiles(data)
-    ! Free all memory allocated by the init_profiles subroutine.
-    use gray_params,  only : profiles_data
-    use coreprofiles, only : unset_profiles_spline
-
-    ! subroutine arguments
-    type(profiles_data), intent(inout) :: data
-
-    ! Free the plasma kinetic profiles arrays
-    if (allocated(data%psrad)) deallocate(data%psrad)
-    if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
-
-    ! Unset global variables of the `coreprofiles` module
-    call unset_profiles_spline
-  end subroutine deinit_profiles
-
-
   subroutine init_antenna(params, err)
     ! Reads the wave launcher file (containing the wave frequency, launcher
     ! position, direction and beam description) and initialises the respective
@@ -390,40 +355,33 @@ contains
   end subroutine init_misc
 
 
-  subroutine deinit_misc
+  subroutine sum_profiles(params, filelist, output_dir, results)
-    ! Free all memory allocated by the init_misc subroutine.
+    ! Combines the EC profiles obtained from multiple gray runs
-    use limiter, only : limiter_unset_globals=>unset_globals
+    ! (of different beams) and recomputes the summary table
-
-    ! Unset the global variables of the `limiter` module
-    call limiter_unset_globals
-  end subroutine deinit_misc
-
-
-  subroutine sum_profiles(params, results)
     use gray_params,  only : gray_parameters
     use magsurf_data, only : compute_flux_averages, dealloc_surfvec
-    use beamdata,     only : init_btr
     use pec,          only : pec_init, postproc_profiles, dealloc_pec, &
                              rhot_tab
 
     ! subroutine arguments
     type(gray_parameters), intent(inout) :: params
+    character(len=*),      intent(in)    :: filelist, output_dir
     type(gray_results),    intent(inout) :: results
 
+    ! local variables
     real(wp_), dimension(:), allocatable :: jphi
     real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
     real(wp_), dimension(:), allocatable :: extracol
-    integer :: i, j, k, l, n, ngam, nextracol, irt
+    integer :: i, j, k, l, nbeams, ncases, nextracols, irt, err
-    character(len=255) :: filename, headerline, fmtstr
+    character(len=255) :: filename, line, extralabels, fmtstr
-    real(wp_), dimension(5) :: f48v
+    real(wp_), dimension(5) :: sums
     real(wp_) :: jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav,  &
                  drhotjava, drhotpav, ratjamx, ratjbmx
-
-    ! local variables
     real(wp_) :: drhotp, drhotj, dpdvmx, jphimx
 
-    ! Initialise the ray variables (beamtracing)
+    ! unit numbers
-    call init_btr(params%raytracing)
+    integer :: list, prof, summ
+    integer, allocatable :: beams(:)
 
     ! Initialise the magsurf_data module
     call compute_flux_averages
@@ -436,80 +394,113 @@ contains
       allocate(results%jcd(nrho), results%dpdv(nrho))
     end associate
 
-    open(100, file=opts%sum_filelist, action='read', status='old', iostat=err)
+    ! Read the list of files to be summed
+    open(newunit=list, file=filelist, action='read', status='old', iostat=err)
     if (err /= 0) then
       call log_message(level=ERROR, mod='main', &
-        msg='opening file list ('//opts%sum_filelist//') failed!')
+        msg='opening file list ('//filelist//') failed!')
       call exit(1)
     end if
 
-    open(unit=100 -1, file='f48sum.txt', action='write', status='unknown')
+    ! Open combined output files
-    open(unit=100 -2, file='f7sum.txt', action='write', status='unknown')
+    open(newunit=prof, file=output_dir//'/'//'sum-ec-profiles.txt', action='write')
+    open(newunit=summ, file=output_dir//'/'//'sum-summary.txt', action='write')
 
-    read(100, *) n, ngam, nextracol
+    ! Move to the filelist directory
-    allocate(extracol(nextracol))
+    ! (the filepaths are assumed relative to it)
-    do i=1,n
+    if (chdir(dirname(filelist)) /= 0) then
-      read(100, *) filename
+      call log_message(level=ERROR, mod='main', &
-      open(100 + i, file=filename, action='read', status='old', iostat=err)
+        msg='chdir to file list directory failed!')
+      call exit(1)
+    end if
+
+    ! Read the filelist header
+    read(list, *) nbeams, ncases, nextracols
+    allocate(beams(nbeams), extracol(nextracols))
+
+    ! Read the list of files
+    do i = 1, nbeams
+      read(list, *) filename
+
+      ! Open the summand file
+      open(newunit=beams(i), file=filename, &
+           action='read', status='old', iostat=err)
       if (err /= 0) then
         call log_message(level=ERROR, mod='main', &
           msg='opening summand file ('//trim(filename)//') failed!')
         call exit(1)
       end if
 
-      do j=1,22
+      ! Copy the Gray header from the ec-profiles file
-        read(100 + i, '(a255)') headerline
+      do j = 1, 21
-        if (i == 1) then
+        read(beams(i), '(a255)') line
-          write(100 -1, '(a)') trim(headerline)
+        if (i == 1) write(prof, '(a)') trim(line)
-          write(100 -2, '(a)') trim(headerline)
+        if (i == 1) write(summ, '(a)') trim(line)
-        end if
       end do
-    end do
 
-    do k=1,ngam
+      ! Copy the column labels
-      jphi         = 0
+      ! Note: for the summary.7 we copy only the first columns
+      read(beams(i), '(a255)') line
+      if (i == 1) then
+        extralabels = line(1:index(line, "ρ_p")-1)
+        write(prof, '(a)') trim(line)
+        write(summ, '(a,x,a,x,a)') trim(extralabels), &
+          "I_cd P_abs J_φ_peak dPdV_peak ρ_max_J ρ_avg_J ρ_max_P", &
+          "ρ_avg_P Δρ_avg_J Δρ_avg_P ratio_Ja_max ratio_Jb_max"
+      end if
+    end do
+    close(list)
+
+    do k = 1, ncases
+      jphi    = 0
+      currins = 0
+      pins    = 0
       results%jcd  = 0
       results%dpdv = 0
-      currins      = 0
+
-      pins         = 0
+      ! Sum the radial profiles
-      do j=1,params%output%nrho
+      do j = 1, params%output%nrho
-        do i=1,n
+        do i = 1, nbeams
-          read(100+i, *) (extracol(l), l=1,nextracol), &
+          read(beams(i), *) &
-              rpin(j), rtin(j), f48v(1:5), irt
+            (extracol(l), l=1,nextracols), rpin(j), rtin(j), sums(1:5), irt
-          jphi(j)         = f48v(1) + jphi(j)
+          jphi(j)         = sums(1) + jphi(j)
-          results%jcd(j)  = f48v(2) + results%jcd(j)
+          results%jcd(j)  = sums(2) + results%jcd(j)
-          results%dpdv(j) = f48v(3) + results%dpdv(j)
+          results%dpdv(j) = sums(3) + results%dpdv(j)
-          currins(j)      = f48v(4) + currins(j)
+          currins(j)      = sums(4) + currins(j)
-          pins(j)         = f48v(5) + pins(j)
+          pins(j)         = sums(5) + pins(j)
         end do
-        write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7
+        ! Store the sums
-        write(100 -1,trim(fmtstr))                           &
+        write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracols + 7
-            (extracol(l), l=1,nextracol), rpin(j), rtin(j),  &
+        write(prof, trim(fmtstr)) &
-            jphi(j), results%jcd(j), results%dpdv(j),        &
+          (extracol(l), l=1,nextracols), rpin(j), rtin(j), jphi(j), &
-            currins(j), pins(j), irt
+          results%jcd(j), results%dpdv(j), currins(j), pins(j), irt
       end do
+
+      ! Combine the total power, current
       results%pabs = pins(params%output%nrho)
       results%icd = currins(params%output%nrho)
-      write(100 -1, *)
+
-      ! Recompute final results
+      ! Recompute the summary values
-      call postproc_profiles(results%pabs, results%icd, rhot_tab,    &
+      call postproc_profiles( &
-                             results%dpdv, jphi,  rhotpav, drhotpav, &
+        results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi,  &
-                             rhotjava, drhotjava, dpdvp, jphip,      &
+        rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,     &
-                             rhotp, drhotp, rhotj, drhotj,           &
+        rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx)
-                             dpdvmx, jphimx, ratjamx, ratjbmx)
+
-      write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12
+      write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12
-      write(100 -2, trim(fmtstr))                                   &
+      write(summ, trim(fmtstr)) &
-           (extracol(l), l=1,nextracol), results%icd, results%pabs, &
+        (extracol(l), l=1,nextracols), results%icd, results%pabs, &
-           jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav,           &
+        jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav,  drhotjava, &
-           drhotjava, drhotpav, ratjamx, ratjbmx
+        drhotpav, ratjamx, ratjbmx
+
+      ! Separate the case
+      write(prof, *)
     end do
 
-    do i=-2,n
+    ! Close all files
-      close(100 + i)
+    close(prof); close(summ)
+    do i = 1, nbeams
+      close(beams(i))
     end do
-    deallocate(jphi, currins, pins, rtin, rpin)
-    deallocate(extracol)
-    deallocate(opts%params_file)
 
     ! Free memory
     call dealloc_surfvec

src/main_convert.f90

@@ -22,6 +22,8 @@ program gray_convert
     integer(kind(config_format))  :: format       ! output file format
   end type
 
+  no_save: block
+
   type(cli_options) :: opts
   type(gray_parameters) :: params
   integer :: err
@@ -47,6 +49,8 @@ program gray_convert
 
   call exit(0)
 
+  end block no_save
+
 contains
 
   subroutine print_help()

src/types.f90

@@ -50,13 +50,6 @@ module types
      procedure :: wrap_init
   end interface
 
-  ! Interface for subroutines that write to file
-  interface
-    subroutine writer(file)
-      integer, intent(in) :: file
-    end subroutine
-  end interface
-
 contains
 
   subroutine queue_put(self, val)
@@ -175,17 +168,17 @@ contains
   end subroutine table_append
 
 
-  subroutine table_save(self, error, filepath, make_header)
+  subroutine table_save(self, error, filepath, header)
     ! Save the table to a file
     !
     ! Note: this operation consumes the table
     use utils, only : get_free_unit, digits
 
     ! suborutine arguments
-    class(table),      intent(inout)         :: self
+    class(table), intent(inout)         :: self
-    integer,           intent(out)           :: error
+    integer,      intent(out)           :: error
-    character(*),      intent(in), optional  :: filepath
+    character(*), intent(in), optional  :: filepath
-    procedure(writer),             optional  :: make_header
+    character(*),             optional  :: header
 
     ! local variables
     integer :: i, file
@@ -209,7 +202,7 @@ contains
     if (error /= 0) return
 
     ! write a custom file header
-    if (present(make_header)) call make_header(file)
+    if (present(header)) write(file, '(a)') trim(header)
 
     ! write table header
     write(file, '(a)', advance='no') '#'