
      program compare_snaps
      implicit none
      integer, parameter :: unit1 = 31
      integer, parameter :: unit2 = 32
      integer ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte

      open(unit=unit1,file="snap_out_001",status="old",form="unformatted",err=22)
      open(unit=unit2,file="snap_gpu_001",status="old",form="unformatted",err=23)

      read(unit=unit1) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
      write(*,*)ids,ide,jds,jde,kds,kde
      write(*,*)ims,ime,jms,jme,kms,kme
      write(*,*)its,ite,jts,jte,kts,kte
      read(unit=unit2) ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
      write(*,*)ids,ide,jds,jde,kds,kde
      write(*,*)ims,ime,jms,jme,kms,kme
      write(*,*)its,ite,jts,jte,kts,kte

      call read_and_output3( 'rthraten', unit1, 88,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte )
      call read_and_output3( 'rthraten', unit2, 98,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte )
      call read_and_output2( 'gsw', unit1, 88, ims,ime,jms,jme,its,ite,jts,jte )
      call read_and_output2( 'gsw', unit2, 98, ims,ime,jms,jme,its,ite,jts,jte )

      close(unit1)
      close(unit2)
      stop
22    write(*,*)'error opening ',unit1
      stop
23    write(*,*)'error opening ',unit2
      stop
      end

      subroutine read_and_output3 ( varname,iu,ou,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte )
      implicit none
      integer iu,ou
      integer ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
      character*(*) varname
      integer i,j,k
      real s3(ims:ime,kms:kme,jms:jme) 
      read(iu)s3
      write(ou,*)ite-its+1,jte-jts+1,' ',trim(varname)
      do j = jts,jte
      do i = its,ite
        write(ou,*)s3(i,1,j)
      enddo
      enddo
      return
      end

      subroutine read_and_output2 ( varname,iu, ou, ims,ime,jms,jme,its,ite,jts,jte )
      implicit none
      integer iu,ou
      integer ims,ime,jms,jme,its,ite,jts,jte
      character*(*) varname
      integer i,j,k
      real s2(ims:ime,jms:jme) 
      read(iu)s2
      write(ou,*)ite-its+1,jte-jts+1,' ',trim(varname)
      do j = jts,jte
      do i = its,ite
        write(ou,*)s2(i,j)
      enddo
      enddo
      return
      end
      


      

      

