filestruct.F90

module filestruct
  use netcdf
  implicit none
  type dim_t
     integer :: dimsize
     integer :: start, kount  ! used for user requested dimension subsetting
     character(len=nf90_MAX_NAME) ::name = ''
  end type dim_t

  type var_t
     integer :: matchid
     integer :: ndims
     integer :: natts
     integer, pointer :: dimids(:)
     integer :: xtype
     character(len=nf90_MAX_NAME) ::name = ''
  end type var_t

  type file_t
     integer :: fh
     integer :: natts
     type(dim_t), pointer :: dim(:)
     type(var_t), pointer :: var(:)
     integer :: unlimdimid
   contains
     procedure :: has_unlimited_dim  ! logical function; returns true if this file has an unlimited dimension
  end type file_t

  logical :: verbose

contains
  logical function has_unlimited_dim(file)
    ! Returns true if this file has an unlimited dimension
    class(file_t), intent(in) :: file

    if (file%unlimdimid == -1) then
       has_unlimited_dim = .false.
    else
       has_unlimited_dim = .true.
    end if
  end function has_unlimited_dim

  subroutine init_file_struct( file, dimoptions )

    type(file_t) :: file
    type(dim_t), optional :: dimoptions(:)
    integer :: ndims, nvars
    integer :: dimids(NF90_MAX_DIMS)
    integer :: i, ierr, docnt, n1, n2
    integer :: j, start, kount
    character(len=NF90_MAX_NAME) :: name, dname
    ierr= nf90_inquire(file%fh, ndims, nvars, file%natts, file%unlimdimid)

    allocate(file%dim(ndims))
    allocate(file%var(nvars))


    do i=1,ndims
       ierr = nf90_inquire_dimension(file%fh, i, file%dim(i)%name, file%dim(i)%dimsize)
       file%dim(i)%start=1
       if(i==file%unlimdimid) then
          file%dim(i)%kount=1
       else
          file%dim(i)%kount=file%dim(i)%dimsize
       end if
    end do

    if(present(dimoptions)) then
       docnt = size(dimoptions)
       do j=1,docnt
          start = dimoptions(j)%start
          kount = dimoptions(j)%kount
          name = dimoptions(j)%name
          n1 = len_trim(name)
          do i=1,ndims
             dname = file%dim(i)%name
             n2 = len_trim(dname)
             if(name(1:n1).eq.dname(1:n2) ) then


                if((start > 0) .and. (start < file%dim(i)%dimsize)) then
                   file%dim(i)%start = start
                else
                   write(6,*) 'Command line start value for dim ',name(1:n1),&
                        ' out of bounds, expected 1-',file%dim(i)%dimsize,' got: ',start
                   stop
                end if
                if(kount > 0 .and. start+kount <= file%dim(i)%dimsize) then
                   file%dim(i)%kount = kount
                else if(kount == -1) then
                   file%dim(i)%kount = file%dim(i)%dimsize-file%dim(i)%start+1
                else
                   write(6,*) 'Command line count value for dim ',name(1:n1),&
                        ' out of bounds, expected 1-',file%dim(i)%dimsize-file%dim(i)%start+1,' got: ',kount
                   stop

                endif
                write(6,*) 'Setting dimension bounds for dim ',name(1:n1),file%dim(i)%start,file%dim(i)%kount

                exit
             end if
          end do
       end do
    end if

    do i=1,nvars
       file%var(i)%matchid=-1
       ierr = nf90_inquire_variable(file%fh, i, file%var(i)%name, file%var(i)%xtype, file%var(i)%ndims, dimids, &
            file%var(i)%natts)
       allocate(file%var(i)%dimids(file%var(i)%ndims))
       file%var(i)%dimids = dimids(1:file%var(i)%ndims)
    end do


  end subroutine init_file_struct


  subroutine compare_metadata(file1, file2, vid)
    type(file_t) :: file1, file2
    integer, optional, intent(in) :: vid

    integer :: id1, id2, natts1, natts2

    integer :: i, ierr
    character(len=NF90_MAX_NAME) :: attname
    integer :: atttype, attlen

    real, pointer :: attreal1(:), attreal2(:)
    double precision, pointer :: attdouble1(:),attdouble2(:)
    integer, pointer :: attint1(:),attint2(:)
    integer, parameter :: maxstrlen=32767
    character(len=maxstrlen) :: attchar1, attchar2
    logical :: found


    if(present(vid)) then
       id1 = vid
       id2 = file1%var(id1)%matchid
       ierr = nf90_inquire_variable(file1%fh, id1, nAtts=natts1)
       ierr = nf90_inquire_variable(file2%fh, id2, nAtts=natts2)
    else
       id1 = NF90_GLOBAL
       id2 = NF90_GLOBAL
       natts1 = file1%natts
       natts2 = file2%natts
    end if

    do i=1,natts1
       found = .true.
       attname = ''
       ierr = nf90_inq_attname(file1%fh, id1, i, attname)
       ierr = nf90_inquire_attribute(file1%fh, id1, trim(attname), atttype, attlen)

       select case(atttype)
       case(nf90_char)
          if (attlen > maxstrlen) then
             stop 'maximum string length exceeded'
          endif
          attchar1=' '
          attchar2=' '

          ierr = nf90_get_att(file1%fh,id1, trim(attname), attchar1)
          ierr = nf90_get_att(file2%fh,id2, trim(attname), attchar2)
          if(ierr==NF90_NOERR) then
             if(trim(attname).ne.'case' .and. attchar1(1:attlen) .ne. attchar2(1:attlen)) then
                print *, 'Attribute ',trim(attname),' from file1: ',attchar1(1:attlen),&
                     ' does not match that found on file2: ',attchar2(1:attlen)
             end if
          else
             print *, 'Attribute ',trim(attname),' from file1: ',attchar1(1:attlen),&
                  ' not found on file2'
          end if
          if(id1==NF90_GLOBAL .and. trim(attname) .eq. 'case') then
             print *, 'CASE 1 : ',trim(attchar1)
             print *, 'CASE 2 : ',trim(attchar2)
          endif
          if(id1==NF90_GLOBAL .and. trim(attname) .eq. 'title') then
             print *, 'TITLE 1 : ',trim(attchar1)
             print *, 'TITLE 2 : ',trim(attchar2)
          end if
       case(nf90_int)
          allocate(attint1(attlen),attint2(attlen))
          ierr = nf90_get_att(file1%fh,id1, trim(attname), attint1)
          ierr = nf90_get_att(file2%fh,id2, trim(attname), attint2)

          if(ierr==NF90_NOERR) then
             if(any(attint1 /= attint2)) then
                print *, 'Attribute ',trim(attname),' from file1: ',attint1,' does not match that found on file2 ',attint2
             end if
          else
             print *, 'Attribute ',trim(attname),' from file1: ',attint1,' not found on file2'
          end if
          deallocate(attint1, attint2)


       case(nf90_float)
          allocate(attreal1(attlen),attreal2(attlen))
          ierr = nf90_get_att(file1%fh,id1, trim(attname), attreal1)
          ierr = nf90_get_att(file2%fh,id2, trim(attname), attreal2)
          if(ierr==NF90_NOERR) then
             if(any(attreal1 /= attreal2)) then
                print *, 'Attribute ',trim(attname),' from file1: ',attreal1,' does not match that found on file2 ',attreal2
             end if
          else
             print *, 'Attribute ',trim(attname),' from file1: ',attreal1,' not found on file2'
          end if
          deallocate(attreal1, attreal2)
       case(nf90_double)
          allocate(attdouble1(attlen), attdouble2(attlen))
          ierr = nf90_get_att(file1%fh,id1, trim(attname), attdouble1)
          ierr = nf90_get_att(file2%fh,id2, trim(attname), attdouble2)
          if(ierr==NF90_NOERR) then
             if(any(attdouble1 /= attdouble2)) then
                print *, 'Attribute ',trim(attname),' from file1: ',attdouble1,' does not match that found on file2 ',attdouble2
             end if
          else
             print *, 'Attribute ',trim(attname),' from file1: ',attdouble1,' not found on file2'
          end if
          deallocate(attdouble1, attdouble2)
       case default
          print *,' Did not recognize attribute with id: ',i,' type: ',atttype, ' name: ',trim(attname), ' len: ',attlen
       end select
    end do

  end subroutine compare_metadata








  subroutine compare_dimensions( dimfile1, dimfile2)
    type(dim_t), intent(in) :: dimfile1(:), dimfile2(:)

    integer :: ds1, ds2
    integer :: i, j
    logical,pointer :: found(:,:)

    ds1 = size(dimfile1)
    ds2 = size(dimfile2)

    allocate(found(2,max(ds1,ds2)))

    found = .false.
    do i=1,ds1
       do j=1,ds2
          if(dimfile1(i)%name .eq. dimfile2(j)%name) then
             if(dimfile1(i)%dimsize == dimfile2(j)%dimsize) then
                print *, 'Dimension ',trim(dimfile1(i)%name), ' matches'
             else
                print *, 'Dimension ',trim(dimfile1(i)%name), ' differs ', dimfile1(i)%dimsize, ' /= ',dimfile2(j)%dimsize
             end if
             found(1,i) = .true.
             found(2,j) = .true.
          end if
       end do
    end do
    do i=1,ds1
       if(.not. found(1,i)) then
          print *, 'Could not find match for file 1 dimension ',trim(dimfile1(i)%name)
       end if
    end do
    do i=1,ds2
       if(.not. found(2,i)) then
          print *, 'Could not find match for file 2 dimension ',trim(dimfile2(i)%name)
       end if
    end do
    deallocate(found)
  end subroutine compare_dimensions


  subroutine match_vars( file1, file2, &
       num_not_found_on_file1, num_not_found_on_file2, &
       num_not_found_on_file1_timeconst, num_not_found_on_file2_timeconst)
    type(file_t), intent(inout) :: file1, file2

    ! Accumulates count of variables on file2 not found on file1; this only considers (a)
    ! fields with an unlimited (time) dimension, and (b) fields without an unlimited
    ! (time) dimension on a file that doesn't have an unlimited dimension.
    integer, intent(inout) :: num_not_found_on_file1

    ! Accumulates count of variables on file1 not found on file2; this only considers (a)
    ! fields with an unlimited (time) dimension, and (b) fields without an unlimited
    ! (time) dimension on a file that doesn't have an unlimited dimension.
    integer, intent(inout) :: num_not_found_on_file2

    ! Accumulates count of variables on file2 not found on file1; this only considers
    ! fields without an unlimited (time) dimension on a file that has an unlimited
    ! dimension.
    integer, intent(inout) :: num_not_found_on_file1_timeconst

    ! Accumulates count of variables on file1 not found on file2; this only considers
    ! fields without an unlimited (time) dimension on a file that has an unlimited
    ! dimension.
    integer, intent(inout) :: num_not_found_on_file2_timeconst

    type(var_t), pointer :: varfile1(:),varfile2(:)

    integer :: vs1, vs2, i, j



    varfile1 => file1%var
    varfile2 => file2%var

    vs1 = size(varfile1)
    vs2 = size(varfile2)

    do i=1,vs1
       do j=1,vs2
          if(varfile1(i)%name .eq. varfile2(j)%name) then
             varfile1(i)%matchid=j
             varfile2(j)%matchid=i
          end if
       end do
    end do
    do i=1,vs1
       if(varfile1(i)%matchid<0) then
          print *, 'Could not find match for file1 variable ',trim(varfile1(i)%name), ' in file2'
          if (file1%has_unlimited_dim() .and. &
               .not. is_time_varying(varfile1(i), file1%has_unlimited_dim(), file1%unlimdimid)) then
             num_not_found_on_file2_timeconst = num_not_found_on_file2_timeconst + 1
          else
             num_not_found_on_file2 = num_not_found_on_file2 + 1
          end if
       end if
    end do
    do i=1,vs2
       if(varfile2(i)%matchid<0) then
          print *, 'Could not find match for file2 variable ',trim(varfile2(i)%name), ' in file1'
          if (file2%has_unlimited_dim() .and. &
               .not. is_time_varying(varfile2(i), file2%has_unlimited_dim(), file2%unlimdimid)) then
             num_not_found_on_file1_timeconst = num_not_found_on_file1_timeconst + 1
          else
             num_not_found_on_file1 = num_not_found_on_file1 + 1
          end if
       end if
    end do
  end subroutine match_vars


  function is_time_varying(var, file_has_unlimited_dim, unlimdimid)
    type(var_t), intent(in) :: var                    ! variable of interest
    logical    , intent(in) :: file_has_unlimited_dim ! true if the file has an unlimited dimension
    integer    , intent(in) :: unlimdimid             ! the file's unlimited dim id (if it has one)

    logical :: is_time_varying     ! true if the given variable is time-varying

    if (file_has_unlimited_dim) then
       is_time_varying = any(var%dimids == unlimdimid)
    else
       is_time_varying = .false.
    end if
  end function is_time_varying


  function vdimsize(dims, dimids)
    type(dim_t), intent(in) :: dims(:)
    integer, intent(in) :: dimids(:)

    integer :: vdimsize
    integer :: i

    vdimsize=1
    do i=1,size(dimids)
       if(verbose) print *,__FILE__,__LINE__,i,dimids(i),size(dims),size(dimids)
       vdimsize = vdimsize*dims(dimids(i))%kount
    end do

  end function vdimsize





  subroutine compare_var_int(f1, f2, i1, i2, t)
    type(file_t) :: f1,f2
    integer, intent(in) :: i1, i2
    integer, optional :: t


    integer :: s1, s2, m1, m2, l1(1), l2(1), i, ierr
    integer, pointer :: v1(:), v2(:), vdiff(:)
    integer :: t1, n1
    integer :: start(NF90_MAX_DIMS), count(NF90_MAX_DIMS)

    if(present(t)) then
       t1 = t
    else
       t1 = 1
    end if

    s1 = vdimsize(f1%dim, f1%var(i1)%dimids)
    s2 = vdimsize(f2%dim, f2%var(i2)%dimids)

    if(s1 /= s2) then
       print *, 'Variable ',f1%var(i)%name,' sizes differ'
    end if

    n1 = size(f1%var(i1)%dimids)
    start = 1
    do i=1,n1
       count(i) = f1%dim(f1%var(i1)%dimids(i))%dimsize
       if(f1%var(i1)%dimids(i) == f1%unlimdimid) then
          count(i)=1
          start(i)=t1
       end if
    end do

    allocate(v1(s1), v2(s2))

    ierr = nf90_get_var(f1%fh, i1, v1, start(1:n1), count(1:n1))
    ierr = nf90_get_var(f2%fh, i2, v2, start(1:n1), count(1:n1))

    if(any(v1 /= v2)) then
       allocate(vdiff(s1))
       vdiff = abs(v1-v2)
       m1 = maxval(vdiff)
       m2 = minval(vdiff)
       l1 = maxloc(vdiff)
       l2 = minloc(vdiff)

       print *,__FILE__,__LINE__,m1,m2,l1,l2
       deallocate(vdiff)
    end if

    deallocate(v1,v2)
  end subroutine compare_var_int

  subroutine compare_var_float(f1, f2, i1, i2, t)
    type(file_t) :: f1,f2
    integer, intent(in) :: i1, i2
    integer, optional :: t


    integer :: s1, s2, m1, m2, l1(1), l2(1), i, ierr
    real, pointer :: v1(:), v2(:), vdiff(:)
    integer :: t1, n1
    integer :: start(NF90_MAX_DIMS), count(NF90_MAX_DIMS)

    if(present(t)) then
       t1 = t
    else
       t1 = 1
    end if

    s1 = vdimsize(f1%dim, f1%var(i1)%dimids)
    s2 = vdimsize(f2%dim, f2%var(i2)%dimids)

    if(s1 /= s2) then
       print *, 'Variable ',f1%var(i)%name,' sizes differ'
    end if

    n1 = size(f1%var(i1)%dimids)
    start = 1
    do i=1,n1
       count(i) =  f1%dim(f1%var(i1)%dimids(i))%dimsize
       if(f1%var(i1)%dimids(i) == f1%unlimdimid) then
          count(i)=1
          start(i)=t1
       end if
    end do

    allocate(v1(s1), v2(s2))

    ierr = nf90_get_var(f1%fh, i1, v1, start(1:n1), count(1:n1))
    ierr = nf90_get_var(f2%fh, i2, v2, start(1:n1), count(1:n1))

    if(any(v1 /= v2)) then
       allocate(vdiff(s1))
       vdiff = abs(v1-v2)
       m1 = maxval(vdiff)
       m2 = minval(vdiff)
       l1 = maxloc(vdiff)
       l2 = minloc(vdiff)

       print *,__FILE__,__LINE__,m1,m2,l1,l2
       deallocate(vdiff)
    end if

    deallocate(v1,v2)
  end subroutine compare_var_float

  subroutine compare_var_double(f1, f2, i1, i2, t)
    type(file_t) :: f1,f2
    integer, intent(in) :: i1, i2
    integer, optional :: t


    integer :: s1, s2, m1, m2, l1(1), l2(1), i, ierr
    double precision, pointer :: v1(:), v2(:), vdiff(:)
    integer :: t1, n1
    integer :: start(NF90_MAX_DIMS), count(NF90_MAX_DIMS)

    if(present(t)) then
       t1 = t
    else
       t1 = 1
    end if

    s1 = vdimsize(f1%dim, f1%var(i1)%dimids)
    s2 = vdimsize(f2%dim, f2%var(i2)%dimids)

    if(s1 /= s2) then
       print *, 'Variable ',f1%var(i)%name,' sizes differ'
    end if

    n1 = size(f1%var(i1)%dimids)
    start = 1
    do i=1,n1
       count(i) =  f1%dim(f1%var(i1)%dimids(i))%dimsize
       if(f1%var(i1)%dimids(i) == f1%unlimdimid) then
          count(i)=1
          start(i)=t1
       end if
    end do

    allocate(v1(s1), v2(s2))

    ierr = nf90_get_var(f1%fh, i1, v1, start(1:n1), count(1:n1))
    ierr = nf90_get_var(f2%fh, i2, v2, start(1:n1), count(1:n1))

    if(any(v1 /= v2)) then
       allocate(vdiff(s1))
       vdiff = abs(v1-v2)
       m1 = maxval(vdiff)
       m2 = minval(vdiff)
       l1 = maxloc(vdiff)
       l2 = minloc(vdiff)

       print *,__FILE__,__LINE__,m1,m2,l1,l2
       deallocate(vdiff)
    end if

    deallocate(v1,v2)
  end subroutine compare_var_double









end module filestruct