PROGRAM GATHER
    
    !USE PARAMETERS
    !USE GLOBAL 
    
    IMPLICIT NONE  
    
    !! ARGUMENTS FROM COMMAND LINE
    CHARACTER (LEN = 80)    :: cmdline      ! generic command line  
    INTEGER                 :: argtot       ! argument count
    INTEGER                 :: Np = 720
    
    INTEGER, PARAMETER      :: N_k = 120, N_z = 7, jmax = 3, ncaptdist = 521
    REAL*8, PARAMETER :: onedp = 1.0D+00, zerodp = 0.0D+00
    
    REAL*8    :: repnum, hoscap, inf_init,                        &
                    beta_adj, beta_e, beta_c, beta_l, pix,          &
                    sui_cbar, hcut, mitdur, hloc, gam, sui_curv
    REAL*8    :: wss, rss, tauc, taus, bigs, bigi, bigr,              &
                    aggkinit, aggcinit, aggyinit, agglinit, agghinit
    
    REAL*8, DIMENSION(ncaptdist)  :: wpath, rpath,                    &
                    aggkpath, aggcpath, aggypath, agglpath, agghpath,   &
                    bigspath, bigipath, bigrpath,                       &
                    dpath, aggtaxpath, aggsuipath, poppath
    
    REAL*8    :: pvdef, deaths, avgwelf, poswelf, peak_inf, yloss
                        
    INTEGER     :: days, lockdown
    
    REAL*8, DIMENSION(jmax,N_z,N_k)       :: muss, welf, welf_nomit
    
    REAL*8      :: dum
    INTEGER     :: i, j, dumi, idx_k, idx_z, t
    
    CHARACTER(LEN=12) string, FORMAT
    LOGICAL     ::  file_exists
    CHARACTER(LEN=19)   :: path
    
     
102 FORMAT(4I16, 12F16.10, I16, 3F16.10)  
499 FORMAT(f18.10,7f18.10,7f18.10,7f18.10)
501 FORMAT (I5, 67f15.10)
    
    OPEN (UNIT=15, FILE='results.csv', POSITION="REWIND")
    WRITE (15,'(*(G0.10,:,","))') 'subsidy','hours threshold','hours cap', &
                'duration','lockdown','Avg welf', &
                'peak inf','deaths','pos welf', &
                'pv deficit','tauc','taus','output loss','config','sui_curv'
   
    ! Get command count
    argtot      = COMMAND_ARGUMENT_COUNT()
    IF (argtot > 0) THEN
        CALL GET_COMMAND_ARGUMENT(NUMBER = 1, VALUE = cmdline)
        READ (cmdline,*) Np
    ENDIF
    
    ! Loop across input pairs and store them    
    DO i=1,Np   
    
		DATA FORMAT /'(I6.6)'/
		WRITE(string,FORMAT) i
	    path = 'output/config'//TRIM(string)

        
        INQUIRE(file=TRIM(path)//'/PARAMS_MATLAB.out',exist=file_exists)
        IF (file_exists) THEN
            
            OPEN (11, FILE=TRIM(path)//'/PARAMS_MATLAB.out', POSITION="REWIND")
            READ(11,*)
            READ (11,102)  days, dumi, dumi, dumi,                          &
                            beta_adj, beta_e, beta_c, beta_l, pix,          &
                            hoscap, dum, dum, dum,                          &
                            sui_cbar, hcut, mitdur, lockdown, hloc,			&
							gam, sui_curv
            CLOSE (11)   
            
            OPEN(5,file=TRIM(path)//'/muss.out', POSITION="REWIND")
            DO idx_k = 1,N_k
                READ (5,499) dum,muss(1,:,idx_k),muss(2,:,idx_k),muss(3,:,idx_k)
            END DO
            CLOSE(5)

            OPEN (UNIT=13, FILE=TRIM(path)//'/equilibrium_tran.out', POSITION="REWIND")
            READ(13,*)
            READ (13,501)  dumi, wss, rss,                                      &
                            aggkinit, aggcinit, aggyinit, agglinit, agghinit,   &
                            bigs, bigi, bigr,                                   &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum, tauc, taus, pvdef,                             &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum,dum,                                        &
                            dum,dum
							
            DO t=1,ncaptdist
                READ (13,501)  dumi, wpath(t), rpath(t),                                        &
							aggkpath(t), aggcpath(t), aggypath(t), agglpath(t), agghpath(t),    &
							bigspath(t), bigipath(t), bigrpath(t),                              &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dpath(t), aggtaxpath(t), aggsuipath(t), dum,                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                                        &
							dum,dum,dum,                                        				&
                            dum
            END DO ! t=1,ncaptdist
            CLOSE (UNIT=13)
            
            IF (aggtaxpath(2)==zerodp) THEN
                tauc = zerodp
                taus = zerodp
            END IF

            poppath = bigipath+bigspath+bigrpath
            deaths = SUM(dpath)
            
			peak_inf = MAXVAL(bigipath)
			yloss = onedp-SUM(aggypath(1:52))/(aggyinit*52.0D+00)
            
            INQUIRE(file=TRIM(path)//'/welf.out',exist=file_exists)
            IF (file_exists) THEN
                OPEN(7,file=TRIM(path)//'/welf.out', POSITION="REWIND")
                DO idx_k = 1,N_k
                    READ(7,499) dum,welf(1,:,idx_k),welf(2,:,idx_k),welf(3,:,idx_k)
                END DO
                CLOSE(7) 
            
                avgwelf = 100.0D+00*SUM(muss*welf)/SUM(muss)
                IF (i>1) THEN
                    poswelf = zerodp
                    DO j=1,jmax
                    DO idx_z=1,N_z
                    DO idx_k=1,N_k
                        IF (welf(j,idx_z,idx_k).GE.welf_nomit(j,idx_z,idx_k)) THEN
                            poswelf = poswelf + muss(j,idx_z,idx_k)
                        END IF
                    END DO
                    END DO
                    END DO
                ELSE
                    poswelf = zerodp
                    welf_nomit = welf
                END IF
                poswelf = poswelf*100.0D+00

                WRITE (15,'(*(G0.10,:,","))') sui_cbar,hcut,hloc, &
                    mitdur,lockdown,avgwelf, &
                    peak_inf,deaths,poswelf, &
                    pvdef,tauc,taus,yloss,i,sui_curv
			ELSE
				WRITE (15,'(*(G0.10,:,","))') sui_cbar,hcut,hloc, &
                    mitdur,lockdown,-1d+2, &
                    peak_inf,deaths,-1d+2, &
                    pvdef,tauc,taus,yloss,i,sui_curv	
            END IF
    
        END IF
    END DO
    CLOSE (UNIT=15)  
    
    
    
END PROGRAM GATHER