From ba8fc001e55452e24413c7a1d9573cfb5b93f8e1 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 16 May 2024 09:19:22 +0200 Subject: [PATCH] fix memory errors and use automatic deallocation --- src/gray_cli.f90 | 17 +-- src/gray_params.f90 | 116 +++++++++++++------ src/main.f90 | 261 +++++++++++++++++++++---------------------- src/main_convert.f90 | 4 + src/types.f90 | 19 +--- 5 files changed, 216 insertions(+), 201 deletions(-) diff --git a/src/gray_cli.f90 b/src/gray_cli.f90 index 497d1c0..0d45c1c 100644 --- a/src/gray_cli.f90 +++ b/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, & - print_version, get_next_command + parse_param_overrides, 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 diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 6ffa6dd..d96800e 100644 --- a/src/gray_params.f90 +++ b/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=10) :: time + character(len=200) :: line + 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(file, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') + write(line, '("# EQL input: 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(file, '("# EQL ssplps ssplf ixp: N/A (analytical)")') + write(line, '("# EQL icocos ipsinorm idesc ifreefmt: 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(file, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') & + write(line, '("# PRF input: ",a)') trim(params%profiles%filenm) + 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(file, '("# PRF iprof iscal factne factte: N/A (vacuum)")') - write(file, '("# PRF irho psnbnd sspld: N/A (vacuum)")') - write(file, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') + write(line, '("# PRF input: N/A (vacuum)")') + call append(header, line) + write(line, '("# PRF iprof iscal factne factte: 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(file, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') & + write(line, '("# ANT input: ",a)') trim(params%antenna%filenm) + 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 - write(file, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') & + call append(header, line) + 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) diff --git a/src/main.f90 b/src/main.f90 index 5dd4e2f..bc1485b 100644 --- a/src/main.f90 +++ b/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) - call deinit_profiles(data%profiles) - call deinit_misc - call deinit_cli_options(opts) - deallocate(results%dpdv, results%jcd) + free_mem: block + use coreprofiles, only : unset_profiles_spline + use equilibrium, only : unset_equil_spline + use limiter, only : limiter_unset_globals=>unset_globals + + 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 - ! Free all memory allocated by the init_misc subroutine. - use limiter, only : limiter_unset_globals=>unset_globals - - ! Unset the global variables of the `limiter` module - call limiter_unset_globals - end subroutine deinit_misc - - - subroutine sum_profiles(params, results) + subroutine sum_profiles(params, filelist, output_dir, results) + ! Combines the EC profiles obtained from multiple gray runs + ! (of different beams) and recomputes the summary table 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 - character(len=255) :: filename, headerline, fmtstr - real(wp_), dimension(5) :: f48v + integer :: i, j, k, l, nbeams, ncases, nextracols, irt, err + character(len=255) :: filename, line, extralabels, fmtstr + 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) - call init_btr(params%raytracing) + ! unit numbers + 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(unit=100 -2, file='f7sum.txt', action='write', status='unknown') + ! Open combined output files + 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 - allocate(extracol(nextracol)) - do i=1,n - read(100, *) filename - open(100 + i, file=filename, action='read', status='old', iostat=err) + ! Move to the filelist directory + ! (the filepaths are assumed relative to it) + if (chdir(dirname(filelist)) /= 0) then + call log_message(level=ERROR, mod='main', & + 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 - read(100 + i, '(a255)') headerline - if (i == 1) then - write(100 -1, '(a)') trim(headerline) - write(100 -2, '(a)') trim(headerline) - end if + ! Copy the Gray header from the ec-profiles file + do j = 1, 21 + read(beams(i), '(a255)') line + if (i == 1) write(prof, '(a)') trim(line) + if (i == 1) write(summ, '(a)') trim(line) end do - end do - do k=1,ngam - jphi = 0 + ! Copy the column labels + ! 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 - do j=1,params%output%nrho - do i=1,n - read(100+i, *) (extracol(l), l=1,nextracol), & - rpin(j), rtin(j), f48v(1:5), irt - jphi(j) = f48v(1) + jphi(j) - results%jcd(j) = f48v(2) + results%jcd(j) - results%dpdv(j) = f48v(3) + results%dpdv(j) - currins(j) = f48v(4) + currins(j) - pins(j) = f48v(5) + pins(j) + + ! Sum the radial profiles + do j = 1, params%output%nrho + do i = 1, nbeams + read(beams(i), *) & + (extracol(l), l=1,nextracols), rpin(j), rtin(j), sums(1:5), irt + jphi(j) = sums(1) + jphi(j) + results%jcd(j) = sums(2) + results%jcd(j) + results%dpdv(j) = sums(3) + results%dpdv(j) + currins(j) = sums(4) + currins(j) + pins(j) = sums(5) + pins(j) end do - write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7 - write(100 -1,trim(fmtstr)) & - (extracol(l), l=1,nextracol), rpin(j), rtin(j), & - jphi(j), results%jcd(j), results%dpdv(j), & - currins(j), pins(j), irt + ! Store the sums + write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracols + 7 + write(prof, trim(fmtstr)) & + (extracol(l), l=1,nextracols), rpin(j), rtin(j), jphi(j), & + 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 - call postproc_profiles(results%pabs, results%icd, rhot_tab, & - results%dpdv, jphi, rhotpav, drhotpav, & - rhotjava, drhotjava, dpdvp, jphip, & - rhotp, drhotp, rhotj, drhotj, & - dpdvmx, jphimx, ratjamx, ratjbmx) - write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12 - write(100 -2, trim(fmtstr)) & - (extracol(l), l=1,nextracol), results%icd, results%pabs, & - jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & - drhotjava, drhotpav, ratjamx, ratjbmx + + ! Recompute the summary values + call postproc_profiles( & + results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi, & + rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & + rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) + + write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12 + write(summ, trim(fmtstr)) & + (extracol(l), l=1,nextracols), results%icd, results%pabs, & + jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, drhotjava, & + drhotpav, ratjamx, ratjbmx + + ! Separate the case + write(prof, *) end do - do i=-2,n - close(100 + i) + ! Close all files + 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 diff --git a/src/main_convert.f90 b/src/main_convert.f90 index aa6f95a..a734a15 100644 --- a/src/main_convert.f90 +++ b/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() diff --git a/src/types.f90 b/src/types.f90 index ea53fc9..2f8c153 100644 --- a/src/types.f90 +++ b/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 - integer, intent(out) :: error - character(*), intent(in), optional :: filepath - procedure(writer), optional :: make_header + class(table), intent(inout) :: self + integer, intent(out) :: error + character(*), intent(in), optional :: filepath + 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') '#'