program nccomp

      use netcdf

      implicit none

      integer, parameter :: r8 = selected_real_kind(12)

      integer :: ret                ! netcdf api return code
      integer :: dim_cnt, dimid
      integer :: negcnt, zerocnt, tot_elements
      integer :: max_dims, max_vars, min_vars
      integer, pointer :: min_ndx
      integer, target  :: minndx(1)
      integer :: m, n, type, var_natts
      integer :: file, unlimid_ndx, unlim_ndx
      integer :: var_ndx
      integer :: perf_cnt
      integer :: var_cnt
      integer :: len_unlimdim
      integer :: scalar_int(2)
      integer :: lev_cnts(8)
      integer :: ncid(2)            ! netcdf file id
      integer :: ndims(2)           ! netcdf file dimension count
      integer :: nvars(2)           ! netcdf file variable  count
      integer :: ngatts(2)          ! netcdf file global attribute count
      integer :: unlimdimid(2)      ! netcdf file unlimited dimension id
      integer :: vndx(2)            ! variable index
      integer :: asizes(10)         ! allocation sizes
      integer :: max_ind(10)        ! max diff indicies
      integer :: start(10)
      integer :: counts(10)

      integer, allocatable :: dim_size(:,:)
      integer, allocatable :: sizes(:,:)
      integer, allocatable :: var_ndims(:,:)
      integer, allocatable :: var_type(:,:)
      integer, allocatable :: var_dimids(:,:,:)

      integer, allocatable :: ivector(:,:)
      integer, allocatable :: imatrix(:,:,:)
      integer, allocatable :: imatrix_3d(:,:,:,:)
      integer, allocatable :: imatrix_4d(:,:,:,:,:)
      integer, allocatable :: imatrix_5d(:,:,:,:,:,:)

      real(r8), allocatable :: vector(:,:)
      real(r8), allocatable :: matrix(:,:,:)
      real(r8), allocatable :: diff_2d(:,:)
      real(r8), allocatable :: matrix_3d(:,:,:,:)
      real(r8), allocatable :: diff_3d(:,:,:)
      real(r8), allocatable :: diff_4d(:,:,:,:)
      real(r8), allocatable :: matrix_4d(:,:,:,:,:)
      real(r8), allocatable :: matrix_5d(:,:,:,:,:,:)

      real(r8) :: eps, rms
      real(r8) :: diff_percents(8)
      real(r8) :: diff_levs(8) = (/ 1.e-9_r8, 1.e-6_r8, 1.e-3_r8, 1.e-2_r8, 1.e-1_r8, 1._r8, 10._r8, 100._r8 /)

      character(len=132) :: filespec(2)                ! netcdf filespec
      character(len=64)  :: varname                    ! variable name
      character(len=64)  :: match_name                 ! variable name
      character(len=NF90_MAX_NAME) :: name
      character(len=NF90_MAX_NAME), allocatable :: var_names(:,:)
      character(len=NF90_MAX_NAME), allocatable :: dim_names(:,:)

      logical :: all_vars = .false.
      logical :: has_unlimited_dim
      logical, allocatable :: perfect_match(:)

      filespec = ' '
!------------------------------------------------------------------------------------------------------------
!	... Get and open the two netcdf files to compare
!------------------------------------------------------------------------------------------------------------
      do file = 1,2
         write(*,*) 'nccomp: Enter netcdf filespec'
         read(*,*) filespec(file)
         ret = nf90_open( trim( filespec(file) ), NF90_NOWRITE, ncid(file) )
         if( ret /= NF90_NOERR ) then
            write(*,'(''nccomp: Failed to open '',a,''; error = '',i8)') trim( filespec(file) ),ret
            stop
         else
            write(*,'(''nccomp: Opened '',a)') trim( filespec(file) )
         end if
         ret = nf90_inquire( ncid(file), ndims(file), nvars(file), ngatts(file), unlimdimid(file) )
         if( ret /= NF90_NOERR ) then
            write(*,'(''nccomp: NF_INQ Failed for file '',a,'' ; error = '',i8)') trim( filespec(file) ), ret
            stop
         end if
         if( ndims(file) <= 0 ) then
            write(*,*) 'nccomp: File ',trim(filespec(file)),' has No dimensions'
            stop
         end if
         if( nvars(file) <= 0 ) then
            write(*,*) 'nccomp: File ',trim(filespec(file)),' has No variables'
            stop
         end if
      end do
!-----------------------------------------------------------------------------------------------------
!	... delineate the dimensions
!-----------------------------------------------------------------------------------------------------
      max_dims = maxval( ndims(:) )
      allocate( dim_size(max_dims,2), sizes(max_dims,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate dimension size array; error = ',ret
         stop
      end if
      write(*,*) ' '
      allocate( dim_names(max_dims,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate dimension names array; error = ',ret
         stop
      end if
      do file = 1,2
        write(*,*) ' '
        do m = 1,ndims(file)
          ret = nf90_inquire_dimension( ncid(file), m, name, dim_size(m,file) )
          if( ret /= NF90_NOERR ) then
            write(*,'(''nccomp: NF_INQ_DIM Failed on dim id = '',i4,''; error = '',i4)') m, ret
            stop
          end if
          write(*,'(''Dimension '',i4,'' is named '',a,'' with size = '',i8)') m,trim(name),dim_size(m,file)
          dim_names(m,file) = name
        end do
      end do

!-----------------------------------------------------------------------------------------------------
!	... delineate the variables
!-----------------------------------------------------------------------------------------------------
      write(*,*) ' '
      write(*,*) 'Variables per file = ',nvars(:)
      write(*,*) ' '
      max_vars = maxval( nvars(:) )
      min_vars = minval( nvars(:) )
      min_ndx => minndx(1)
      minndx  = minloc( nvars(:) )

      allocate( var_type(max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate variable type array; error = ',ret
	 stop
      end if
      allocate( var_ndims(max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate variable dimension size array; error = ',ret
	 stop
      end if
      allocate( var_dimids(max_dims,max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate variable dimension id array; error = ',ret
	 stop
      end if
      allocate( var_names(max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate variable names array; error = ',ret
	 stop
      end if
      allocate( perfect_match(max_vars), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'nccomp: Failed to allocate perfect_match array; error = ',ret
	 stop
      end if

      do file = 1,2
        do m = 1,nvars(file)
          ret = nf90_inquire_variable( ncid(file), m, name, type, var_ndims(m,file), var_dimids(:,m,file), var_natts )
          if( ret /= NF90_NOERR ) then
            write(*,'(''nccomp: NF_INQ_VAR Failed on varaible id = '',i8,''; error = '',i8)') m, ret
	    stop
          end if
	  var_names(m,file) = name
	  var_type(m,file)  = type
        end do
      end do

!------------------------------------------------------------------------------------------------------------
!	... check for unlimited time dimension
!------------------------------------------------------------------------------------------------------------
      ret = nf90_inquire_dimension( ncid(1), unlimdimid(1), len=len_unlimdim )
      if( ret == 0 ) then
        if( len_unlimdim > 1 ) then
          write(*,*) ' '
          write(*,*) 'nccomp: There are ',len_unlimdim,' indicies for the unlimited dimension'
          write(*,*) 'nccomp: Enter desired index number'
          read(*,*) unlim_ndx
          if( unlim_ndx < 1 .or. unlim_ndx > len_unlimdim ) then
            write(*,*) 'nccomp: Unlimited index number is out of range'
	    stop
	  else
            write(*,*) 'nccomp: Unlimited index number = ',unlim_ndx
	  end if
	else
          unlim_ndx = len_unlimdim
	end if
      end if

      eps = 100._r8 * EPSILON( eps )
!------------------------------------------------------------------------------------------------------------
!	... Get variable to compare; name = quit will exit the program
!------------------------------------------------------------------------------------------------------------
variable_loop : &
      do
sng_var : &
         if( .not. all_vars ) then
            write(*,*) ' '
            write(*,*) 'nccomp: Enter variable to compare'
            read(*,*) varname
	    if( trim( varname ) == 'quit' ) then
	       exit
	    else if( trim(varname) == 'ALL' .or. trim(varname) == 'all' ) then
               all_vars = .true. 
               var_ndx  = 1
               perf_cnt = 0
	    else
               write(*,*) 'nccomp: Request comparison for ',trim( varname )
	    end if
            if( .not. all_vars ) then
              match_name = varname
	    else
              match_name = var_names(var_ndx,min_ndx)
	    endif
            call get_var_ndx
            if( all_vars ) then
               var_cnt   = 0
               perfect_match(:) = .false.
	    end if
         else sng_var
            var_ndx = var_ndx + 1
            if( var_ndx <= min_vars ) then
              match_name = var_names(var_ndx,min_ndx)
              call get_var_ndx
            else
               all_vars = .false.
               write(*,*) ' '
               write(*,'(1x,i5,'' matched variables out of '',i5)') perf_cnt,var_cnt
               if( perf_cnt /= var_cnt ) then
                  write(*,*) 'Following variables do not perfectly match'
                  do m = 1,min_vars
                     if( .not. perfect_match(m) ) then
                        write(*,'(a)') trim( var_names(m,min_ndx) )
                     end if
                  end do
               end if
               write(*,*) ' '
               cycle variable_loop
            end if
	 end if sng_var
	 if( any( vndx(:) == 0 ) ) then
            perfect_match(var_ndx) = .true.
            cycle variable_loop
	 end if
!------------------------------------------------------------------------------------------------------------
!	... Check some basics
!------------------------------------------------------------------------------------------------------------
	 name = var_names(vndx(1),1)
	 if( var_type(vndx(1),1) /= var_type(vndx(2),2) ) then
            write(*,*) ' '
            write(*,*) 'nccomp: ',trim( name ),' has different type'
            write(*,*) ' '
	    cycle variable_loop
	 end if
	 if( var_ndims(vndx(1),1) /= var_ndims(vndx(2),2) ) then
            write(*,*) ' '
            write(*,*) 'nccomp: ',trim( name ),' has different dimensionality'
            write(*,*) ' '
	    cycle variable_loop
	 end if
	 do m = 1,var_ndims(vndx(1),1)
	    if( var_dimids(m,vndx(1),1) /= var_dimids(m,vndx(2),2) ) then
               write(*,*) ' '
               write(*,*) 'nccomp: ',trim( name ),' has different dimension order'
               write(*,*) ' '
	       cycle variable_loop
	    end if
	 end do
	 type = var_type(vndx(1),1)
	 has_unlimited_dim = any( var_dimids(:var_ndims(vndx(1),1),vndx(1),1) == unlimdimid(1) )
	 if( has_unlimited_dim ) then
	    do m = 1,var_ndims(vndx(1),1)
	       if( var_dimids(m,vndx(1),1) == unlimdimid(1) ) then
		  unlimid_ndx = m
	          exit
	       end if
	    end do
	    dim_cnt = var_ndims(vndx(1),1) - 1
	 else
	    dim_cnt = var_ndims(vndx(1),1)
	 end if

         if( dim_cnt < 2  ) then
            perfect_match(var_ndx) = .true.
            cycle variable_loop
         end if
         var_cnt = var_cnt + 1

         select case( dim_cnt )
	    case( 0 )
	       select case( type )
	          case( nf90_int )
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), scalar_int(file:file), start=(/1/),count=(/1/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
                        end if
                     end do
	          case( nf90_float,nf90_double )
		     allocate( vector(1,2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate vector; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), vector(:,file), count=(/1/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
                        end if
                     end do
	       end select
	    case( 1 )
	       select case( type )
	          case( nf90_int )
		     dimid = var_dimids(1,vndx(1),1)
		     allocate( ivector(dim_size(dimid,1),2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate ivector; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), ivector(:,file), count=(/dim_size(dimid,1)/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
                        end if
                     end do
	          case( nf90_float,nf90_double )
		     dimid = var_dimids(1,vndx(1),1)
		     allocate( vector(dim_size(dimid,1),2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate vector; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), vector(:,file), count=(/dim_size(dimid,1)/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
                        end if
		     end do
	       end select
	    case( 2 )
	       if( .not. has_unlimited_dim ) then
	          do n = 1,var_ndims(vndx(1),1)
	             asizes(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       else
		  m = 0
	          do n = 1,var_ndims(vndx(1),1)
		     if( n /= unlimid_ndx ) then
			m = m + 1
	                asizes(m) = dim_size(var_dimids(n,vndx(1),1),1)
	             end if
	             counts(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       end if
	       select case( type )
	          case( nf90_int )
		     allocate( imatrix(asizes(1),asizes(2),2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate imatrix; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
			if( .not. has_unlimited_dim ) then
		          ret = nf90_get_var( ncid(file), vndx(file), imatrix(:,:,file), count=(/asizes(1),asizes(2)/) )
			else
			  start(:3) = 1
			  start(unlimid_ndx)  = unlim_ndx
			  counts(unlimid_ndx) = 1
		          ret = nf90_get_var( ncid(file), vndx(file), start=start(:3), count=counts(:3), values=imatrix(:,:,file) )
			end if
                        if( ret /= NF90_NOERR ) then
                          write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                  stop
		        else
                          write(*,*) ' '
                          write(*,'(''variable '',a,'' Min,Max val = '',2i10)') &
                               trim(name), MINVAL(imatrix(:,:,file)), maxval(imatrix(:,:,file))
			  negcnt = count( imatrix(:,:,file) < 0 )
                          zerocnt = count( imatrix(:,:,file) == 0 )
			  tot_elements = asizes(1)*asizes(2)
                          write(*,'(''                 Neg,zero %  = '',2f7.2)') &
                            100.*real(negcnt)/real(tot_elements),100.*real(zerocnt)/real(tot_elements)
                        end if
		     end do
		     if( all(imatrix(:,:,1) == imatrix(:,:,2)) ) then
		       write(*,'('' Matricies are a perfect match'')')
                       perf_cnt = perf_cnt + 1
                       perfect_match(var_ndx) = .true.
                     end if
		     deallocate( imatrix )
	          case( nf90_float,nf90_double )
		     allocate( matrix(asizes(1),asizes(2),2), &
		               diff_2d(asizes(1),asizes(2)), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate matrix; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
			if( .not. has_unlimited_dim ) then
		           ret = nf90_get_var( ncid(file), vndx(file), matrix(:,:,file), count=(/asizes(1),asizes(2)/) )
			else
			   start(:3) = 1
			   start(unlimid_ndx)  = unlim_ndx
			   counts(unlimid_ndx) = 1
		           ret = nf90_get_var( ncid(file), vndx(file), start=start(:3), count=counts(:3), values=matrix(:,:,file) )
			end if
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
		        else
                           write(*,*) ' '
                           write(*,'(''variable '',a,'' Min,Max val = '',1p,2e21.13)') &
                                trim(name), MINVAL(matrix(:,:,file)), maxval(matrix(:,:,file))
			   negcnt = count( matrix(:,:,file) < 0. )
                           zerocnt = count( matrix(:,:,file) == 0. )
			   tot_elements = asizes(1)*asizes(2)
                           write(*,'(''                 Neg,zero %  = '',2f7.2)') &
                             100.*real(negcnt)/real(tot_elements),100.*real(zerocnt)/real(tot_elements)
                        end if
                     end do
		     where( matrix(:,:,1) /= 0._r8 )
			diff_2d(:,:) = 100._r8* abs( (matrix(:,:,2) - matrix(:,:,1))/matrix(:,:,1) )
		     elsewhere
			diff_2d(:,:) = 0._r8
		     endwhere
	             lev_cnts(1) = count( diff_2d(:,:) <= diff_levs(1) )
		     do n = 2,8
			lev_cnts(n) = count( diff_2d(:,:) <= diff_levs(n) .and. diff_2d(:,:) > diff_levs(n-1) )
		     end do
		     diff_percents(:) = 100._r8 * real( lev_cnts(:) ) / real( tot_elements )
		     write(*,*) ' '
		     write(*,'('' Diff bin counts'')')
		     write(*,'(8i10)') lev_cnts(:)
		     write(*,'(8f10.2)') diff_percents(:)
                     max_ind(1:2) = maxloc( diff_2d(:,:) )
		     if( diff_2d(max_ind(1),max_ind(2)) == 0._r8 ) then
		        write(*,'('' Matricies are a perfect match'')')
                        perf_cnt = perf_cnt + 1
                        perfect_match(var_ndx) = .true.
		     else if( all( diff_2d(:,:) <= eps ) ) then
		        write(*,'('' Matricies are a nearly perfect match'')')
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_2d == 0._r8 ),tot_elements
			rms = sqrt( sum( diff_2d(:,:)**2 ) )/real( count( diff_2d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     else
		        write(*,'('' Matricies do not match'')')
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_2d == 0._r8 ),tot_elements
		        write(*,'('' Max diff at '',2i5,'' = '',1p,e22.15)') max_ind(1:2),diff_2d(max_ind(1),max_ind(2))
!	                write(*,'('' Values = '',z16,1x,z16)') matrix(max_ind(1),max_ind(2),1) , &
!                                                              matrix(max_ind(1),max_ind(2),2)
		        write(*,'('' Values = '',1p,g22.15,1x,g22.15)') matrix(max_ind(1),max_ind(2),1) , &
                                                               matrix(max_ind(1),max_ind(2),2)
			rms = sqrt( sum( diff_2d(:,:)**2 ) )/real( count( diff_2d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     end if
		     deallocate( matrix, diff_2d )
	       end select
	    case( 3 )
	       if( .not. has_unlimited_dim ) then
	          do n = 1,var_ndims(vndx(1),1)
	             asizes(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       else
		  m = 0
	          do n = 1,var_ndims(vndx(1),1)
		     if( n /= unlimid_ndx ) then
			m = m + 1
	                asizes(m) = dim_size(var_dimids(n,vndx(1),1),1)
	             end if
	             counts(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       end if
	       select case( type )
	          case( nf90_int )
		     allocate( imatrix_3d(asizes(1),asizes(2),asizes(3),2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate imatrix_3d; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), imatrix_3d(:,:,:,file), count=(/asizes(1),asizes(2),asizes(3)/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
		        else
                           write(*,*) ' '
                           write(*,'(''variable '',a,'' Min,Max val = '',2i8)') &
                               trim(name),MINVAL(imatrix_3d(:,:,:,file)),maxval(imatrix_3d(:,:,:,file))
                        end if
		     end do
		     deallocate( imatrix_3d )
	          case( nf90_float,nf90_double )
		     allocate( matrix_3d(asizes(1),asizes(2),asizes(3),2), diff_3d(asizes(1),asizes(2),asizes(3)),stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate matrix_3d; error = '',i8)') ret
	                stop
		     end if
		     do file = 1,2
			if( .not. has_unlimited_dim ) then
		           ret = nf90_get_var( ncid(file), vndx(file), matrix_3d(:,:,:,file), count=(/asizes(1),asizes(2),asizes(3)/) )
			else
			   start(:4) = 1
			   start(unlimid_ndx) = unlim_ndx
			   counts(unlimid_ndx) = 1
		           ret = nf90_get_var( ncid(file), vndx(file), start=start(:4), count=counts(:4), values=matrix_3d(:,:,:,file) )
			end if
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
		        else
                           write(*,*) ' '
                           write(*,'(''variable '',a,'' Min,Max val = '',1p,2e21.13)') &
                                   trim(name),minval(matrix_3d(:,:,:,file)),maxval(matrix_3d(:,:,:,file))
			   negcnt = count( matrix_3d(:,:,:,file) < 0. )
			   zerocnt = count( matrix_3d(:,:,:,file) == 0. )
			   tot_elements = asizes(1)*asizes(2)*asizes(3)
                           write(*,'(''                 Neg,zero count = '',2i8)') negcnt,zerocnt
                           write(*,'(''                 Neg,zero %  = '',2f7.2)') &
                             100._r8*real(negcnt)/real(tot_elements),100._r8*real(zerocnt)/real(tot_elements)
                        end if
                     end do
		     where( matrix_3d(:,:,:,1) /= 0._r8 )
			diff_3d(:,:,:) = 100._r8* abs( (matrix_3d(:,:,:,2) - matrix_3d(:,:,:,1))/matrix_3d(:,:,:,1) )
		     elsewhere
			diff_3d(:,:,:) = 0._r8
		     endwhere
	             lev_cnts(1) = count( diff_3d(:,:,:) <= diff_levs(1) )
		     do n = 2,8
			lev_cnts(n) = count( diff_3d(:,:,:) <= diff_levs(n) .and. diff_3d(:,:,:) > diff_levs(n-1) )
		     end do
		     write(*,'('' count > 100% = '',i8)') count( diff_3d(:,:,:) > diff_levs(8) )
		     write(*,'('' count > min err = '',i8)') count( diff_3d(:,:,:) >= diff_levs(1) )
		     diff_percents(:) = 100._r8 * real( lev_cnts(:) ) / real( tot_elements )
		     write(*,*) ' '
		     write(*,'('' Diff bin counts'')')
		     write(*,'(8(1p,e10.3))') diff_levs(:)
		     write(*,'(8i10)') lev_cnts(:)
		     write(*,'(8f10.2)') diff_percents(:)
                     max_ind(1:3) = maxloc( diff_3d(:,:,:) )
		     if( diff_3d(max_ind(1),max_ind(2),max_ind(3)) == 0._r8 ) then
		        write(*,'('' Matricies are a perfect match'')')
                        perf_cnt = perf_cnt + 1
                        perfect_match(var_ndx) = .true.
		     else if( all( diff_3d(:,:,:) <= eps ) ) then
		        write(*,'('' Matricies are a nearly perfect match'')')
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_3d == 0._r8 ),tot_elements
			rms = sqrt( sum( diff_3d(:,:,:)**2 ) )/real( count( diff_3d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     else
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_3d == 0._r8 ),tot_elements
		        write(*,'('' Max diff at '',3i5,'' = '',1p,e22.15)') max_ind(1:3),diff_3d(max_ind(1),max_ind(2),max_ind(3))
!	                write(*,'('' Values = '',z16,1x,z16)') matrix_3d(max_ind(1),max_ind(2),max_ind(3),1) , &
!                                                              matrix_3d(max_ind(1),max_ind(2),max_ind(3),2)
		        write(*,'('' Values = '',1p,g22.15,1x,g22.15)') matrix_3d(max_ind(1),max_ind(2),max_ind(3),1) , &
                                                               matrix_3d(max_ind(1),max_ind(2),max_ind(3),2)
			rms = sqrt( sum( diff_3d(:,:,:)**2 ) )/real( count( diff_3d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     end if
		     deallocate( matrix_3d, diff_3d )
	       end select
	    case( 4 )
	       if( .not. has_unlimited_dim ) then
	          do n = 1,var_ndims(vndx(1),1)
	             asizes(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       else
		  m = 0
	          do n = 1,var_ndims(vndx(1),1)
		     if( n /= unlimid_ndx ) then
			m = m + 1
	                asizes(m) = dim_size(var_dimids(n,vndx(1),1),1)
	             end if
	             counts(n) = dim_size(var_dimids(n,vndx(1),1),1)
	          end do
	       end if
	       select case( type )
	          case( nf90_int )
		     allocate( imatrix_4d(asizes(1),asizes(2),asizes(3),asizes(4),2), stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate imatrix_4d; error = '',i8)') ret
	                stop 'Alloc err'
		     end if
		     do file = 1,2
		        ret = nf90_get_var( ncid(file), vndx(file), imatrix_4d(:,:,:,:,file), &
                                            count=(/asizes(1),asizes(2),asizes(3),asizes(4)/) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop 'NCD err'
		        else
                           write(*,*) ' '
                           write(*,'(''variable '',a,'' Min,Max val = '',2i8)') &
                               trim(name),minval(imatrix_4d(:,:,:,:,file)),maxval(imatrix_4d(:,:,:,:,file))
                        end if
		     end do
		     deallocate( imatrix_4d )
	          case( nf90_float,nf90_double )
		     allocate( matrix_4d(asizes(1),asizes(2),asizes(3),asizes(4),2), &
                               diff_4d(asizes(1),asizes(2),asizes(3),asizes(4)),stat=ret )
                     if( ret /= 0 ) then
                        write(*,'(''nccomp: Failed to allocate matrix_4d; error = '',i8)') ret
	                stop 'Alloc err'
		     end if
		     do file = 1,2
			if( .not. has_unlimited_dim ) then
		           ret = nf90_get_var( ncid(file), vndx(file), matrix_4d(:,:,:,:,file), &
                                               count=(/asizes(1),asizes(2),asizes(3),asizes(4)/) )
			else
			   start(:5) = 1
			   start(unlimid_ndx) = unlim_ndx
			   counts(unlimid_ndx) = 1
		           ret = nf90_get_var( ncid(file), vndx(file), start=start(:5), count=counts(:5), values=matrix_4d(:,:,:,:,file) )
			end if
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''nccomp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop 'NCD err'
		        else
                           write(*,*) ' '
                           write(*,'(''variable '',a,'' Min,Max val = '',1p,2e21.13)') &
                                   trim(name),minval(matrix_4d(:,:,:,:,file)),maxval(matrix_4d(:,:,:,:,file))
			   negcnt  = count( matrix_4d(:,:,:,:,file) < 0. )
			   zerocnt = count( matrix_4d(:,:,:,:,file) == 0. )
			   tot_elements = product( asizes(1:4) )
                           write(*,'(''                 Neg,zero count = '',2i8)') negcnt,zerocnt
                           write(*,'(''                 Neg,zero %  = '',2f7.2)') &
                             100._r8*real(negcnt)/real(tot_elements),100._r8*real(zerocnt)/real(tot_elements)
                        end if
                     end do
		     where( matrix_4d(:,:,:,:,1) /= 0._r8 )
			diff_4d(:,:,:,:) = 100._r8* abs( (matrix_4d(:,:,:,:,2) - matrix_4d(:,:,:,:,1))/matrix_4d(:,:,:,:,1) )
		     elsewhere
			diff_4d(:,:,:,:) = 0._r8
		     endwhere
	             lev_cnts(1) = count( diff_4d(:,:,:,:) <= diff_levs(1) )
		     do n = 2,8
			lev_cnts(n) = count( diff_4d(:,:,:,:) <= diff_levs(n) .and. diff_4d(:,:,:,:) > diff_levs(n-1) )
		     end do
		     write(*,'('' count > 100% = '',i8)') count( diff_4d(:,:,:,:) > diff_levs(8) )
		     write(*,'('' count > min err = '',i8)') count( diff_4d(:,:,:,:) >= diff_levs(1) )
		     diff_percents(:) = 100._r8 * real( lev_cnts(:) ) / real( tot_elements )
		     write(*,*) ' '
		     write(*,'('' Diff bin counts'')')
		     write(*,'(8(1p,e10.3))') diff_levs(:)
		     write(*,'(8i10)') lev_cnts(:)
		     write(*,'(8f10.2)') diff_percents(:)
                     max_ind(1:4) = maxloc( diff_4d(:,:,:,:) )
		     if( diff_4d(max_ind(1),max_ind(2),max_ind(3),max_ind(4)) == 0._r8 ) then
		        write(*,'('' Matricies are a perfect match'')')
                        perf_cnt = perf_cnt + 1
                        perfect_match(var_ndx) = .true.
		     else if( all( diff_4d(:,:,:,:) <= eps ) ) then
		        write(*,'('' Matricies are a nearly perfect match'')')
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_4d == 0._r8 ),tot_elements
			rms = sqrt( sum( diff_4d(:,:,:,:)**2 ) )/real( count( diff_4d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     else
		        write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_4d == 0._r8 ),tot_elements
		        write(*,'('' Max diff at '',4i5,'' = '',1p,e22.15)') max_ind(1:4),diff_4d(max_ind(1),max_ind(2),max_ind(3),max_ind(4))
		        write(*,'('' Values = '',1p,g22.15,1x,g22.15)') matrix_4d(max_ind(1),max_ind(2),max_ind(3),max_ind(4),1) , &
                                                               matrix_4d(max_ind(1),max_ind(2),max_ind(3),max_ind(4),2)
			rms = sqrt( sum( diff_3d(:,:,:)**2 ) )/real( count( diff_4d /= 0._r8 ) )
		        write(*,*) 'Rms % error = ',rms
		     end if
		     deallocate( matrix_4d, diff_4d )
	       end select
	 end select
      end do variable_loop

      if( allocated( dim_size ) ) then
         deallocate( dim_size )
      end if
      if( allocated( dim_names ) ) then
         deallocate( dim_names )
      end if
      if( allocated( var_names ) ) then
         deallocate( var_names )
      end if
      if( allocated( var_ndims ) ) then
         deallocate( var_ndims )
      end if
      if( allocated( var_dimids ) ) then
         deallocate( var_dimids )
      end if
      if( allocated( perfect_match ) ) then
         deallocate( perfect_match )
      end if

      do m = 1,2
         ret = nf90_close( ncid(m) )
      end do

      CONTAINS

      subroutine get_var_ndx

      vndx(:) = 0
      do file = 1,2
        do m = 1,nvars(file)
	  if( trim( match_name ) == trim( var_names(m,file) ) ) then
            vndx(file) = m
            exit
	  end if
	end do
      end do

      end subroutine get_var_ndx

      end program nccomp