From e6d864e63b20ff2d255fcf4d726e86d053f193f2 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Thu, 24 Nov 2022 23:37:29 +0100
Subject: [PATCH] Fix the summation mode + temporary fix for file unit numbers

- Add missing array allocations
- Add parameter for varying number of columns in input files
- Change output unit numbers (dirty fix. Original units created an empty
  named file, but wrote in default named fort.* files)
---
 src/main.f90 | 84 +++++++++++++++++++++++++++++-----------------------
 1 file changed, 47 insertions(+), 37 deletions(-)

diff --git a/src/main.f90 b/src/main.f90
index 275b586..5ba6d69 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -81,13 +81,13 @@ program main
     call log_message(level=INFO, mod='main', msg='summing profiles')
 
     sum: block
-      real(wp_) :: pabs, icd, pec
-      real(wp_), dimension(:), allocatable :: dpdv, jcd, jphi
+      real(wp_), dimension(:), allocatable :: jphi
       real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
-      integer :: i, j, k, n, ngam, irt
-      character(len=255) :: filename
+      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
-      real(wp_) :: gam,alp,bet, jphip,dpdvp, &
+      real(wp_) :: jphip,dpdvp, &
                    rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
       allocate(jphi(params%output%nrho), currins(params%output%nrho), &
                pins(params%output%nrho), rtin(params%output%nrho),    &
@@ -100,7 +100,11 @@ program main
         call exit(1)
       end if
 
-      read(100, *) n, ngam
+      open(unit=100 -1, file='f48sum.txt', action='write', status='unknown')
+      open(unit=100 -2, file='f7sum.txt', action='write', status='unknown')
+
+      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)
@@ -111,49 +115,55 @@ program main
         end if
 
         do j=1,22
-          read(100 + i, *)
+          read(100 + i, '(a255)') headerline
+          if (i == 1) then
+            write(100 -1, '(a)') trim(headerline)
+            write(100 -2, '(a)') trim(headerline)
+          end if
         end do
       end do
-      close(100)
-
-      open(100 + n+1, file='f48sum.txt', action='write', status='unknown')
-      open(100 + n+2, file='f7sum.txt', action='write', status='unknown')
 
+      allocate(results%jcd(params%output%nrho), results%dpdv(params%output%nrho))
       do k=1,ngam
-        jphi    = zero
-        jcd     = zero
-        dpdv    = zero
-        currins = zero
-        pins    = zero
+        jphi         = zero
+        results%jcd  = zero
+        results%dpdv = zero
+        currins      = zero
+        pins         = zero
         do j=1,params%output%nrho
           do i=1,n
-            read(100+i, *) gam, alp, bet, rpin(j), rtin(j), f48v(1:5), irt
-            jphi(j)    = f48v(1) + jphi(j)
-            jcd(j)     = f48v(2) + jcd(j)
-            dpdv(j)    = f48v(3) + dpdv(j)
-            currins(j) = f48v(4) + currins(j)
-            pins(j)    = f48v(5) + pins(j)
+            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)
           end do
-          write(100 + n+1,'(10(1x,e16.8e3),i5)')                  &
-              gam, alp, bet, rpin(j), rtin(j),                    &
-              jphi(j), jcd(j), dpdv(j), currins(j), pins(j), irt
+          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
         end do
-        pec = pins(params%output%nrho)
-        icd = currins(params%output%nrho)
-        write(100 + n+1, *)
-        call sum_profiles(params, jphi, jcd, dpdv, currins,     &
-                          pins, pabs, icd, jphip, dpdvp, rhotj, &
-                          rhotjava, rhotp, rhotpav, drhotjava,  &
-                          drhotpav, ratjamx, ratjbmx)
-        write(100 + n+2, '(15(1x,e12.5),i5,4(1x,e12.5))') &
-             gam, alp, bet, icd, pabs, jphip, dpdvp,      &
-             rhotj, rhotjava, rhotp, rhotpav,             &
+        results%pabs = pins(params%output%nrho)
+        results%icd = currins(params%output%nrho)
+        write(100 -1, *)
+        call sum_profiles(params, jphi, results%jcd, results%dpdv,       &
+                          currins, pins, results%pabs, results%icd,      &
+                          jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
+                          drhotjava, drhotpav, 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
       end do
-      do i=1,n+2
+      do i=-2,n
         close(100 + i)
       end do
-     deallocate(dpdv, jcd, jphi, currins, pins, rtin, rpin)
+     deallocate(jphi, currins, pins, rtin, rpin)
+     deallocate(extracol)
      deallocate(opts%params_file)
    end block sum
   else