module compare_vars_mod use filestruct, only : file_t, var_t, vdimsize, dim_t, verbose use prec, only : r4, r8, i4 use netcdf, only : nf90_char, nf90_int, nf90_double, nf90_float, nf90_get_var, nf90_max_dims, & nf90_inq_varid, nf90_get_att, nf90_noerr use utils, only : checknf90, get_dim_str, get_dimname_str implicit none logical :: ignoretime interface compute_var_stats ! TYPE real,double,int module procedure compute_var_stats_{TYPE} end interface interface get_rdiff_stats ! TYPE real,double module procedure get_rdiff_stats_{TYPE} end interface contains subroutine compare_vars(n,file, vtotal, ndiffs, nfilldiffs, vfailed, vnotfound ) integer, intent(in) :: n ! number of files to analyze (1 or 2) integer, intent(out) :: vtotal integer, intent(out) :: ndiffs ! number of fields with differences (not counting fields that differ only in the fill pattern) integer, intent(out) :: nfilldiffs ! number of fields with differences in fill pattern integer, intent(out) :: vfailed integer, intent(out) :: vnotfound type(file_t) :: file(n) double precision, pointer :: time(:,:) double precision :: tdiff ! start by making sure that times match integer :: ierr, ns1, ns2, vid1 type(var_t), pointer :: v1 integer :: i, t, nvars, t1, t2, nt integer :: idiff, ifilldiff, inotfound, vidnsteph integer, allocatable :: nsteph(:) character(len=132) :: dimstr type(dim_t), pointer :: udim real(r8), parameter :: timeepsilon = 1.e-9 ! time diff less than this considered insignificant vtotal = 0 vfailed = 0 vnotfound=0 if(n==2 .and. .not.ignoretime) then call checknf90(nf90_inq_varid(file(1)%fh, 'time', vid1)) ns1 = file(1)%dim(file(1)%unlimdimid)%dimsize if(n==2) then ns2 = file(2)%dim(file(2)%unlimdimid)%dimsize else ns2=1 end if allocate(time(max(ns1,ns2),2)) call checknf90(nf90_get_var(file(1)%fh, vid1, time(1:ns1,1))) if(n==2) then call checknf90(nf90_get_var(file(2)%fh, file(1)%var(vid1)%matchid, time(1:ns2,2))) end if if(verbose) then print *,'File 1 time: ', time(1:ns1,1) print *,'File 2 time: ', time(1:ns2,2) end if end if nvars = size(file(1)%var) udim => file(1)%dim(file(1)%unlimdimid) ndiffs = 0 nfilldiffs = 0 ! First look at variables which do not have unlimdim do i=1,nvars v1 => file(1)%var(i) if(.not. any(v1%dimids == file(1)%unlimdimid)) then call get_dimname_str(v1%ndims,v1%dimids,file(1)%dim,dimstr) write(6,140) trim(v1%name),trim(dimstr) vtotal = vtotal+1 ! initialize output arguments of compare_var in case compare_var doesn't get called idiff = 0 ifilldiff = 0 ierr = 0 inotfound = 0 select case(v1%xtype) case(nf90_int) call compare_var_int(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound) case(nf90_float) call compare_var_real(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound) case(nf90_double) call compare_var_double(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound) case(nf90_char) ierr = 1 ! call compare_var_char(file1,file2,i,v1%matchid) case default print *,'Type not recognized for variable: ', v1%name end select vnotfound = vnotfound+inotfound vfailed = vfailed+ierr ndiffs = ndiffs+idiff nfilldiffs = nfilldiffs + ifilldiff end if end do ierr = nf90_inq_varid(file(1)%fh, 'nsteph', vidnsteph) if(ierr == NF90_NOERR) then allocate(nsteph(udim%kount)) call checknf90(nf90_get_var(file(1)%fh, vidnsteph, nsteph)) end if do t=1,udim%dimsize,udim%kount t1 = t ! need to find mathing times - assumed for now t2 = t if(.not. ignoretime) then tdiff = abs(time(t1,1) - time(t2,2)) do while(t1<=ns1 .and. t2<= ns2 .and. tdiff > timeepsilon) if(time(t1,1) < time(t2,2)) then Write(6,*) 'Skipping a time sample on file 1' t1=t1+1 else if(time(t1,1) > time(t2,2)) then Write(6,*) 'Skipping a time sample on file 2' t2=t2+1 end if end do if(verbose) print *,__FILE__,__LINE__,tdiff,timeepsilon, t1, t2 if(tdiff< timeepsilon .and. t1/=t2) then Write(6,*) 'Found common timesteps:', t1, t2 else if(tdiff > timeepsilon) then Write(6,*) 'No matching time found.' vfailed = nvars return end if if(verbose) print *,__FILE__,__LINE__,time(t1,1),time(t2,2), t1, t2, time(:,:) end if if(allocated(nsteph)) then print *,'NSTEPH: ',nsteph(t) deallocate(nsteph) end if do i=1,nvars v1 => file(1)%var(i) if( any(v1%dimids == file(1)%unlimdimid)) then call get_dimname_str(v1%ndims,v1%dimids,file(1)%dim,dimstr) vtotal = vtotal+1 write(6,145) trim(v1%name),trim(dimstr), t1, t2 ! initialize output arguments of compare_var in case compare_var doesn't get called idiff = 0 ifilldiff = 0 ierr = 0 inotfound = 0 select case(v1%xtype) case(nf90_int) call compare_var_int(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound, (/t1, t2/)) case(nf90_float) call compare_var_real(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound, (/t1, t2/)) case(nf90_double) call compare_var_double(n,file,(/i,v1%matchid/), idiff, ifilldiff, ierr, inotfound, (/t1, t2/)) case(nf90_char) ierr = 1 ! call compare_var_char(file1,file2,i,v1(i)%matchid) case default print *,'Type not recognized for variable: ', v1%name end select vnotfound=vnotfound+inotfound vfailed = vfailed+ierr ndiffs = ndiffs+idiff nfilldiffs = nfilldiffs + ifilldiff end if end do end do 140 format(1x,a,3x,a) 145 format(1x,a,3x,a,' t_index = ',2i6) end subroutine compare_vars ! TYPE real,int,double subroutine compare_var_{TYPE}(n,file, vid, idiff, ifilldiff, ifail, inotfound, tindex) integer, intent(in) :: n type(file_t) :: file(2) integer, intent(in) :: vid(2) integer, intent(out) :: idiff ! 1 if diffs in field, 0 otherwise (0 if only diffs are in fill pattern) integer, intent(out) :: ifilldiff ! 1 if diffs in fill pattern, 0 otherwise ! (idiff & ifilldiff are both 1 if there are diffs in both the fill pattern and the valid values) integer, intent(out) :: ifail ! 1 if variable sizes differ, 0 otherwise integer, intent(out) :: inotfound ! 1 if variable on file 1 not found on file 2, 0 otherwise integer, optional :: tindex(2) integer :: s1, s2, l1(1), i, ierr {VTYPE}, pointer :: buf(:,:), vdiff(:) {VTYPE} :: fv1, fv2, rms, min_val(2), max_val(2), avgval(2), m1, rdmax {VTYPE} :: rdsum ! sum of relative differences {VTYPE} :: rdlogsum ! sum of negative log10 of relative differences {VTYPE} :: rdbar ! average of relative differences {VTYPE} :: rdlogbar ! rdlogsum normalized by number of non-zero differences integer :: t(2), n1, n2, min_loc(2), max_loc(2), spacelen integer :: start(NF90_MAX_DIMS,2), kount(NF90_MAX_DIMS,2), dsizes(NF90_MAX_DIMS,2) logical, pointer :: mask(:,:) integer :: diffcnt, rdmaxloc character(len=80) :: min_str(2), max_str(2), dmax_str, rdmax_str, space logical :: compare2 if(present(tindex)) then t = tindex else t = 1 end if compare2 = (n==2 .and. vid(2)>0) inotfound = 0 ifail = 0 ifilldiff = 0 idiff = 0 s1 = vdimsize(file(1)%dim, file(1)%var(vid(1))%dimids) if(verbose) print *,__FILE__,__LINE__,s1,file(1)%var(vid(1))%name if(compare2) then s2 = vdimsize(file(2)%dim, file(2)%var(vid(2))%dimids) if(s1 /= s2) then print *, 'Variable ',trim(file(1)%var(vid(1))%name),' sizes differ' ifail = 1 return end if end if n1 = size(file(1)%var(vid(1))%dimids) do i=1,n1 start(i,1) = file(1)%dim(file(1)%var(vid(1))%dimids(i))%start kount(i,1) = file(1)%dim(file(1)%var(vid(1))%dimids(i))%kount dsizes(i,1) = file(1)%dim(file(1)%var(vid(1))%dimids(i))%dimsize if(file(1)%var(vid(1))%dimids(i) == file(1)%unlimdimid) then start(i,1)=t(1) dsizes(i,1) = kount(i,1) end if end do allocate(buf(s1,n)) call checknf90(nf90_get_var(file(1)%fh, vid(1), buf(:,1), start(1:n1,1), kount(1:n1,1))) allocate(mask(s1,n)) ierr = nf90_get_att(file(1)%fh, vid(1), '_FillValue', fv1) if(ierr == NF90_NOERR) then mask(:,1) = (buf(:,1)/=fv1) else mask(:,1) = .true. end if if(n1>0) then call compute_var_stats(buf(:,1), s2, mask(:,1), min_loc(1), max_loc(1), min_val(1), max_val(1), avgval(1)) call get_dim_str(n1,translate_loc(n1,min_loc(1),start(1:n1,1),kount(1:n1,1),dsizes(1:n1,1)),min_str(1)) call get_dim_str(n1,translate_loc(n1,max_loc(1),start(1:n1,1),kount(1:n1,1),dsizes(1:n1,1)),max_str(1)) end if space = ' ' spacelen=1 if(n1>3) spacelen=(n1-3)*8 ! adjusts the output format if(compare2) then n2 = size(file(2)%var(vid(2))%dimids) if(n2/=n1) then print *,'WARNING variable ',trim(file(1)%var(vid(1))%name),' dims differ but total size is the same, will try to compare anyway' endif do i=1,n2 start(i,2) = file(2)%dim(file(2)%var(vid(2))%dimids(i))%start kount(i,2) = file(2)%dim(file(2)%var(vid(2))%dimids(i))%kount dsizes(i,2) = file(2)%dim(file(2)%var(vid(2))%dimids(i))%dimsize if(file(2)%var(vid(2))%dimids(i) == file(2)%unlimdimid) then start(i,2)=t(2) dsizes(i,2) = kount(i,2) end if end do call checknf90(nf90_get_var(file(2)%fh, vid(2), buf(:,2), start(1:n2,2), kount(1:n2,2))) ierr = nf90_get_att(file(2)%fh, vid(2), '_FillValue', fv2) if(ierr == NF90_NOERR) then mask(:,2) = (buf(:,2)/=fv2) else mask(:,2) = .true. end if if(n2>0) then call compute_var_stats(buf(:,2), s2, mask(:,2), min_loc(2), max_loc(2), min_val(2), max_val(2), avgval(2)) call get_dim_str(n2,translate_loc(n2,min_loc(2),start(1:n2,2),kount(1:n2,2),dsizes(1:n2,2)),min_str(2)) call get_dim_str(n2,translate_loc(n2,max_loc(2),start(1:n2,2),kount(1:n2,2),dsizes(1:n2,2)),max_str(2)) end if diffcnt=0 if(any(buf(:,1) /= buf(:,2))) then allocate(vdiff(s1)) ! Use the union of mask1 and mask2 if(any(mask(:,1) .neqv. mask(:,2))) then write(6,*) 'WARNING: Fill patterns differ between files' write(6,'(a,a8)') ' FILLDIFF ', file(1)%var(vid(1))%name ifilldiff = 1 mask(:,1) = (mask(:,1) .and. mask(:,2)) end if s2 = count(mask(:,1)) vdiff = abs(buf(:,1)-buf(:,2)) rms = sqrt(sum(vdiff**2,mask(:,1))/real(s2)) diffcnt = count((vdiff>0) .and. mask(:,1)) ! diffcnt==0 implies only diffs are in missing values if(diffcnt>0) then idiff = 1 m1 = maxval(vdiff, mask=mask(:,1)) l1 = maxloc(vdiff, mask=mask(:,1)) call get_dim_str(n1,translate_loc(n1,l1(1),start(1:n1,1),kount(1:n1,1),dsizes(1:n1,1)),dmax_str) #if ({ITYPE} != TYPEINT) call get_rdiff_stats(s1,buf(:,1),buf(:,2),vdiff,mask(:,1),rdsum, rdlogsum, rdmax, rdmaxloc) call get_dim_str(n1,translate_loc(n1,rdmaxloc,start(1:n1,1),kount(1:n1,1),dsizes(1:n1,1)),rdmax_str) #endif deallocate(vdiff) rdbar = rdsum / real(s2) rdlogbar = rdlogsum / real(diffcnt) if(n1==0) then write(6,803) s2, buf(1,1), buf(2,1) else write(6,800) diffcnt, s1, trim(max_str(1)),trim(min_str(1)), trim(dmax_str), trim(rdmax_str) #if ({ITYPE} == TYPEINT) write(6,903) s2, max_val(1), space(1:spacelen),min_val(1), m1, buf(l1(1),1), rdbar, buf(rdmaxloc,1), & count(mask(:,2)), max_val(2), space(1:spacelen),min_val(2), buf(l1(1),2), buf(rdmaxloc,2) write(6,810) s1, trim(max_str(2)), trim(min_str(2)) ! write(6,905) avgval(1), rms, rdbar, avgval(2), rdlogbar, log10(1./rdmax) write(6,'(a,a8,i10,/)') ' RMS ', file(1)%var(vid(1))%name, rms #else write(6,803) s2, max_val(1), space(1:spacelen),min_val(1), m1, buf(l1(1),1), rdbar, buf(rdmaxloc,1), & count(mask(:,2)), max_val(2), space(1:spacelen),min_val(2), buf(l1(1),2), buf(rdmaxloc,2) write(6,810) s1, trim(max_str(2)), trim(min_str(2)) write(6,805) avgval(1), rms, rdbar, avgval(2), rdlogbar, log10(1./rdmax) write(6,'(a,a8,1pe11.4,/)') ' RMS ', file(1)%var(vid(1))%name, rms #endif endif endif end if if(diffcnt==0) then ! no differences found if(n1>0) then write(6,810) s1, trim(max_str(1)),trim(min_str(1)) #if ({ITYPE} == TYPEINT) write(6,914) s2, max_val(1), space(1:spacelen),min_val(1), count(mask(:,2)),& max_val(2),space(1:spacelen),min_val(2) write(6,810) s1, trim(max_str(2)), trim(min_str(2)) write(6,915) avgval(1), avgval(2) #else write(6,814) s2, max_val(1), space(1:spacelen),min_val(1), count(mask(:,2)),& max_val(2),space(1:spacelen),min_val(2) write(6,810) s1, trim(max_str(2)), trim(min_str(2)) write(6,815) avgval(1), avgval(2) #endif endif end if else ! Single file analysis output if(n==2 ) then inotfound=1 write(6,*) 'Variable on file1: ',trim(file(1)%var(vid(1))%name),' not found on file2' end if write(6,810 ) s1, trim(max_str(1)),trim(min_str(1)) #if ({ITYPE} == TYPEINT) write(6, 820) s2, max_val(1),min_val(1) write(6, 821) avgval(1) #else write(6, 825) s2, max_val(1),min_val(1) write(6, 826) avgval(1) #endif end if deallocate(buf, mask) 800 format(3x,i8,1x,i8,2x,a,1x,a,1x,a,1x,a) 803 format(12x, i8,1x,1pe23.15,a,e23.15,e8.1,e23.15,e8.1,e23.15,/, & 12x, i8,1x, e23.15,a,e23.15,8x, e23.15,8x, e23.15) 805 format(10x,'avg abs field values: ',1pe23.15,4x,'rms diff:',e8.1, & 3x,'avg rel diff(npos): ',e8.1,/, & 10x,' ', e23.15,24x, & 'avg decimal digits(ndif): ',0p,f4.1,' worst: ',f4.1) 810 format(12x,i8,2x,a,1x,a) 811 format(12x, i8,2x,a,1x,a,11x,a,9x,a) 814 format(12x, i8,1x,e23.15,a,e23.15,/, & 12x, i8,1x,e23.15,a,e23.15) 815 format(10x,'avg abs field values: ',1pe23.15,/, & 10x,' ', e23.15) 820 format(12x,i8,1x,2i12) 821 format(12x,'avg abs field values: ',i12,/) 825 format(12x,i8,1x,1p2e23.15) 826 format(12x,'avg abs field values: ',1pe23.15,/) 903 format(12x, i8,1x,i8,a,5i8,/, & 12x, i8,1x, i8,a,i8,8x, i8,8x, i8) 905 format(10x,'avg abs field values: ',i8,4x,'rms diff:',i8, & 3x,'avg rel diff(npos): ',i8,/, & 10x,' ', i8,24x, & 'avg decimal digits(ndif): ',i8,' worst: ',i8) 914 format(12x, i8,1x,i8,a,i8,/, & 12x, i8,1x,i8,a,i8) 915 format(10x,'avg abs field values: ',i8,/, & 10x,' ', i8) end subroutine compare_var_{TYPE} ! TYPE real,double subroutine get_rdiff_stats_{TYPE} (s1, v1, v2, vdiff, mask, rdsum, rdlogsum, rdmax, loc) integer, intent(in) :: s1 {VTYPE}, intent(in) :: v1(:), v2(:), vdiff(:) logical, intent(in) :: mask(:) {VTYPE}, intent(out) :: rdsum, rdlogsum, rdmax integer, intent(out) :: loc {VTYPE} :: denom, rdiff(s1) integer :: i, iloc(1) rdiff=0 rdsum=0 rdlogsum=0 do i=1,s1 if(vdiff(i)>0) then denom = max(abs(v1(i)), abs(v2(i))) rdiff(i) = vdiff(i)/denom rdsum = rdsum+rdiff(i) rdlogsum = rdlogsum - log10(rdiff(i)) end if end do rdmax = maxval(rdiff) iloc = maxloc(rdiff) loc = iloc(1) end subroutine get_rdiff_stats_{TYPE} ! TYPE real,int,double subroutine compute_var_stats_{TYPE} (buf, nvalid, mask, min_loc, max_loc, min_val, max_val, avgval) {VTYPE}, intent(in) :: buf(:) logical, intent(in) :: mask(:) integer, intent(out) :: nvalid, min_loc, max_loc {VTYPE}, intent(out) :: min_val, max_val, avgval integer :: loc(2) nvalid = count(mask) if(nvalid>0) then loc(1:1) = maxloc(buf, mask=mask) loc(2:2) = minloc(buf, mask=mask) max_loc = loc(1) min_loc = loc(2) max_val = maxval(buf, mask=mask) min_val = minval(buf, mask=mask) avgval = sum(abs(buf),mask=mask)/real(nvalid) else max_loc=0 min_loc=0 max_val=0 min_val=0 avgval=0 end if end subroutine compute_var_stats_{TYPE} function translate_loc(ndims, loc, start, kount, dsize) integer, intent(in) :: ndims, loc, start(:), kount(:), dsize(:) integer :: translate_loc(ndims) integer :: i, tloc, tprod tprod = product(kount) if(loc>tprod) then write(6,*) 'ERROR in translate_loc: location ',loc,' exceeds array size',tprod stop end if if(ndims<1) then stop '0D array in translate_loc' endif translate_loc = 1 if(ndims==1) then translate_loc = loc else if(loc<=dsize(1)) then translate_loc(1) = loc else tloc = loc if(verbose) print *,__LINE__,loc,ndims,dsize(1:ndims) do i=ndims,1,-1 tprod = tprod/dsize(i) if(tloc>=tprod) then translate_loc(i) = tloc/tprod + start(i) tloc = tloc - (tloc/tprod)*tprod end if end do translate_loc(1) = translate_loc(1)-1 if(verbose) print *,__LINE__,translate_loc(1:ndims) end if end function translate_loc end module compare_vars_mod