program wrf_comp

      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, max_dims, max_vars
      integer :: m, n, type, var_natts
      integer :: file, unlimid_ndx, unlim_ndx
      integer :: fs, fe
      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) = -1 ! 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 :: var_rms(:)
      real(r8), allocatable :: vector(:,:)
      real(r8), allocatable :: diff_1d(:)
      real(r8), allocatable :: matrix(:,:,:)
      real(r8), allocatable :: diff_2d(:,:)
      real(r8), allocatable :: matrix_3d(:,:,:,:)
      real(r8), allocatable :: diff_3d(:,:,:)
      real(r8), allocatable :: matrix_4d(:,:,:,:,:)
      real(r8), allocatable :: diff_4d(:,:,:,:)
      real(r8), allocatable :: matrix_5d(:,:,:,:,:,:)

      real(r8) :: eps, rms
      real(r8) :: max_val, thresh, threshold = 1.e-6_r8
      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=64) :: filespec(2)                ! netcdf filespec
      character(len=64) :: varname                    ! variable name
      character(len=64) :: vname(2)                   ! 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(:)

      write(*,*) 'wrf_comp: Enter number files to analysis'
      read(*,*) m
      fe = max( min(m,2),1 )
      fs = max( fe-1,1 )

      filespec = ' '
!------------------------------------------------------------------------------------------------------------
!	... Get and open the netcdf file
!------------------------------------------------------------------------------------------------------------
      do file = fs,fe
         write(*,*) 'wrf_comp: Enter netcdf filespec'
         read(*,*) filespec(file)
         ret = nf90_open( trim( filespec(file) ), NF90_NOWRITE, ncid(file) )
         if( ret /= NF90_NOERR ) then
            write(*,'(''wrf_comp: Failed to open '',a,''; error = '',i8)') trim( filespec(file) ),ret
            stop
         else
            write(*,'(''wrf_comp: 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(*,'(''wrf_comp: NF_INQ Failed for file '',a,'' ; error = '',i8)') trim( filespec(file) ), ret
	    stop
         end if
         if( ndims(file) <= 0 ) then
            write(*,*) 'wrf_comp: File ',trim(filespec(file)),' has No dimensions'
	    stop
         end if
         if( nvars(file) <= 0 ) then
            write(*,*) 'wrf_comp: File ',trim(filespec(file)),' has No variables'
	    stop
         end if
      end do
!------------------------------------------------------------------------------------------------------------
!	... Read threshold value
!------------------------------------------------------------------------------------------------------------
      write(*,*) 'wrf_comp: Enter threshold value'
      read(*,*) threshold
      write(*,*) 'wrf_comp: threshold value = ',threshold
!-----------------------------------------------------------------------------------------------------
!	... Now loop through and delineate the dimensions
!-----------------------------------------------------------------------------------------------------
      max_dims = maxval( ndims(:) )
      allocate( dim_size(max_dims,2), sizes(max_dims,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'wrf_comp: 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(*,*) 'wrf_comp: Failed to allocate dimension names array; error = ',ret
	 stop
      end if
      do file = fs,fe
         do m = 1,ndims(file)
            ret = nf90_inquire_dimension( ncid(file), m, name, dim_size(m,file) )
            if( ret /= NF90_NOERR ) then
               write(*,'(''wrf_comp: 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

!-----------------------------------------------------------------------------------------------------
!	... allocate arrays for variable properties
!-----------------------------------------------------------------------------------------------------
      write(*,*) ' '
      write(*,*) 'Variables per file = ',nvars(:)
      write(*,*) ' '

      max_vars = maxval( nvars(:) )
      allocate( var_type(max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'wrf_comp: Failed to allocate variable type array; error = ',ret
	 stop
      end if
      allocate( var_ndims(max_vars,2), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'wrf_comp: 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(*,*) 'wrf_comp: 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(*,*) 'wrf_comp: Failed to allocate variable names array; error = ',ret
	 stop
      end if
      allocate( perfect_match(max_vars), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'wrf_comp: Failed to allocate perfect_match array; error = ',ret
	 stop
      end if
      allocate( var_rms(max_vars), stat=ret )
      if( ret /= 0 ) then
         write(*,*) 'wrf_comp: Failed to allocate var_rms array; error = ',ret
	 stop
      end if
      var_rms(:) = 0._r8

      do file = fs,fe
         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(*,'(''wrf_comp: 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(*,*) 'wrf_comp: There are ',len_unlimdim,' indicies for the unlimited dimension'
          write(*,*) 'wrf_comp: Enter desired index number'
          read(*,*) unlim_ndx
          if( unlim_ndx < 1 .or. unlim_ndx > len_unlimdim ) then
            write(*,*) 'wrf_comp: Unlimited index number is out of range'
            stop
          else
            write(*,*) 'wrf_comp: Unlimited index number = ',unlim_ndx
          end if
        else
          unlim_ndx = len_unlimdim
        end if
      end if

      eps = 100._r8 * EPSILON( eps )
!------------------------------------------------------------------------------------------------------------
!	... Get variables to compare; name = stop or quit will exit the program
!------------------------------------------------------------------------------------------------------------
variable_loop : &
      do
	vndx(:) = 0
        do n = fs,fe
            write(*,*) ' '
            write(*,*) 'wrf_comp: Enter variable to compare'
            read(*,*) varname
	    if( trim( varname ) == 'quit' ) then
	       exit variable_loop
	    else
               var_ndx  = 1
               write(*,*) 'wrf_comp: Request comparison for ',trim( varname )
	    end if
	    do m = 1,nvars(1)
	      if( trim( varname ) == trim( var_names(m,1) ) ) then
	        vndx(n) = m
                vname(n) = trim( varname )
	        exit
	      end if
	    end do
        end do
!------------------------------------------------------------------------------------------------------------
!	... Check some basics
!------------------------------------------------------------------------------------------------------------
	 if( any( vndx(fs:fe) == 0 ) ) then
            write(*,*) ' '
            write(*,*) 'wrf_comp: ',trim( varname ),' not found in one dataset'
            write(*,*) ' '
	    exit
	 end if
	 if( var_type(vndx(fs),1) /= var_type(vndx(fe),1) ) then
            write(*,*) ' '
            write(*,*) 'wrf_comp: ',trim( varname ),' has different type'
            write(*,*) ' '
	    exit
	 end if
	 if( var_ndims(vndx(fs),1) /= var_ndims(vndx(fe),1) ) then
            write(*,*) ' '
            write(*,*) 'wrf_comp: ',trim( varname ),' has different dimensionality'
            write(*,*) ' '
	    exit
	 end if
	 do m = 1,var_ndims(vndx(1),1)
	    if( var_dimids(m,vndx(fs),1) /= var_dimids(m,vndx(fe),1) ) then
               write(*,*) ' '
               write(*,*) 'wrf_comp: ',trim( varname ),' has different dimension order'
               write(*,*) ' '
	       exit
	    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
             endif
           end do
!          dim_cnt = var_ndims(vndx(1),1) - 1
           dim_cnt = var_ndims(vndx(1),1)
         else
           dim_cnt = var_ndims(vndx(1),1)
         end if

! type = var_type(vndx(1),1)
! name = var_names(vndx(1),1)
! has_unlimited_dim = .false.
! dim_cnt = var_ndims(vndx(1),1)

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

         write(*,*) ' '
         write(*,*) 'dimension count = ',dim_cnt

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

               write(*,*) ' '
               write(*,*) 'asizes = ',asizes(:4)
	       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(*,'(''wrf_comp: Failed to allocate imatrix_4d; error = '',i8)') ret
	                stop
		     end if
                     write(*,*) '----------------------------------------------------------------------------'
		     do file = fs,fe
		        ret = nf90_get_var( ncid(file), vndx(file), imatrix_4d(:,:,:,:,file), count=asizes(1:4) )
                        if( ret /= NF90_NOERR ) then
                           write(*,'(''wrf_comp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
		        else
                           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(*,'(''wrf_comp: Failed to allocate matrix_4d; error = '',i8)') ret
	                stop
		     end if
                     write(*,*) '----------------------------------------------------------------------------'
		     do file = fs,fe
                        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(*,'(''wrf_comp: Failed to get variable '',a,'' ; error = '',i8)') trim(name),ret
	                   stop
		        else
                           write(*,'(''variable '',a,'' Min,Max val = '',1p,2e21.13)') &
                                trim(vname(file)),minval(matrix_4d(:,:,:,:,file)),maxval(matrix_4d(:,:,:,:,file))
			   negcnt  = count( matrix_4d(:,:,:,:,file) < 0._r8 )
			   zerocnt = count( matrix_4d(:,:,:,:,file) == 0._r8 )
			   tot_elements = PRODUCT( asizes(1:4) )
                           write(*,'(''                 Neg,zero %  = '',2f7.2)') &
                             100._r8*real(negcnt,r8)/real(tot_elements,r8),100.*real(zerocnt,r8)/real(tot_elements,r8)
                        end if
                     end do
                     if( fe == 1 ) then
                       matrix_4d(:,:,:,:,2) = 1._r8
                     endif
                     max_val = maxval(abs(matrix_4d(:,:,:,:,:)))
                     thresh = threshold*max_val
!                    where( matrix_4d(:,:,:,:,1) /= 0._r8 )
                     where( matrix_4d(:,:,:,:,1) >= thresh )
                       diff_4d(:,:,:,:) = 100._r8* abs( (matrix_4d(:,:,:,:,2) - matrix_4d(:,:,:,:,1))/matrix_4d(:,:,:,:,1) )
                     elsewhere
                       diff_4d(:,:,:,:) = 0._r8
                     endwhere
                     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
	                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(:),r8 ) / real( tot_elements,r8 )
		        write(*,*) ' '
		        write(*,'('' Diff bin counts'')')
		        write(*,'(8(1p,e10.3))') diff_levs(:)
		        write(*,'(8i10)') lev_cnts(:)
		        write(*,'(8f10.2)') diff_percents(:)
		     end if
		     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_3d == 0._r8 ),tot_elements
		        else
		           write(*,'(1x,i8,'' Zero elements out of '',i8)') count( diff_4d == 0._r8 ),tot_elements
		           write(*,'(a,'' Max diff at '',4i5,'' = '',1p,e22.15)') trim(vname(1)),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)
		        end if
			rms = sqrt( sum( diff_4d(:,:,:,:)**2 ) )/real( count( diff_4d /= 0._r8 ),r8 )
                        var_rms(var_ndx) = rms
		        write(*,*) 'Rms % error = ',rms
		     endif
		     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

      end program wrf_comp