!***********************************************************************

      program convert_POP_T42

!-----------------------------------------------------------------------
!
!     this program remaps dummy atmosphere data from a T42 grid to 
!     the x3prime ocean grid
!
!-----------------------------------------------------------------------

      implicit none

      include 'netcdf.inc'

!-----------------------------------------------------------------------
!
!     variables to be mapped
!
!-----------------------------------------------------------------------

      integer :: nx, ny, nz, i, j, k, nbits, rec_length

      integer, dimension(4) ::
     &  beg,       ! beginning point for netCDF array
     &  cnt        ! count for each axis of netCDF array 

      real, dimension(:,:,:), allocatable ::
     &      field  ! remapped field to be written
      real, dimension(:,:), allocatable ::
     &      field2d ! temporary for checking NaNs

      double precision, dimension(:,:), allocatable :: work

      character*80 var_name, in_file, out_file

      real :: land_val, new_land_val

      logical :: sp_output
!-----------------------------------------------------------------------
!
!     variables required for mapping
!
!-----------------------------------------------------------------------

      integer :: ! netCDF ids
     &         ncstat, nc_file_id, nc_field_id

!-----------------------------------------------------------------------
!    get netcdf file info
!-----------------------------------------------------------------------

      write(*,*)' input nx,ny,nz'
      read(*,*)nx,ny,nz
      write(*,*)' input name of netcdf file'
      read(*,'(a80)')in_file
      write(*,*)' input name of output file'
      read(*,'(a80)')out_file
      write(*,*)' input variable name'
      read(*,'(a80)')var_name
      write(*,*)' input 1 for single precision output'
      read(*,*)i
      if (i == 1) then 
        sp_output = .true.
      else
        sp_output = .false.
      endif

!-----------------------------------------------------------------------
!
!     read in data and remap
!
!-----------------------------------------------------------------------

      print *,'reading file ',in_file

      ncstat = nf_open(in_file,NF_NOWRITE,nc_file_id)
      call netcdf_error_handler(ncstat)
 
      ncstat = nf_inq_varid(nc_file_id,var_name,nc_field_id)
      call netcdf_error_handler(ncstat)

!-----------------------------------------------------------------------
!   allocate 3d fields
!-----------------------------------------------------------------------

      allocate( field(nx, ny, nz), field2d(nx,ny) )

      beg(1) = 1
      beg(2) = 1
      beg(3) = 1
      beg(4) = 1
      cnt(1) = nx
      cnt(2) = ny
      cnt(3) = 1
      cnt(4) = nz

      ncstat = nf_get_vara_real(nc_file_id,nc_field_id,beg,cnt,
     &                                field)
      call netcdf_error_handler(ncstat)

      ncstat = nf_close(nc_file_id)

!-----------------------------------------------------------------------

!     where ( abs(field) > 1.e15) field = 0.
      land_val = field(1,1,1)
      write(*,*)' land value appears to be ',land_val
      write(*,*)' input 1 for new land value'
      read(*,*) i
      if(i == 1) then
         write(*,*)' input new land value'
         read(*,*) new_land_val
         land_val = new_land_val
      endif
      do k = 1, nz
        field2d = land_val
        do j = 1, ny
          do i = 1, nx
            if ( abs(field(i,j,k)) < 1.e30) field2d(i,j) = field(i,j,k)
          enddo
        enddo
        field(:,:,k) = field2d
      enddo
      write(*,*)'min, max = ',minval(field),maxval(field)

      if ( sp_output ) then
         inquire (iolength=rec_length) field(:,:,1)
      else
         allocate( work(nx, ny) )
         inquire (iolength=rec_length) work
      endif
      open(11, file=out_file, status='unknown',
     &     form='unformatted', access='direct', recl=rec_length)
      do k= 1,nz
        if ( sp_output ) then
          write(11,rec=k)field(:,:,k)
        else
          work = field(:,:,k)
          write(11,rec=k)work
        endif
      enddo

      close(11)
      write(*,*)' wrote file: ',out_file

!-----------------------------------------------------------------------

      end program convert_POP_T42

!***********************************************************************

      subroutine netcdf_error_handler(istat)

!-----------------------------------------------------------------------
!
!     This routine provides a simple interface to netCDF error message
!     routine.
!
!-----------------------------------------------------------------------

      implicit none

      include 'netcdf.inc'

      integer, intent(in) :: 
     &    istat   ! integer status returned by netCDF function call

!-----------------------------------------------------------------------

      if (istat /= NF_NOERR) then
        print *,'Error in netCDF: ',nf_strerror(istat)
        stop
      endif

!-----------------------------------------------------------------------

      end subroutine netcdf_error_handler

!***********************************************************************