program calclim

  IMPLICIT NONE

  CHARACTER(LEN=400) :: climmeanfile, climstdfile, climfd, climnfile
  CHARACTER(LEN=400) :: buffer

  INTEGER :: NCOL 
  INTEGER :: NROW
  integer :: i, j, n, iNum, ios
  CHARACTER(LEN=400) :: tempfile
  CHARACTER(LEN=4) :: yyyystr
  real :: diff
  character(len=1024) :: filename
  character(len=1024) :: filepile(150)
  integer :: iunit

  REAL, allocatable :: input(:,:,:)
  REAL, allocatable :: etmean(:,:)
  REAL, allocatable :: etstd(:,:)
  real, allocatable :: etsum(:,:)
  real, allocatable :: etsum_sqr(:,:)
  real, allocatable :: etcount(:,:)

  CALL GETARG(1, buffer)
  read(buffer, *) NCOL
  CALL GETARG(2, buffer)
  read(buffer, *) NROW
  CALL GETARG(3, climfd)
  CALL GETARG(4, climmeanfile)
  CALL GETARG(5, climstdfile)
  CALL GETARG(6, climnfile)

!------------
!Get all files within the range
    ! Open a temporary file to store the list of files
    iunit = 10 
    open(unit=iunit, file='filelist.txt', status='replace', action='write')
    call system('ls /tmp > filelist.txt')
    close(iunit)

    ! Reopen the file for reading
    open(unit=iunit, file='filelist.txt', status='old', action='read')

    ! Read the file names from the file
    n = 0
    filepile(:) = " "
    do
        read(iunit, '(A)', iostat=ios) filename
        if (ios /= 0) exit ! Exit the loop if end of file or error
        print *, trim(filename)
        n = n+1
        filepile(n) = trim(climfd) // '/' // trim(filename)
    end do
    ! Close the file
    close(iunit)

    ! Optionally, delete the temporary file
    call system('rm filelist.txt')

  print*, "totalN = ", n

  allocate (input(n, NCOL,NROW))
  allocate (etmean(NCOL,NROW))
  allocate (etstd(NCOL,NROW))
  allocate (etsum(NCOL,NROW))
  allocate (etsum_sqr(NCOL,NROW))
  allocate (etcount(NCOL,NROW))

  etcount(:,:)=0.0
  etsum(:,:)=0.0
  etsum_sqr(:,:)=0.0  
  etmean(:,:)=-9999.;
  etstd(:,:)=-9999.

  input(:,:,:) = -9999;

  do iNum = 1, n

  tempfile = filepile(iNum)
  print*, trim(tempfile)
  open(12,file=tempfile,action='READ', status='OLD',form='unformatted', &
  convert='little_endian', access='direct',recl=NCOL*NROW*4,iostat=ios)

  if(ios == 0) then
  read(12,rec=1) input(iNum,:,:)
  close(12)

  do j=1, NROW
  do i=1, NCOL
     if ( (input(iNum,i,j) >= 0.) .and. (input(iNum,i,j) <= 5.)) then
     etsum(i,j)=etsum(i,j)+input(iNum,i,j)
     etcount(i,j)=etcount(i,j)+1
     endif
  enddo
  enddo
  else
  print*, 'missing: ', trim(tempfile)
  endif

  enddo  !enddo iNum

  do j=1, NROW
  do i=1, NCOL

    if (etcount(i,j) .gt. 1) then
    etmean(i,j)=etsum(i,j)/etcount(i,j)
    endif
  enddo
  enddo

  do j=1, NROW
  do i=1, NCOL

    diff = 0.
    do iNum = 1, n
      if ( (etcount(i,j) .gt. 1) .and. (input(iNum,i,j) >= 0.) .and.  &
           (input(iNum,i,j) <= 5) ) then
          diff = input(iNum,i,j) - etmean(i,j)
          etsum_sqr(i,j) = etsum_sqr(i,j) + diff*diff
      endif
    enddo

    if (etcount(i,j) .gt. 1) then
       etstd(i,j) = sqrt( etsum_sqr(i,j)/etcount(i,j) )
    endif

  enddo
  enddo 

  open(50,file=climmeanfile,action='WRITE', status='REPLACE',form='unformatted', &
        convert='little_endian',access='direct',recl=NCOL*NROW*4)
  write(50,rec=1) etmean
  close(50)

  open(51,file=climstdfile,action='WRITE', status='REPLACE',form='unformatted', &
      convert='little_endian',access='direct',recl=NCOL*NROW*4)
  write(51,rec=1) etstd
  close(51)

  open(52,file=climnfile,action='WRITE', status='REPLACE',form='unformatted', &
      convert='little_endian',access='direct',recl=NCOL*NROW*4)
  write(52,rec=1) etcount
  close(52)


end program