c#######################################################################
c     Reads input netcdf file, looks for deformation variable
c     Returns defmax=0.0 if no deformation var found
c#######################################################################
      subroutine ReadDeform(ncid1in,h1a0g,h2a0g,
     &                            dx0,dy0,nx0g,ny0g,
     &                            ma,na,h1a,h2a,
     &                            ddeform,defmax)
      implicit none
      include 'netcdf.inc'

      integer*4 ncid1in, iret
      integer*4 DEF_id, GLON_id, GLAT_id
      integer*4 dx0, dy0
      integer*4 nx0g, ny0g, ma, na
      integer defdimids(20)
      character*50 defdimname
      real*8 h1a(ma),h2a(na)
      real*8 h1a0g(nx0g),h2a0g(ny0g)
      real*8 ddeform(nx0g,ny0g)
      real*8 defmax
      logical deformfound
      integer i,j

      defmax=0.0
      deformfound = .false.
      if(nf_inq_varid(ncid1in,'DEFORMATION',DEF_id).eq.NF_NOERR ) then
        deformfound = .true.
      endif
      if(.not.deformfound.AND.
     &   (nf_inq_varid(ncid1in,'deformation',DEF_id).eq.NF_NOERR)) then
        deformfound = .true.
      endif
      if (.not.deformfound.AND.
     &    (nf_inq_varid(ncid1in,'DEFLC',DEF_id).eq.NF_NOERR)) then
         deformfound = .true.
      endif
      if (.not.deformfound.AND.
     &    (nf_inq_varid(ncid1in,'deflc',DEF_id).eq.NF_NOERR)) then
         deformfound = .true.
      endif
      if (deformfound) then
        write(3,*) ' '
        write(3,*) ' Found deformation in IC source file'
        iret = nf_inq_vardimid(ncid1in,DEF_id,defdimids)
        iret = nf_inq_dimname(ncid1in,defdimids(1),defdimname)
        if (iret .ne. NF_NOERR) then
           write(3,*) 'cant get defdimname'
           return
        endif
        if ( defdimname(1:8).eq.'GRID_LON' .or.
     &       defdimname(1:8).eq.'grid_lon') then
          GLON_id=defdimids(1)
          GLAT_id=defdimids(2)
          iret = nf_inq_dimlen(ncid1in,GLON_id,dx0)
          iret = nf_inq_dimlen(ncid1in,GLAT_id,dy0)
          iret = nf_inq_varid(ncid1in,defdimname,GLON_id)
c.........Ferret appends stuff to the dim var name
          iret = nf_inq_dimname(ncid1in,defdimids(2),defdimname)
          iret = nf_inq_varid(ncid1in,defdimname,GLAT_id)
          if (iret .ne. NF_NOERR) then
             write(3,*) 'cant read def grid axis id'
             return
          endif
        else
          GLON_id=defdimids(2)
          GLAT_id=defdimids(1)
          iret = nf_inq_dimlen(ncid1in,GLON_id,dx0)
          iret = nf_inq_dimlen(ncid1in,GLAT_id,dy0)
          iret = nf_inq_varid(ncid1in,defdimname,GLAT_id)
c.........Ferret appends stuff to the dim var name
          iret = nf_inq_dimname(ncid1in,defdimids(1),defdimname)
          iret = nf_inq_varid(ncid1in,defdimname,GLAT_id)
          if (iret .ne. NF_NOERR) then
             write(3,*) 'cant read def grid axis id'
             return
          endif
        endif
        call readdeformnc(ncid1in, DEF_id,dx0,dy0,nx0g,ny0g,ddeform)
        iret = nf_get_var_double(ncid1in,GLON_id,h1a0g)
        call check_err(iret)
        iret = nf_get_var_double(ncid1in,GLAT_id,h2a0g)
        call check_err(iret)
         do i=1,dx0
          do j=1,dy0
c              write(3,*) i, j, h1a0g(i), h2a0g(j), ddeform(i,j)
c              call flush(3)
            if (h1a0g(i).ge.h1a(1) .AND. h1a0g(i).le.h1a(ma) .AND.
     &          h2a0g(j).ge.h2a(1) .AND. h2a0g(j).le.h2a(na)) then
              if (abs(ddeform(i,j)).gt.defmax) then
                defmax=ddeform(i,j)
              endif
            endif
          enddo
        enddo
      endif

      end

c#######################################################################
c     ReadDeformNC reads ddeform array from netcdf forcing file 
c     deformation variable, using dummy array of proper size
c#######################################################################
      subroutine ReadDeformNC(ncid1in,DEF_id,dx0,dy0,nx0g,ny0g,ddeform)
      implicit none
      include 'netcdf.inc'

      integer ncid1in, iret, i, j, dx0, dy0, nx0g, ny0g, DEF_id
      real*8 deftmp(dx0,dy0)
      real*8 ddeform(nx0g,ny0g)

      iret = nf_get_var_double(ncid1in,DEF_id,deftmp)

      do j=1,dy0
         do i=1,dx0
            ddeform(i,j) = deftmp(i,j)
         end do
      end do

      end
