! # given a WEPS run directory, parse the dec#.btmp files into unique crop files

program agedecompbycrop

! for command line processing
integer :: iargc
character*256 :: rundirectory, argv

! others
integer :: maxfile ! maximum number of residue age pool input files that can be used
integer :: maxname ! maximum number of unique resiude name that can be used
parameter ( maxfile = 10 ) ! don't set greater than 99
parameter ( maxname = 20 ) ! don't set greater than 99

integer :: idx, jdx ! loop counters
character*15 :: decfile(maxfile) ! array of input file names
integer :: tfile ! maximum number of input files found
integer :: ifile ! index set to the residue being selected from the file
integer :: lundec(maxfile) ! array of logical unit number used for input files
integer :: ios(maxname) ! input output status return for each file
character*512 :: inline(maxfile) ! input line read for each file
integer :: tname ! total number of unique crop names
integer :: iname ! index set to the residue being selected from the file
integer :: nstart, nend ! start and end character count for name field
character*30 :: linename(maxfile) ! residue name for present line in each residue file
character*30 :: uniquename(maxname) ! to hold unique residue names
character*15 :: resfile(maxname) ! array of output file names
integer :: lunres(maxname) ! logical unit number for residue output file
character*15 :: restatfile(maxname) ! array of output file names
integer :: lunrestat(maxname) ! logical unit number for output files containing min, max, mean, sd of each resday for each unique residue
integer :: iosres ! input output status return for output file
integer :: skip ! flag to indicate name is not unique

integer :: daysim(maxfile), resday(maxfile), resyear(maxfile), doy, year, pool
real :: cumddysta(2), cumddyflat(2), cumddybg10(2)
integer :: prevday ! check that residue days are continuous
integer :: prevyear ! residue year check

integer :: lundecall ! logical unit number for output file containing all pools in residue number order
integer :: iresyr ! residue year counter
integer :: iosall ! input output status for decall.btmp file

integer :: idata(6)
real :: rdata(30)

integer :: maxresday(maxname)
integer, allocatable :: cntresyr(:) ! count number of values for each residue day (n value for across year statistics)
real, allocatable :: resmin(:,:), resmax(:,:), resmean(:,:), resvar(:,:), ressd(:,:), resmnmsd(:,:), resmnpsd(:,:)
integer :: oldmincnt(30), oldmaxcnt(30), newmincnt(30), newmaxcnt(30) ! counts of exceedences on residue curve
real :: delta

! set character fields for name parsing
nstart = 361
nend = 390

! create age pool input file names
! assign first three characters to file name
decfile(1:maxfile) (1:3) = 'dec'
! assign last four characters to file name
decfile(1:maxfile) (6:10) = '.btmp'
do idx = 1, maxfile
    ! assign 4 character of file name to be the number character corresponsing to tens place of idx
    decfile(idx) (4:4) = char(48+(idx/10))
    ! assign 5 character of file name to be the number character corresponsing to ones place of idx
    decfile(idx) (5:5) = char(48+(idx-(idx/10)*10))
    ! display created file
    !write(*,*) 'File name created: ', decfile(idx)
    ! assign logical unit number of opening file to array
    lundec(idx) = 100+idx
end do

!single output file
lundecall = 100

! create residue output file names
! assign first three characters to file name
resfile(1:maxname) (1:3) = 'res'
restatfile(1:maxname) (1:3) = 'rst'
! assign last four characters to file name
resfile(1:maxname) (6:10) = '.btmp'
restatfile(1:maxname) (6:10) = '.btmp'
do idx = 1, maxname
    ! assign 4 character of file name to be the number character corresponsing to tens place of idx
    resfile(idx) (4:4) = char(48+(idx/10))
    restatfile(idx) (4:4) = char(48+(idx/10))
    ! assign 5 character of file name to be the number character corresponsing to ones place of idx
    resfile(idx) (5:5) = char(48+(idx-(idx/10)*10))
    restatfile(idx) (5:5) = char(48+(idx-(idx/10)*10))
    ! display created file
    !write(*,*) 'File name created: ', resfile(idx)
    ! assign logical unit number of opening file to array
    lunres(idx) = 200+idx
    lunrestat(idx) = 250+idx
end do

! process command line
if (iargc() .eq. 1) then
    call getarg(1,argv)
    rundirectory = argv
else
    write(*,*) 'provide a run directory (ending in .wjr) after the command'
    stop
endif

! check the directory for the the decomposition age files
do idx = 1, maxfile
    open(lundec(idx), file=rundirectory(1:len_trim(rundirectory)) // &
         decfile(idx)(1:len_trim(decfile(idx))), status="old", iostat=ios(idx))
    if( ios(idx) .ne. 0 ) then
        ! write(*,*) 'ERROR, unable to open output file named: ', decfile(idx)
        tfile = idx-1
        exit
    end if
end do

write(*,*) 'Found ', tfile, 'residue age output files'

! create single file ordered by resyear (contains all crops, skips resyear 0 which is "No Crop")
! open single output file for residue
open(lundecall, file=rundirectory(1:len_trim(rundirectory)) // 'decall.btmp', status="unknown", iostat=iosres)

! get all input files started
do ifile = 1, tfile
    ! read first line (header) and discard
    read(lundec(ifile),'(a)',iostat=ios(ifile)) inline(ifile)
    if( ifile .eq. 1 ) then
        ! write header to output file
        write(lundecall,'(a)') inline(ifile)(1:len_trim(inline(ifile)))
    end if
    ! read next line
    read(lundec(ifile),'(a)',iostat=ios(ifile)) inline(ifile)
end do

ifile = 1
iresyr = 1
do while( ios(ifile) .eq. 0 )
    ! parse daysim, resday and resyear
    read(inline(ifile),*) daysim(ifile), resday(ifile), resyear(ifile)
    if( iresyr .eq. resyear(ifile) ) then 
        ! correct year, write line to file
        write(lundecall,'(a)') inline(ifile)(1:nend)
        ! read next line
        read(lundec(ifile),'(a)',iostat=ios(ifile)) inline(ifile)
    else if( iresyr .gt. resyear(ifile) ) then
        ! file needs to catch up, keep reading lines
        read(lundec(ifile),'(a)',iostat=ios(ifile)) inline(ifile)
    else ! if( iresyr .lt. resyear(ifile) ) then
        ! file is ahead, move to next file
        ifile = ifile + 1
        if( ifile .gt. tfile ) then
            ! we just finished the last file, this resyear is complete, jump back to first file
            ifile = 1
            ! insert blank lines into output for graphing
            write(lundecall,*)
            write(lundecall,*)
            ! increment resyear
            iresyr = iresyr + 1
            ! show year
            ! write(*,*) 'processing residue year ', iresyr
        end if
    end if

    ! check for this file at end of file, keep reading until all files complete
    do while( ios(ifile) .ne. 0 )
        ! this file is at end, go to next file that is not
        if( ifile .eq. tfile ) then
            ! this is last file, go back to 1
            ifile = 1
            ! since last file is done, reduce number of files
            tfile = tfile - 1
        else if( ifile .eq. 1 ) then
            ! first file is finished, we are really done
            exit
        else
            ! keep going
            ifile = ifile + 1
        end if
    end do
end do

close(lundecall)

write(*,*) 'Found and wrote ', iresyr, ' residue years'

!close the residue files (dec??.btmp)
do idx = 1, tfile
    close( lundec(idx) )
end do

! open decall.btmp for reading (this file has all names, others may not)
open(lundecall, file=rundirectory(1:len_trim(rundirectory)) // & 
     "decall.btmp", status="old", position="rewind", iostat=iosall)

! read first file (shorter than decall.btmp) and find number of unique residue types by name
tname=0
skip = 0
! read header line and discard
read(lundecall,'(a)',iostat=ios(1)) inline(1)
! read first data line
read(lundecall,'(a)',iostat=ios(1)) inline(1)
! loop until done reading file
do while( ios(1) .eq. 0 )
    ! residue name substring
    linename(1) = inline(1)(nstart:nend)
    ! check for No Crop and blank line
    if( (linename(1)(1:7) .ne. 'No Crop') .and. (linename(1)(1:7) .ne. '                    ') ) then
        ! check for unique name
        if( tname .eq. 0 ) then
            tname = 1
            uniquename(tname) = linename(1)
            write(*,*) 'Residue Name: ', tname, ' is ', uniquename(tname)
        else
            do idx = 1, tname
                if( linename(1) .eq. uniquename(idx) ) then
                    skip = 1
                    exit
                end if
            end do
            if( skip .eq. 0 ) then
                ! name did not match any other in array
                ! increment number of unique names
                tname = tname + 1
                ! assign new unique name to array
                uniquename(tname) = linename(1)
                ! show unique names found
                write(*,*) 'Residue Name: ', tname, ' is ', uniquename(tname)
            end if
        end if
    end if
    ! reset flag for next check
    skip = 0
    read(lundecall,'(a)',iostat=ios(1)) inline(1)
end do

! reads the file with all residue and breaks it into one file for each unique residue
! this also finds the maximum residue days for each unique residue
! to be used when finding the max, average, and min values

! rewind decall.btmp for reading again
rewind(lundecall)

! read header line
read(lundecall,'(a)',iostat=iosall) inline(1)

! open file for each named residue and place header line and two blank line at beginning of each file
do iname = 1, tname
    ! open output file for residue
    open(lunres(iname), file=rundirectory(1:len_trim(rundirectory)) // &
         resfile(iname)(1:len_trim(resfile(iname))), status="unknown", iostat=ios(iname))
    ! write header line
    write(lunres(iname),'(a)') inline(ifile)(1:nend)

    ! intialize maxresday
    maxresday(iname) = 0
end do

! read next line
read(lundecall,'(a)',iostat=iosall) inline(1)

write(*,*) 'Creating single residue files'

iname = 0 ! no name has been found yet
do while( iosall .eq. 0 )
    
    ! parse residue name
    linename(1) = inline(1)(nstart:nend)

    ! sort lines into each respective file
    if( iname .eq. 0 ) then
        ! no name has been found yet
        ! check against unique names
        do idx = 1, tname
            if( linename(1) .eq. uniquename(idx) ) then
                ! the line matches unique name, set index
                iname = idx
                ! no need to check more, exit loop
                exit
            end if
        end do
        if( iname .eq. 0 ) then
            ! no name matches, read another line
            read(lundecall,'(a)',iostat=iosall) inline(1)
        end if
    else
        ! a named residue block is already being output
        if( linename(1) .eq. uniquename(iname) ) then
            ! this name matches, so wrote line to numbered output file
            write(lunres(iname),'(a)') inline(ifile)(1:nend)
            ! this is a line with a name, parse daysim and resday
            read( inline(1),*) daysim(1), resday(1)
            ! check for maximum resday
            maxresday(iname) = max(maxresday(iname), resday(1))
            ! read another line
            read(lundecall,'(a)',iostat=iosall) inline(1)
        else
            ! the residue name does not match, could be a new name or blank line
            ! regardless, this is end of block so write two blank line for graphing
            write(lunres(iname),*)
            write(lunres(iname),*)

            !write(*,*) 'Wrote block ',maxresday(iname),' of ',uniquename(iname),' to ',resfile(iname)

            ! no residue name is selected
            iname = 0
        end if
    end if

end do

! output format statement from bpools.for in weps
2355       format (i6,1x,i5,1x,i4,1x,i3,1x,i4,1x,i2,30(1x,f10.5),1x,a30)

! for each single residue file, create a file with the minimum, average and maximum for each residue day
! and a minimum, median, maximum residue curve, where the curve has more days in this value than any other.

! note: minimum is the curve with the most values greater than another, not the maximum value for each day
do iname = 1, tname
    ! close for writing
    close( lunres(iname) )
    ! open for reading
    open(lunres(iname), file=rundirectory(1:len_trim(rundirectory)) // &
         resfile(iname)(1:len_trim(resfile(iname))), status="old", position="rewind", iostat=iosres)

    ! read header line
    read(lunres(iname),'(a)',iostat=iosres) inline(1)

    ! open stat file for writing
    open(lunrestat(iname), file=rundirectory(1:len_trim(rundirectory)) // &
         restatfile(iname)(1:len_trim(resfile(iname))), status="unknown", iostat=iosres)

    ! write header line
    write(lunrestat(iname),'(a)') inline(1)

    ! dynamically allocate min, max, mean and variance arrays of maxresday length
    if( allocated(resmin) ) then
        ! array is allocated, so deallocate before giving new dimension
        deallocate( resmin, resmax, resmean, resvar,ressd, resmnmsd, resmnpsd, cntresyr )
    end if
    allocate( resmin(30,maxresday(iname)), resmax(30,maxresday(iname)), &
              resmean(30,maxresday(iname)), resvar(30,maxresday(iname)), &
              ressd(30,maxresday(iname)), resmnmsd(30,maxresday(iname)), &
              resmnpsd(30,maxresday(iname)), cntresyr(maxresday(iname)) )

    ! intialize values
    do idx = 1, maxresday(iname)
        ! count of years
        cntresyr(idx) = 0
        ! for each data column
        do jdx = 1, 30
            ! set min to a high value so first value enters
            resmin(jdx,idx) = huge(resmin)
            ! set max to a low value so first value enters
            resmax(jdx,idx) = 0.0
            ! mean and variance values to zero
            resmean(jdx,idx) = 0.0
            resvar(jdx,idx) = 0.0
        end do
    end do

    do while( iosres .eq. 0 )
        ! read lines until selected residue name matches, return needed sub values
        call read_to_name ( lunres(iname), nstart, nend, uniquename(iname), inline(1),idata, rdata, iosres )

        ! idata(2) is the residue calendar day
        ! the data count increases each time a different resday is processed
        cntresyr(idata(2)) = cntresyr(idata(2)) + 1
        ! for each data column
        do jdx = 1, 30
            ! minimum value
            if( rdata(jdx) .lt. resmin(jdx,idata(2)) ) then
                ! new data value less less than old, make minimum
                resmin(jdx,idata(2)) = rdata(jdx)
            end if
            ! maximum value
            if( rdata(jdx) .gt. resmax(jdx,idata(2)) ) then
                ! new data value less less than old, make minimum
                resmax(jdx,idata(2)) = rdata(jdx)
            end if
            ! online variance algorithm (wikipedia cites Knuth who cites Welford)
            ! http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
            delta = rdata(jdx) - resmean(jdx,idata(2))
            resmean(jdx,idata(2)) = resmean(jdx,idata(2)) + delta / cntresyr(idata(2))
            ! this is not yet variance, just the accumulator
            resvar(jdx,idata(2)) = resvar(jdx,idata(2)) + delta * (rdata(jdx) - resmean(jdx,idata(2)))
        end do
    end do
    ! completed reading file, create output
    ! close for reading
    close( lunres(iname) )

    ! write stat output file
    ! write minimum block first
    do idx = 1, maxresday(iname)
        write(lunrestat(iname),2355) 0,idx,-4,0,0,0,(resmin(jdx,idx),jdx=1,30)
    end do

    ! write two blank lines for graphing purposes
    write(lunrestat(iname),*)
    write(lunrestat(iname),*)

    ! write maximum block
    do idx = 1, maxresday(iname)
        write(lunrestat(iname),2355) 0,idx,-3,0,0,0,(resmax(jdx,idx),jdx=1,30)
    end do

    ! write two blank lines for graphing purposes
    write(lunrestat(iname),*)
    write(lunrestat(iname),*)

    ! write mean block
    do idx = 1, maxresday(iname)
        write(lunrestat(iname),2355) 0,idx,-2,0,0,0,(resmean(jdx,idx),jdx=1,30)
    end do

    ! write two blank lines for graphing purposes
    write(lunrestat(iname),*)
    write(lunrestat(iname),*)

    ! compute the standard deviation from variance
    do idx = 1, maxresday(iname)
        ! for each data column
        do jdx = 1, 30
            ! here we divide by n-1 to get variance, then sqrt to get sd
            ressd(jdx,idx) = sqrt(resvar(jdx,idx)/(cntresyr(idx)-1))
            ! find the mean minus the standard deviation
            resmnmsd(jdx,idx) = resmean(jdx,idx) - ressd(jdx,idx)
            ! find the mean plus the standard deviation
            resmnpsd(jdx,idx) = resmean(jdx,idx) + ressd(jdx,idx)
        end do
    end do


    ! write mean minus standard deviation block
    do idx = 1, maxresday(iname)
        write(lunrestat(iname),2355) 0,idx,-1,0,0,0,(resmnmsd(jdx,idx),jdx=1,30)
    end do

    ! write two blank lines for graphing purposes
    write(lunrestat(iname),*)
    write(lunrestat(iname),*)

    ! write mean plus standard deviation block
    do idx = 1, maxresday(iname)
        write(lunrestat(iname),2355) 0,idx,0,0,0,0,(resmnpsd(jdx,idx),jdx=1,30)
    end do

    ! close for writing
    close( lunrestat(iname) )

    ! status output
    write(*,*) 'Wrote min, max, mean, SD for ', uniquename(iname)
end do

stop
end

! subroutine to read a given file until a line with a name is reached
subroutine read_to_name ( lun, nstart, nend, name, inline, idata, rdata, ios )

! arguments
integer :: lun ! logical unit number of file to be read
integer :: nstart, nend ! starting and ending character positions of name string
character*30 :: name ! name of residue to look for in line
character*512 :: inline ! input line read
integer :: idata(6) ! integer values on data line
real :: rdata(3) ! real values on data line
integer :: ios

! local values
character*30 :: linename ! name of residue found in line
integer :: idx

! read next line from file
read(lun,'(a)',iostat=ios) inline
! parse residue name
linename = inline(nstart:nend)
do while( (linename .ne. name) .and. (ios .eq. 0) )
    ! read next line from file
    read(lun,'(a)',iostat=ios) inline
    ! parse residue name
    linename = inline(nstart:nend)
end do

if( ios .eq. 0 ) then
    ! parse idata and rdata
    read(inline,*) (idata(idx),idx=1,6), (rdata(idx),idx=1,30)
end if

return
end