! given a WEPS run directory, parse the dec#.btmp files into unique crop files
! this will not work given files generated by a calibration run.
! details should only be generated after calibration is complete

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*10 :: 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*10 :: resfile(maxname) ! array of output file names
integer :: lunres(maxname) ! logical unit number for residue output file
character*10 :: 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)

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), previdata(6)
real :: rdata(30), prevrdata(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(:,:)
real :: delta

integer :: cntresblock(maxname), iblock  ! count of residue blocks for each name, counter to place values in array

real, allocatable :: endDDstand(:), endDDflat(:), endDDbgL10(:)
integer, allocatable :: endresyrstand(:), endresyrflat(:), endresyrbgL10(:)

! 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), 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"
! terminates with residue year that is in last (catch all) pool when simulation ends.
! incomplete residue years are not put into this file
! 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)

! 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
! and the number of "residue years" or "blocks" for each residue
! to be used to create arrays of the last day of the block for Decomposition Days

! 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), status="unknown", iostat=ios(iname))
    ! write header line
    write(lunres(iname),'(a)') inline(1)(1:nend)

    ! intialize maxresday
    maxresday(iname) = 0

    ! zero out block count array
    cntresblock(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(1)(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),*)

            ! update block count for this residue
            cntresblock(iname) = cntresblock(iname) + 1
            !write(*,*) 'Wrote block ',cntresblock(iname),' of ',uniquename(iname),' to ',resfile(iname)

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

    if( (iosall .ne. 0) .and. (iname .ne. 0) ) then
        ! reached end of file while name block was being output
        ! regardless, this is end of block so write two blank line for graphing
        write(lunres(iname),*)
        write(lunres(iname),*)

        ! update block count for this residue
        cntresblock(iname) = cntresblock(iname) + 1
    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-SD, average, average+SD and maximum for each residue 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), 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), status="unknown", iostat=iosres)

    ! write header line
    write(lunrestat(iname),'(a)') inline(1)(1:len_trim(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, &
                    endDDstand, endDDflat, endDDbgL10, endresyrstand, endresyrflat, endresyrbgL10 )
    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)), &
              endDDstand(cntresblock(iname)), endDDflat(cntresblock(iname)), &
              endDDbgL10(cntresblock(iname)), endresyrstand(cntresblock(iname)), &
              endresyrflat(cntresblock(iname)), endresyrbgL10(cntresblock(iname)) )

    ! intialize values
    previdata(2) = 0
    iblock = 0
    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
        if( idata(2) .gt. 0 ) then
            ! 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 if

        ! create array of last day values for Decomp Days: stand, flat, below ground Layer 10
        if( idata(2) .lt. previdata(2) ) then
            ! calendar day age this line is less than calendar day age previous line
            ! this means bump into a new age block so previous value was last block value
            iblock = iblock + 1
            ! Decomp Days standing
            endDDstand(iblock) = prevrdata(1)
            ! the residue year for this value
            endresyrstand(iblock) = previdata(3)
        
            ! Decomp Days flat
            endDDflat(iblock) = prevrdata(2)
            ! the residue year for this value
            endresyrflat(iblock) = previdata(3)

            ! Decomp Days below ground layer 10
            endDDbgL10(iblock) = prevrdata(3)
            ! the residue year for this value
            endresyrbgL10(iblock) = previdata(3)
        end if
        
        ! roll data into prevday values
        previdata = idata
        prevrdata = rdata

    end do
    ! completed reading file, create output

    ! 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
            if( cntresyr(idx) .gt. 1 ) then
                ressd(jdx,idx) = sqrt(resvar(jdx,idx)/(cntresyr(idx)-1))
            else
                ressd(jdx,idx) = 0.0
            end if
            ! 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

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

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

    ! find the min year, max year, median year, and quartiles
    ! sort array for Decomp Days standing
    call sort2( cntresblock(iname), endDDstand, endresyrstand )
    ! search in file and output block for min year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrstand(1) )
    ! search in file and output block for max year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrstand(cntresblock(iname)) )
    ! search in file and output block for median year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrstand(cntresblock(iname)/2) )
    ! search in file and output block for lower quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrstand(cntresblock(iname)/4) )
    ! search in file and output block for upper quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrstand((cntresblock(iname)*3)/4) )

    ! sort array for Decomp Days flat
    call sort2( cntresblock(iname), endDDflat, endresyrflat )
    ! search in file and output block for min year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrflat(1) )
    ! search in file and output block for max year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrflat(cntresblock(iname)) )
    ! search in file and output block for median year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrflat(cntresblock(iname)/2) )
    ! search in file and output block for lower quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrflat(cntresblock(iname)/4) )
    ! search in file and output block for upper quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrflat((cntresblock(iname)*3)/4) )

    ! sort array for Decomp Days bgL10
    call sort2( cntresblock(iname), endDDbgL10, endresyrbgL10 )
    ! search in file and output block for min year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrbgL10(1) )
    ! search in file and output block for max year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrbgL10(cntresblock(iname)) )
    ! search in file and output block for median year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrbgL10(cntresblock(iname)/2) )
    ! search in file and output block for lower quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrbgL10(cntresblock(iname)/4) )
    ! search in file and output block for upper quartile year
    call read_write_resyr_block( lunres(iname), lunrestat(iname), endresyrbgL10((cntresblock(iname)*3)/4) )

    ! close for reading
    close( lunres(iname) )
    ! close for writing
    close( lunrestat(iname) )

    ! status output
    write(*,*) 'Wrote min, max, median, quartile curves: DD(stand, flat, bgL10) 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(30) ! 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)
else
    ! triggers end for block count
    idata(2) = 0
end if

return
end

! subroutine to sort two arrays together, based on values in first one
! from Numerical Recipes, 1986
subroutine sort2( n, rval, ival )

! arguments
integer :: n  ! dimension of arrays
integer :: ival(n)
real :: rval(n)

! local values
integer :: level, retire
integer :: iival
real :: rrval
integer :: idx, jdx

level = n/2 + 1
retire = n

10 continue
if( level .gt. 1 ) then
    level = level - 1
    rrval = rval(level)
    iival = ival(level)
else
    rrval = rval(retire)
    iival = ival(retire)
    rval(retire) = rval(1)
    ival(retire) = ival(1)
    retire = retire - 1
    if( retire .eq. 1 ) then
        rval(1) = rrval
        ival(1) = iival
        return
    end if
end if
idx = level
jdx = level + level

20 if( jdx .le. retire ) then
    if( jdx .lt. retire ) then
        if( rval(jdx) .lt. rval(jdx+1) ) then
            jdx = jdx + 1
        end if
    end if
    if( rrval .lt. rval(jdx) )then
        rval(idx) = rval(jdx)
        ival(idx) = ival(jdx)
        idx = jdx
        jdx = jdx + jdx
    else
        jdx = retire + 1
    end if
    go to 20
end if
rval(idx) = rrval
ival(idx) = iival
go to 10
end

subroutine read_write_resyr_block( lunin, lunout, selresyr )

! arguments
integer :: lunin   ! logical unit for input file
integer :: lunout  ! logical unit for output file
integer :: selresyr ! selected residue year block to be copied from input to output file

! local values
integer :: iosin, iosout  ! input and output file status
character*512 :: inline   ! input line read
integer :: daysim, resday, resyr ! first three values in data line

 
! start at beginning of file
rewind( lunin )
! read header line
read( lunin,'(a)',iostat=iosin ) inline
! read second line
read( lunin,'(a)',iostat=iosin ) inline

do while( iosin .eq. 0 )

    !check for data value
    if( inline(1:20) .ne. '                    ' ) then
        ! parse data for line
        read(inline,*) daysim, resday, resyr
        if( resyr .eq. selresyr ) then
            ! this is the residue year we want, write line to file
            write( lunout,'(a)', iostat=iosout ) inline(1:len_trim(inline))
        end if
    end if

    if( resyr .gt. selresyr ) then
        iosin = 1
    else
        ! read another line
        read( lunin,'(a)',iostat=iosin ) inline
    end if

end do

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

return
end