program compare ! ! compares all of the double precision arrays in two netcdf files. ! ! this is hopefully a generically useful utility, and does not ! link in any gryffin-specific code. ! ! this code won't notice if there are new variables in the second netcdf ! file that aren't in the first. Also, it doesn't check integer variables. ! ! Note: Some variables can appear to have a large relative ! error, but are actually zero to within roundoff. ! This compare program accounts for this in many cases, such as ! arrays calculated by FFT's, by using the avg value of an ! array as a lower bound in calculating the relative error. ! But there are still other cases where a large relative error ! might be reported. An example of this is in runs with adiabatic ! electrons, where the particle flux is by definition zero ! (though small roundoff errors may make it slightly non-zero). ! Thus, although the relative error for fluxi and fluximn might ! be reported as order unity in some cases, fluxi and ! fluximn should still be negligible in comparison with ! the ion heat flux Q_i (=wpfx). Use ncdump to check this. ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! use command_line implicit none include 'netcdf.inc' character filename_nc1*80 character filename_nc2*80 integer ncid1, ncid2, ilen character var_name*80 character var_name2*80 character, dimension(:), allocatable :: var_names2*80 character dim_name*80 integer var_type integer ndims integer dimids(NF_MAX_VAR_DIMS) integer dim_lengths(NF_MAX_VAR_DIMS) integer dim_lengths2(NF_MAX_VAR_DIMS) integer natts integer var_id, var_id2 integer status integer nvars, nvars2, dim_i, var_size real*8, dimension(:), allocatable :: var1, var2, err_abs_a, err_rel_a real avg_var, err_abs, err_abs_all integer err_loc(1) real err_rel,tol_rel, tol_rel_default real tol_abs, tol_abs_default real err_rel_all character var_name_all*80 character tol_word*80 integer err_loc_all(1) integer icount, ierr ! errors must be larger than larger tol_rel as a relative error and ! larger than tol_abs in order to be printed. tol_rel_default=1.e-10 tol_abs_default=0.0 tol_rel=tol_rel_default tol_abs=tol_abs_default err_rel_all=0.0 ! determine the names of the two netcdf files: icount=iargc() if(icount >= 2) then call cl_getarg(1,filename_nc1,ilen,ierr) call cl_getarg(2,filename_nc2,ilen,ierr) else goto 99 endif ! If the third argument exists, get the relative error tolerance if(icount >= 3) then call cl_getarg(3,tol_word,ilen,ierr) read(tol_word,*) tol_rel write(*,*) "using relative error tolerance ", tol_rel endif ! If the 4th argument exists, get the absolute error tolerance if(icount >= 4) then call cl_getarg(4,tol_word,ilen,ierr) read(tol_word,*) tol_abs write(*,*) "using absolute error tolerance ", tol_abs endif status = nf_open(filename_nc1, nf_nowrite, ncid1) if(status /= nf_noerr) then write(*,*) "ERROR: can't open file ", trim(filename_nc1) goto 99 endif status = nf_open(filename_nc2, nf_nowrite, ncid2) if(status /= nf_noerr) then write(*,*) "ERROR: can't open file ", trim(filename_nc2) goto 99 endif ! Get a list of all variables in file2: status=nf_inq_nvars(ncid2,nvars2) allocate(var_names2(nvars2)) do var_id=1,nvars2 status = nf_inq_varname(ncid2,var_id, var_names2(var_id)) enddo ! Loop over every variable in file1: status=nf_inq_nvars(ncid1,nvars) do var_id=1,nvars status = nf_inq_var(ncid1, var_id, var_name, var_type, & ndims, dimids, natts) if(var_type == nf_double) then do dim_i = 1, ndims status=nf_inq_dim(ncid1, dimids(dim_i), dim_name, & dim_lengths(dim_i)) enddo var_size=1 do dim_i = 1, ndims var_size=var_size*dim_lengths(dim_i) enddo ! Find file1's var_name in file2: do var_id2=1, nvars2 if(var_name == var_names2(var_id2)) goto 55 enddo write(*,*) "ERROR: variable ", trim(var_name) write(*,*) "can't be found in ",trim(filename_nc2) goto 85 55 continue status = nf_inq_var(ncid2, var_id2, var_name2, var_type, & ndims, dimids, natts) do dim_i = 1, ndims status=nf_inq_dim(ncid2, dimids(dim_i), dim_name, & dim_lengths2(dim_i)) if(dim_lengths(dim_i) /= dim_lengths2(dim_i)) then write(*,*) "ERROR: variable ", trim(var_name) write(*,*) " is not the same size in the two files." goto 85 endif enddo ! Read the data as one long vector. It would be too much work ! right now to add the needed code for various rank matrices. ! IDL could be used to write a more complete and robust compare ! routine quickly. allocate(var1(var_size),var2(var_size)) allocate(err_abs_a(var_size),err_rel_a(var_size)) status=nf_get_var_double(ncid1,var_id,var1) status=nf_get_var_double(ncid2,var_id2,var2) ! calculate avg value of the variables: avg_var=sum(abs(var1)+abs(var2))/(2*size(var1)) avg_var=max(avg_var,tiny(avg_var)) ! calculate an array of the errors, relative to the ! local value of the array, or to the average value of the ! array. This is because FFT's can introduce an error that ! is small in the RMS sense, though locally for some value of k ! it might appear large. ! err_abs_a=abs(var1-var2) err_rel_a=err_abs_a/max(abs(var1),abs(var2),avg_var) err_loc=maxloc(err_rel_a) err_rel=err_rel_a(err_loc(1)) err_abs=err_abs_a(err_loc(1)) if(err_rel >= err_rel_all) then err_rel_all=err_rel err_abs_all=err_abs var_name_all=var_name err_loc_all=err_loc endif if(err_rel > tol_rel .and. err_abs > tol_abs) then write(*,*) 'Variable ',trim(var_name) write(*,*) 'Maximum relative error =', err_rel write(*,*) 'corresponding abs err =', err_abs write(*,*) "Occurs at ", err_loc, "'th position" write(*,*) ' ' endif deallocate(var1, var2, err_abs_a, err_rel_a) 85 continue endif enddo if(err_rel_all .eq. 0.0) then write(*,*) "No differences found" stop endif write(*,*) "Maximum relative error =", err_rel_all write(*,*) "Found in variable ", trim(var_name_all) write(*,*) " at ", err_loc_all, "'th position" stop 99 write(*,*) & 'Usage: compare file1.nc file2.nc [rel_tolerance] [abs_tolerance]' write(*,*) ' will compare all double precision arrays in two netcdf files' write(*,*) ' and report all arrays with relative error larger than' write(*,*) ' rel_tolerance (default value = ', tol_rel_default, ')' write(*,*) ' and absolute error larger than' write(*,*) ' abs_tolerance (default value = ', tol_abs_default, ')' write(*,*) ' ' write(*,*) ' see compare.f90 source code for more info.' stop end program compare