program remap

integer, parameter :: gx=4000, gy=2000
real, parameter :: glon=-179.955, glat=89.955,  gres=0.09

!Parameter for CONUS 1km
!integer, parameter :: dx=5824, dy=2500
!real, parameter :: dlon=-124.995, dlat=24.805,  dres=0.01

!Parameter for CONUS 4km
!integer, parameter :: dx=1456, dy=625
!real, parameter :: dlon=-124.98, dlat=24.82,  dres=0.04

!Parameter for NA 8km
!integer, parameter :: dx=1050, dy=700
!real, parameter :: dlon=-136.46, dlat=4.54,  dres=0.08

!Parameter for CONUS 2km
integer, parameter :: dx=2912, dy=1250
real, parameter :: dlon=-124.99, dlat=24.83,  dres=0.02
  
integer, parameter :: window=1

integer :: startx, starty
real:: clat, clon
real :: sum, samplen
real, allocatable:: inputmap(:,:)
real :: inputvalue

real :: outputmap(dx,dy)
real :: outputmap_rev(dx,dy)
real :: lcpixel

integer :: status, i, j, m, n

  CHARACTER(LEN=400) :: input_filename     ! Input file name with full path
  CHARACTER(LEN=400) :: output_filename     ! Input file name with full path

  CALL GETARG(1, input_filename)
  
  CALL GETARG(2, output_filename)

allocate(inputmap(gx,gy),stat=status)

open(12,file=input_filename,form='unformatted',access='direct', &
      convert='little_endian', recl=gx*4)
do i=1, gy
read(12,rec=i) inputmap(:,i)
enddo
close(12)

outputmap(:,:) = -9999.
outputmap_rev(:,:) = -9999.

  do j=1, dy
  clat = dlat + (j-1)*dres
  do i=1, dx
  clon = dlon + (i-1)*dres

  startx = NINT((clon-glon)/gres + 1.)
  starty = NINT((glat-clat)/gres + 1.)

  inputvalue=inputmap(startx,starty)
  
  if( (inputvalue >= 200) .and. (inputvalue <= 370) ) then

   outputmap_rev(i,j)=inputvalue

  else
  ! filling gap
  sum = 0.
  samplen = 0.

  do n=0-window, window
  do m=0-window, window
  if( (startx+m >=1) .and. (startx+m <= dx) .and. &
      (starty+n >=1) .and. (starty+n <= dy)) then          
	
	inputvalue=inputmap(startx+m,starty+n)
	if( (inputvalue >= 200) .and. (inputvalue <= 370) ) then
	sum=sum+inputvalue
	samplen = samplen + 1
	endif
  endif
  enddo
  enddo

  if( samplen >= 1) then
	outputmap_rev(i,j) = sum / samplen 
  else
   outputmap_rev(i,j)= -9999
  endif
   
  endif 
  !end filling gap

  enddo
  enddo

!  do j=1, dy
!  do i=1, dx

!      outputmap(i,j) = outputmap_rev(i, dy-j+1)

!  enddo
!  enddo

open(50,file=output_filename,form='unformatted',access='direct', &
     convert='little_endian', recl=dx*dy*4)
write(50,rec=1) outputmap_rev
close(50)

!for test
!open(55,file='testreadin.bin',form='unformatted',access='direct', &
!     convert='little_endian', recl=gx*4)
!do i=1, gy
!write(55,rec=i) inputmap(:,i)
!enddo
!close(55)

DEALLOCATE(inputmap)

end program