Changeset 4489 for LMDZ6/trunk/libf/phylmd/ecrad/easy_netcdf.F90

Timestamp:

Mar 31, 2023, 8:42:57 PM (15 months ago)

Author:

lguez

Message:

Merge LMDZ_ECRad branch back into trunk!

Location:

LMDZ6/trunk

Files:

• . (modified) (1 prop)
• libf/phylmd/ecrad/easy_netcdf.F90 (modified) (44 diffs)

LMDZ6/trunk/libf/phylmd/ecrad/easy_netcdf.F90

@@ -35 +35 @@
   ! a NetCDF file
   type netcdf_file
-    integer :: ncid             ! NetCDF file ID
+    integer :: ncid = -1! NetCDF file ID
     integer :: iverbose ! Verbosity: 0 = report only fatal errors,
                         !            1 = ...and warnings,
@@ -55 +55 @@
     procedure :: create => create_netcdf_file
     procedure :: close => close_netcdf_file
+    procedure :: is_open
     procedure :: get_real_scalar
+    procedure :: get_int_scalar
     procedure :: get_real_vector
-    procedure :: get_integer_vector
+    procedure :: get_int_vector
     procedure :: get_real_matrix
     procedure :: get_real_array3
@@ -65 +67 @@
     procedure :: get_real_array3_indexed
     procedure :: get_real_array4
-    generic   :: get => get_real_scalar, get_real_vector, &
+    procedure :: get_char_vector
+    procedure :: get_char_matrix
+    generic   :: get => get_real_scalar, get_int_scalar, &
+         &              get_real_vector, get_int_vector, &
          &              get_real_matrix, get_real_array3, &
-         &              get_real_array4, get_integer_vector, &
+         &              get_real_array4, &
          &              get_real_scalar_indexed, get_real_vector_indexed, &
-         &              get_real_matrix_indexed, get_real_array3_indexed
+         &              get_real_matrix_indexed, get_real_array3_indexed, &
+         &              get_char_vector, get_char_matrix
     procedure :: get_real_scalar_attribute
     procedure :: get_string_attribute
@@ -83 +89 @@
     procedure :: put_real_scalar
     procedure :: put_real_vector
+    procedure :: put_int_vector
     procedure :: put_real_matrix
     procedure :: put_real_array3
@@ -91 +98 @@
          &              put_real_matrix, put_real_array3, &
          &              put_real_scalar_indexed, put_real_vector_indexed, &
-         &              put_real_matrix_indexed
+         &              put_real_matrix_indexed, put_int_vector
     procedure :: set_verbose
     procedure :: transpose_matrices
@@ -101 +108 @@
     procedure :: attribute_exists
     procedure :: global_attribute_exists
+#ifdef NC_NETCDF4
+    procedure :: copy_dimensions
+#endif
+    procedure :: copy_variable_definition
+    procedure :: copy_variable
     procedure, private :: get_array_dimensions
     procedure, private :: get_variable_id
@@ -257 +269 @@
     end if
 
+    this%ncid = -1
+
   end subroutine close_netcdf_file
 
@@ -367 +381 @@
 
     integer                        :: j, istatus
-    integer                        :: dimids(NF90_MAX_VAR_DIMS)
+    integer                        :: idimids(NF90_MAX_VAR_DIMS)
 
     istatus = nf90_inquire_variable(this%ncid, ivarid, &
-         &                          ndims=ndims, dimids=dimids)
+         &                          ndims=ndims, dimids=idimids)
     if (istatus /= NF90_NOERR) then
       write(nulerr,'(a,i0,a,a)') '*** Error inquiring about NetCDF variable with id ', &
@@ -379 +393 @@
     ndimlens(:) = 0
     do j = 1,ndims
-      istatus = nf90_inquire_dimension(this%ncid, dimids(j), len=ndimlens(j))
+      istatus = nf90_inquire_dimension(this%ncid, idimids(j), len=ndimlens(j))
       if (istatus /= NF90_NOERR) then
         write(nulerr,'(a,i0,a,i0,a,a)') '*** Error reading length of dimension ', &
@@ -420 +434 @@
 
   !---------------------------------------------------------------------
+  ! Return true if file is open, false otherwise
+  function is_open(this)
+    class(netcdf_file) :: this
+    logical            :: is_open
+    is_open = (this%ncid >= 0)
+  end function is_open
+
+  !---------------------------------------------------------------------
   ! Return the number of dimensions of variable with name var_name, or
   ! -1 if the variable is not found
@@ -619 +641 @@
 
   !---------------------------------------------------------------------
+  ! Read an integer scalar
+  subroutine get_int_scalar(this, var_name, scalar)
+    class(netcdf_file)           :: this
+    character(len=*), intent(in) :: var_name
+    integer,          intent(out):: scalar
+
+    integer                      :: istatus
+    integer                      :: ivarid, ndims
+    integer                      :: ndimlens(NF90_MAX_VAR_DIMS)
+    integer                      :: j, ntotal
+
+    ! Inquire the ID, shape & size of the variable
+    call this%get_variable_id(var_name, ivarid)
+    call this%get_array_dimensions(ivarid, ndims, ndimlens)
+
+    ! Compute number of elements of the variable in the file
+    ntotal = 1
+    do j = 1, ndims
+      ntotal = ntotal * ndimlens(j)
+    end do
+
+    if (this%iverbose >= 3) then
+      write(nulout,'(a,a)',advance='no') '  Reading ', var_name
+      call this%print_variable_attributes(ivarid,nulout)
+    end if
+
+    ! Abort if the number of elements is anything other than 1
+    if (ntotal /= 1) then
+      write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
+           &    var_name, ' with total length ', ntotal, ' as a scalar'
+      call my_abort('Error reading NetCDF file')
+    end if
+
+    ! Read variable
+    istatus = nf90_get_var(this%ncid, ivarid, scalar)
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
+           &  var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+  end subroutine get_int_scalar
+
+
+  !---------------------------------------------------------------------
   ! Read a scalar from a larger array, where "index" indexes the most
   ! slowly varying dimension
-  subroutine get_real_scalar_indexed(this, var_name, index, scalar)
+  subroutine get_real_scalar_indexed(this, var_name, scalar, index)
     class(netcdf_file)           :: this
     character(len=*), intent(in) :: var_name
@@ -676 +743 @@
 
   !---------------------------------------------------------------------
-  ! Read a 1D array into "vector", which must be allocatable and will
-  ! be reallocated if necessary
+  ! Read a 1D real array into "vector", which must be allocatable and
+  ! will be reallocated if necessary
   subroutine get_real_vector(this, var_name, vector)
     class(netcdf_file)           :: this
@@ -734 +801 @@
 
   !---------------------------------------------------------------------
-  ! Read a 1D integer array into "vector", which must be allocatable
+  ! Read a 1D character array into "vector", which must be allocatable
   ! and will be reallocated if necessary
-  subroutine get_integer_vector(this, var_name, vector)
+  subroutine get_char_vector(this, var_name, vector)
     class(netcdf_file)           :: this
     character(len=*), intent(in) :: var_name
-    integer, allocatable, intent(out) :: vector(:)
+    character(len=1), allocatable, intent(out) :: vector(:)
 
     integer                      :: n  ! Length of vector
@@ -784 +851 @@
     if (istatus /= NF90_NOERR) then
       write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
+           &  var_name, ' as a vector of chars: ', trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+  end subroutine get_char_vector
+
+
+  !---------------------------------------------------------------------
+  ! Read a 1D integer array into "vector", which must be allocatable
+  ! and will be reallocated if necessary
+  subroutine get_int_vector(this, var_name, vector)
+
+    class(netcdf_file)           :: this
+    character(len=*), intent(in) :: var_name
+    integer, allocatable, intent(out) :: vector(:)
+
+    integer                      :: n  ! Length of vector
+    integer                      :: istatus
+    integer                      :: ivarid, ndims
+    integer                      :: ndimlens(NF90_MAX_VAR_DIMS)
+    integer                      :: j
+
+    call this%get_variable_id(var_name, ivarid)
+    call this%get_array_dimensions(ivarid, ndims, ndimlens)
+
+    ! Ensure variable has only one dimension in the file
+    n = 1
+    do j = 1, ndims
+      n = n * ndimlens(j)
+      if (j > 1 .and. ndimlens(j) > 1) then
+        write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
+             & var_name, &
+             & ' as a vector: all dimensions above the first must be singletons'
+        call my_abort('Error reading NetCDF file')
+      end if
+    end do
+
+    ! Reallocate if necessary
+    if (allocated(vector)) then
+      if (size(vector) /= n) then
+        if (this%iverbose >= 1) then
+          write(nulout,'(a,a)') '  Warning: resizing vector to read ', var_name
+        end if
+        deallocate(vector)
+        allocate(vector(n))
+      end if
+    else
+      allocate(vector(n))
+    end if
+
+    if (this%iverbose >= 3) then
+      write(nulout,'(a,a,a,i0,a)',advance='no') '  Reading ', var_name, '(', n, ')'
+      call this%print_variable_attributes(ivarid,nulout)
+    end if
+
+    ! Read variable
+    istatus = nf90_get_var(this%ncid, ivarid, vector)
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
            &  var_name, ' as an integer vector: ', trim(nf90_strerror(istatus))
       call my_abort('Error reading NetCDF file')
     end if
 
-  end subroutine get_integer_vector
-
+  end subroutine get_int_vector
 
   !---------------------------------------------------------------------
   ! Read a vector of data from a larger array; the vector must be
   ! allocatable and will be reallocated if necessary
-  subroutine get_real_vector_indexed(this, var_name, index, vector)
+  subroutine get_real_vector_indexed(this, var_name, vector, index)
     class(netcdf_file)           :: this
     character(len=*), intent(in) :: var_name
@@ -975 +1100 @@
 
   !---------------------------------------------------------------------
+  ! Read 2D array of characters into "matrix", which must be
+  ! allocatable and will be reallocated if necessary.  Whether to
+  ! transpose is specifed by the final optional argument, but can also
+  ! be specified by the do_transpose_2d class data member.
+  subroutine get_char_matrix(this, var_name, matrix, do_transp)
+    class(netcdf_file)           :: this
+    character(len=*), intent(in) :: var_name
+    character(len=1), allocatable, intent(inout) :: matrix(:,:)
+    logical, optional, intent(in):: do_transp ! Transpose data?
+
+    character(len=1), allocatable:: tmp_matrix(:,:)
+    integer                      :: ndimlen1, ndimlen2
+    integer                      :: istatus
+    integer                      :: ivarid, ndims
+    integer                      :: ndimlens(NF90_MAX_VAR_DIMS)
+    integer                      :: vstart(NF90_MAX_VAR_DIMS)
+    integer                      :: vcount(NF90_MAX_VAR_DIMS)
+    integer                      :: j, ntotal
+    logical                      :: do_transpose
+
+    ! Decide whether to transpose the array
+    if (present(do_transp)) then
+      do_transpose = do_transp
+    else
+      do_transpose = this%do_transpose_2d
+    end if
+
+    call this%get_variable_id(var_name, ivarid)
+    call this%get_array_dimensions(ivarid, ndims, ndimlens)
+
+    ! Ensure the variable has no more than two non-singleton
+    ! dimensions
+    ntotal = 1
+    do j = 1, ndims
+      ntotal = ntotal * ndimlens(j)
+      if (j > 2 .and. ndimlens(j) > 1) then
+        write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
+           & var_name, &
+           & ' as a matrix: all dimensions above the second must be singletons'
+        call my_abort('Error reading NetCDF file')
+      end if
+    end do
+
+    ! Work out dimension lengths
+    if (ndims >= 2) then
+      ndimlen1 = ndimlens(1)
+      ndimlen2 = ntotal/ndimlen1
+    else
+      ndimlen1 = ntotal
+      ndimlen2 = 1
+    end if
+
+    if (do_transpose) then
+      ! Read and transpose
+      allocate(tmp_matrix(ndimlen1, ndimlen2))
+
+      ! Reallocate if necessary
+      if (allocated(matrix)) then
+        if (size(matrix,1) /= ndimlen2 .or. size(matrix,2) /= ndimlen1) then
+          if (this%iverbose >= 1) then
+            write(nulout,'(a,a)') '  Warning: resizing matrix to read ', var_name
+          end if
+          deallocate(matrix)
+          allocate(matrix(ndimlen2, ndimlen1))
+        end if
+      else
+        allocate(matrix(ndimlen2, ndimlen1))
+      end if
+
+      if (this%iverbose >= 3) then
+        write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') '  Reading ', var_name, '(', &
+             &                            ndimlen2, ',', ndimlen1, ')'
+        call this%print_variable_attributes(ivarid,nulout)
+      end if
+
+      istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix)
+      matrix = transpose(tmp_matrix)
+      deallocate(tmp_matrix)
+    else
+      ! Read data without transposition
+
+      ! Reallocate if necessary
+      if (allocated(matrix)) then
+        if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
+          if (this%iverbose >= 1) then
+            write(nulout,'(a,a)') '  Warning: resizing matrix to read ', var_name
+          end if
+          allocate(matrix(ndimlen1, ndimlen2))
+        end if
+      else
+        allocate(matrix(ndimlen1, ndimlen2))
+      end if
+
+      if (this%iverbose >= 3) then
+        write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') '  Reading ', var_name, '(', &
+             &                            ndimlen1, ',', ndimlen2, ')'
+        call this%print_variable_attributes(ivarid,nulout)
+      end if
+
+      vstart = 1
+      vcount(1:2) = [ndimlen1,1]
+      
+      do j = 1,ndimlen2
+        vstart(2) = j
+        istatus = nf90_get_var(this%ncid, ivarid, matrix(:,j), start=vstart, count=vcount)
+      end do
+    end if
+
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
+           &    var_name, ' as a matrix of characters: ', trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+  end subroutine get_char_matrix
+
+
+  !---------------------------------------------------------------------
   ! Read matrix of data from a larger array, which must be allocatable
   ! and will be reallocated if necessary.  Whether to transpose is
   ! specifed by the final optional argument, but can also be specified
   ! by the do_transpose_2d class data member.
-  subroutine get_real_matrix_indexed(this, var_name, index, matrix, do_transp)
+  subroutine get_real_matrix_indexed(this, var_name, matrix, index, do_transp)
     class(netcdf_file)           :: this
     character(len=*), intent(in) :: var_name
@@ -1187 +1430 @@
       if (this%iverbose >= 3) then
         write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') '  Reading ', var_name, &
-             & ' (permuted dimensions ', i_permute_3d, ')'
+             & ' (permuting dimensions ', i_permute_3d, ')'
         call this%print_variable_attributes(ivarid,nulout)
       end if
@@ -1235 +1478 @@
   ! be allocatable and will be reallocated if necessary.  Whether to
   ! pemute is specifed by the final optional argument
-  subroutine get_real_array3_indexed(this, var_name, index, var, ipermute)
+  subroutine get_real_array3_indexed(this, var_name, var, index, ipermute)
     class(netcdf_file)                   :: this
     character(len=*), intent(in)         :: var_name
@@ -1331 +1574 @@
         write(nulout,'(a,i0,a,a,a,i0,i0,i0,a)') '  Reading slice ', index, &
              &  ' of ', var_name, &
-             & ' (permuted dimensions ', i_permute_3d, ')'
+             & ' (permuting dimensions ', i_permute_3d, ')'
       end if
 
@@ -1471 +1714 @@
       if (this%iverbose >= 3) then
         write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') '  Reading ', var_name, &
-             & ' (permuted dimensions ', i_permute_4d, ')'
+             & ' (permuting dimensions ', i_permute_4d, ')'
         call this%print_variable_attributes(ivarid,nulout)
       end if
@@ -1549 +1792 @@
     !    allocate(character(len=i_attr_len) :: attr_str)
     if (len(attr_str) < i_attr_len) then
-      write(nulerr,'(a,a)') '*** Error: not enough space to read attribute ', attr_name
+      write(nulerr,'(a,a)') '*** Not enough space to read attribute ', attr_name
       call my_abort('Error reading NetCDF file')
     end if
@@ -1577 +1820 @@
     real(jprb),       intent(out) :: attr
 
-    integer :: i_attr_len, ivarid
+    integer :: ivarid
     integer :: istatus
-    integer :: j
 
     istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
@@ -1624 +1866 @@
     !    allocate(character(len=i_attr_len) :: attr_str)
     if (len(attr_str) < i_attr_len) then
-      write(nulerr,'(a,a)') '*** Error: not enough space to read global attribute ', attr_name
+      write(nulerr,'(a,a)') '*** Not enough space to read global attribute ', attr_name
       call my_abort('Error reading NetCDF file')
     end if
@@ -1652 +1894 @@
 
     character(len=4000) :: attr_str
-    integer :: i_attr_len
     integer :: istatus
-    integer :: j
 
     if (this%iverbose >= 4) then
       istatus = nf90_get_att(this%ncid, ivarid, 'long_name', attr_str)
       if (istatus == NF90_NOERR) then
-        write(iunit, '(a)') ':'
-        write(iunit, '(a,a)', advance='no') '    ', trim(attr_str)
+        write(iunit, '(a,a,a)', advance='no') ': "', trim(attr_str), '"'
         istatus = nf90_get_att(this%ncid, ivarid, 'units', attr_str)
         if (istatus == NF90_NOERR) then
@@ -1721 +1960 @@
   ! names.
   subroutine define_variable(this, var_name, dim1_name, dim2_name, dim3_name, &
-       &                     long_name, units_str, comment_str, standard_name, is_double, &
-       &                     data_type_name, fill_value, deflate_level, shuffle, chunksizes)
+       &                     dim4_name, long_name, units_str, comment_str, &
+       &                     standard_name, is_double, data_type_name, fill_value, &
+       &                     deflate_level, shuffle, chunksizes, ndims)
     class(netcdf_file)                     :: this
     character(len=*), intent(in)           :: var_name
     character(len=*), intent(in), optional :: long_name, units_str, comment_str, standard_name
-    character(len=*), intent(in), optional :: dim1_name, dim2_name, dim3_name
+    character(len=*), intent(in), optional :: dim1_name, dim2_name, dim3_name, dim4_name
     logical,          intent(in), optional :: is_double
     character(len=*), intent(in), optional :: data_type_name
@@ -1733 +1973 @@
     logical,          intent(in), optional :: shuffle ! Shuffle bytes before compression
     integer, dimension(:), intent(in), optional :: chunksizes
-
-    integer :: istatus, ndims, ivarid
+    integer,          intent(in), optional :: ndims
+
+    integer :: istatus, ndims_local, ndims_input, ivarid
     integer, dimension(NF90_MAX_VAR_DIMS) :: idimids
     integer :: data_type
 
-    if (present(dim1_name)) then
+    ! Sometimes a program may not know at compile time the exact
+    ! dimensions of a variable - if ndims is present then only up to
+    ! that many dimensions will be defined
+    ndims_input = 4
+    if (present(ndims)) then
+      ndims_input = ndims
+    end if
+
+    if (present(dim1_name) .and. ndims_input >= 1) then
       ! Variable is at least one dimensional
-      ndims = 1
+      ndims_local = 1
       istatus = nf90_inq_dimid(this%ncid, dim1_name, idimids(1))
       if (istatus /= NF90_NOERR) then
@@ -1747 +1996 @@
         call my_abort('Error writing NetCDF file')
       end if
-      if (present(dim2_name)) then
+      if (present(dim2_name) .and. ndims_input >= 2) then
         ! Variable is at least two dimensional
-        ndims = 2
+        ndims_local = 2
         istatus = nf90_inq_dimid(this%ncid, dim2_name, idimids(2))
         if (istatus /= NF90_NOERR) then
@@ -1756 +2005 @@
           call my_abort('Error writing NetCDF file')
         end if
-        if (present(dim3_name)) then
+        if (present(dim3_name) .and. ndims_input >= 3) then
           ! Variable is at least three dimensional
-          ndims = 3
+          ndims_local = 3
           istatus = nf90_inq_dimid(this%ncid, dim3_name, idimids(3))
           if (istatus /= NF90_NOERR) then
@@ -1765 +2014 @@
             call my_abort('Error writing NetCDF file')
           end if
+          if (present(dim4_name) .and. ndims_input >= 4) then
+            ! Variable is at least three dimensional
+            ndims_local = 4
+            istatus = nf90_inq_dimid(this%ncid, dim4_name, idimids(4))
+            if (istatus /= NF90_NOERR) then
+              write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
+                   &             dim4_name, ': ', trim(nf90_strerror(istatus))
+              call my_abort('Error writing NetCDF file')
+            end if
+          end if
         end if
       end if
     else
       ! Variable is a scalar
-      ndims = 0
+      ndims_local = 0
     end if
 
     ! Read output precision from optional argument "is_double" if
     ! present, otherwise from default output precision for this file
-    data_type = NF90_FLOAT ! Default
     if (present(data_type_name)) then
       if (data_type_name == 'double') then
@@ -1787 +2045 @@
         data_type = NF90_FLOAT
       else
-        write(nulerr,'(a,a,a)') '*** Error: netCDF data type "', data_type_name, '" not supported'
+        write(nulerr,'(a,a,a)') '*** NetCDF data type "', data_type_name, '" not supported'
         call my_abort('Error writing NetCDF file')
       end if
     else if (present(is_double)) then
+      if (is_double) then
+        data_type = NF90_DOUBLE
+      else
+        data_type = NF90_FLOAT
+      end if
+    else if (this%is_double_precision) then
       data_type = NF90_DOUBLE
+    else
+      data_type = NF90_FLOAT
     end if
 
     ! Define variable
 #ifdef NC_NETCDF4
-    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims), &
+    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims_local), &
          & ivarid, deflate_level=deflate_level, shuffle=shuffle, chunksizes=chunksizes)
 #else
-    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims), ivarid)
+    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims_local), ivarid)
 #endif
     if (istatus /= NF90_NOERR) then
@@ -1811 +2077 @@
       istatus = nf90_put_att(this%ncid, ivarid, "long_name", long_name)
       if (this%iverbose >= 4) then
-        write(nulout,'(a,a,a,a)') '  Defining ',trim(var_name),': ',long_name
+        write(nulout,'(a,a,a,a,a)', advance='no') '  Defining ',trim(var_name), &
+             &  ': "', long_name, '"'
       end if
     else
       if (this%iverbose >= 4) then
-        write(nulout,'(a,a)') '  Defining ',trim(var_name)
-      end if
-    end if
+        write(nulout,'(a,a)', advance='no') '  Defining ',trim(var_name)
+      end if
+    end if
+
     if (present(units_str)) then
       istatus = nf90_put_att(this%ncid, ivarid, "units", units_str)
-    end if
+      if (this%iverbose >= 4) then
+        if (trim(units_str) == '1') then
+          write(nulout, '(a)') ' (dimensionless)'
+        else
+          write(nulout, '(a,a,a)') ' (', trim(units_str), ')'
+        end if
+      end if
+    else
+      if (this%iverbose >= 4) then
+        write(nulout, '(1x)')
+      end if
+    end if
+
     if (present(standard_name)) then
       istatus = nf90_put_att(this%ncid, ivarid, "standard_name", standard_name)
@@ -1863 +2143 @@
   subroutine put_global_attributes(this, title_str, inst_str, source_str, &
        &  comment_str, references_str, creator_name, creator_email_str, &
-       &  contributor_name, project_str, conventions_str)
+       &  contributor_name, project_str, conventions_str, prior_history_str)
     class(netcdf_file)                     :: this
 
@@ -1872 +2152 @@
     character(len=*), intent(in), optional :: contributor_name, project_str
     character(len=*), intent(in), optional :: comment_str, conventions_str
-    character(len=*), intent(in), optional :: references_str
+    character(len=*), intent(in), optional :: references_str, prior_history_str
 
     character(len=32)   :: date_time_str
@@ -1887 +2167 @@
          &   time_vals(1), time_vals(2), time_vals(3), time_vals(5), time_vals(6), time_vals(7)
 
-    history_str = trim(date_time_str) // ': ' // trim(command_line_str)
+    if (present(prior_history_str)) then
+      history_str = trim(prior_history_str) // new_line('a') &
+           &  // trim(date_time_str) // ': ' // trim(command_line_str)
+    else
+      history_str = trim(date_time_str) // ': ' // trim(command_line_str)
+    end if
 
     if (present(title_str))   i=nf90_put_att(this%ncid, NF90_GLOBAL, "title", title_str)
@@ -2057 +2342 @@
 
   !---------------------------------------------------------------------
+  ! Save an integer vector with name var_name in the file
+  subroutine put_int_vector(this, var_name, var)
+    class(netcdf_file)             :: this
+    character(len=*), intent(in)   :: var_name
+    integer,          intent(in)   :: var(:)
+
+    integer :: ivarid, ndims, istatus
+    integer(kind=jpib) :: ntotal
+    integer :: ndimlens(NF90_MAX_VAR_DIMS)
+
+    call this%end_define_mode()
+
+    ! Check the vector is of the right length
+    call this%get_variable_id(var_name, ivarid)
+    call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
+    if (ntotal /= size(var,kind=jpib)) then
+      write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
+           & size(var), ' to ', var_name, ' which has total length ', ntotal
+      call my_abort('Error writing NetCDF file')
+    end if
+
+    ! Save the vector
+    istatus = nf90_put_var(this%ncid, ivarid, var)
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error writing vector ', var_name, ': ', &
+           &                    trim(nf90_strerror(istatus))
+      call my_abort('Error writing NetCDF file')
+    end if
+
+  end subroutine put_int_vector
+
+
+  !---------------------------------------------------------------------
   ! Save a vector slice with name var_name in the file
-  subroutine put_real_vector_indexed(this, var_name, index, var)
+  subroutine put_real_vector_indexed(this, var_name, var, index2, index3)
     class(netcdf_file)             :: this
     character(len=*), intent(in)   :: var_name
     real(jprb), intent(in)         :: var(:)
-    integer, intent(in)            :: index
+    integer, intent(in)            :: index2
+    integer, intent(in), optional  :: index3
 
     integer :: ivarid, ndims, istatus
@@ -2070 +2389 @@
     integer :: vcount(NF90_MAX_VAR_DIMS)
 
+    character(len=512) :: var_slice_name
+    integer :: index_last
+
     call this%end_define_mode()
 
@@ -2076 +2398 @@
     call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
     ntotal = ntotal / ndimlens(ndims)
+    if (present(index3)) then
+      ntotal = ntotal / ndimlens(ndims-1)
+      index_last = index3
+      write(var_slice_name,'(a,a,i0,a,i0,a)') var_name, '(:,', index2, ',', index3, ')'
+    else
+      index_last = index2
+      write(var_slice_name,'(a,a,i0,a)') var_name, '(:,', index2, ')'
+    end if
+
     if (ntotal /= size(var,kind=jpib)) then
-      write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
-           & size(var), ' to slice of ', var_name, ' which has length ', ntotal
+      write(nulerr,'(a,i0,a,a,i0)') '*** Error: attempt to write vector of length ', &
+           & size(var), ' to ', trim(var_slice_name), ' which has length ', ntotal
       call my_abort('Error writing NetCDF file')
     end if
-    if (index < 1 .or. index > ndimlens(ndims)) then
-      write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector to slice ', &
-           &  index, ' of ', var_name, ' which has outer dimension  ', ndimlens(ndims)
+    if (index_last < 1 .or. index_last > ndimlens(ndims)) then
+      write(nulerr,'(a,a,a,i0)') '*** Error: attempt to write vector to ', &
+           &  trim(var_slice_name), ' which has outer dimension  ', ndimlens(ndims)
       call my_abort('Error writing NetCDF file')
     end if
@@ -2089 +2420 @@
     ! Save the vector
     vstart(1:ndims-1) = 1
-    vstart(ndims)     = index
     vcount(1:ndims-1) = ndimlens(1:ndims-1)
     vcount(ndims)     = 1
+    if (present(index3)) then
+      vstart(ndims)   = index3
+      vstart(ndims-1) = index2
+      vcount(ndims-1) = 1
+    else
+      vstart(ndims)   = index2
+    end if
+
     istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
     if (istatus /= NF90_NOERR) then
-      write(nulerr,'(a,a,a,a)') '*** Error writing vector to ', var_name, ': ', &
-           &                    trim(nf90_strerror(istatus))
+      write(nulerr,'(a,a,a,a)') '*** Error writing vector to ', trim(var_slice_name), &
+           &  ': ', trim(nf90_strerror(istatus))
       call my_abort('Error writing NetCDF file')
     end if
@@ -2171 +2509 @@
   ! dimensions if either optional argument transp is .true., or the
   ! transpose_matrices method has already been called.
-  subroutine put_real_matrix_indexed(this, var_name, index, var, do_transp)
+  subroutine put_real_matrix_indexed(this, var_name, var, index3, index4, do_transp)
     class(netcdf_file)             :: this
     character(len=*), intent(in)   :: var_name
     real(jprb), intent(in)         :: var(:,:)
-    integer, intent(in)            :: index
+    integer, intent(in)            :: index3
+    integer, intent(in), optional  :: index4
 
     real(jprb), allocatable        :: var_transpose(:,:)
-    logical, optional, intent(in):: do_transp
+    logical, optional, intent(in)  :: do_transp
 
     integer :: ivarid, ndims, nvarlen, istatus
@@ -2186 +2525 @@
     integer :: vcount(NF90_MAX_VAR_DIMS)
 
+    character(len=512) :: var_slice_name
+
     logical :: do_transpose
 
@@ -2204 +2545 @@
     ! ntotal is zero then there must be an unlimited dimension)
     ntotal = ntotal / ndimlens(ndims)
+    if (present(index4)) then
+      ntotal = ntotal / ndimlens(ndims-1)
+      write(var_slice_name,'(a,a,i0,a,i0,a)') var_name, '(:,:,', index3, ',', index4, ')'
+    else
+      write(var_slice_name,'(a,a,i0,a)') var_name, '(:,:,', index3, ')'
+    end if
     if (ntotal /= size(var,kind=jpib) .and. ntotal /= 0) then
       write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write matrix of total size ', &
-           & nvarlen, ' to ', var_name, ' which has total size ', ntotal
+           & nvarlen, ' to ', trim(var_slice_name), ' which has total size ', ntotal
       call my_abort('Error writing NetCDF file')
     end if
 
     vstart(1:ndims-1) = 1
-    vstart(ndims)     = index
     vcount(1:ndims-1) = ndimlens(1:ndims-1)
     vcount(ndims)     = 1
+    if (present(index4)) then
+      vstart(ndims)   = index4
+      vstart(ndims-1) = index3
+      vcount(ndims-1) = 1
+    else
+      vstart(ndims)   = index3
+    end if
 
     if (do_transpose) then
       ! Save the matrix with transposition
       if (this%iverbose >= 3) then
-        write(nulout,'(a,i0,a,a,a)') '  Writing slice ', index, ' of ', var_name, &
+        write(nulout,'(a,a,a)') '  Writing ', trim(var_slice_name), &
              & ' (transposing dimensions)'
       end if
@@ -2228 +2581 @@
       ! Save the matrix without transposition
       if (this%iverbose >= 3) then
-        write(nulout,'(a,i0,a,a)') '  Writing slice ', index, ' of ', var_name
+        write(nulout,'(a,a)') '  Writing ', trim(var_slice_name)
       end if
       istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
@@ -2234 +2587 @@
 
     if (istatus /= NF90_NOERR) then
-      write(nulerr,'(a,a,a,a)') '*** Error writing matrix ', var_name, &
+      write(nulerr,'(a,a,a)') '*** Error writing ', trim(var_slice_name), &
            &                    ': ', trim(nf90_strerror(istatus))
       call my_abort('Error writing NetCDF file')
@@ -2291 +2644 @@
       if (this%iverbose >= 3) then
         write(nulout,'(a,a,a,i0,i0,i0,a)') '  Writing ', var_name, &
-             & ' (permuted dimensions: ', i_permute_3d, ')'
+             & ' (permuting dimensions: ', i_permute_3d, ')'
       end if
       n_dimlens_permuted = (/ size(var,i_permute_3d(1)), &
            &                  size(var,i_permute_3d(2)), &
            &                  size(var,i_permute_3d(3))  /)
-      if (this%iverbose >= 4) then
-        write(nulout,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') '    (', &
-             &  n_dimlens_permuted(1), ',', n_dimlens_permuted(2), &
-             &  ',', n_dimlens_permuted(3), ') -> (', ndimlens(1), &
-             &  ',', ndimlens(2), ',', ndimlens(3), ')'
-      end if
+      !! FIX: This makes it look like the dimensions have stayed the same
+      ! if (this%iverbose >= 4) then
+      !   write(nulout,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') '    (', &
+      !        &  n_dimlens_permuted(1), ',', n_dimlens_permuted(2), &
+      !        &  ',', n_dimlens_permuted(3), ') -> (', ndimlens(1), &
+      !        &  ',', ndimlens(2), ',', ndimlens(3), ')'
+      ! end if
       allocate(var_permute(n_dimlens_permuted(1), &
            &   n_dimlens_permuted(2), n_dimlens_permuted(3)))
@@ -2326 +2680 @@
   end subroutine put_real_array3
 
+
+#ifdef NC_NETCDF4
+  !---------------------------------------------------------------------
+  ! Copy dimensions from "infile" to "this"
+  subroutine copy_dimensions(this, infile)
+    class(netcdf_file)            :: this
+    type(netcdf_file), intent(in) :: infile
+
+    integer :: jdim
+    integer :: ndims
+    integer :: idimids(1024)
+    integer :: dimlen
+    character(len=512) :: dimname
+    integer :: istatus
+    integer :: include_parents
+    
+    include_parents = 0
+
+    istatus = nf90_inq_dimids(infile%ncid, ndims, idimids, include_parents)
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a)') '*** Error reading dimensions of NetCDF file: ', &
+           trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+    do jdim = 1,ndims
+      istatus = nf90_inquire_dimension(infile%ncid, idimids(jdim), &
+           &  name=dimname, len=dimlen)
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a)') '*** Error reading NetCDF dimension properties: ', &
+             trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+      call this%define_dimension(trim(dimname), dimlen)
+    end do
+
+  end subroutine copy_dimensions
+#endif
+
+  !---------------------------------------------------------------------
+  ! Copy variable definition and attributes from "infile" to "this"
+  subroutine copy_variable_definition(this, infile, var_name)
+    class(netcdf_file)            :: this
+    type(netcdf_file), intent(in) :: infile
+    character(len=*),  intent(in) :: var_name
+
+#ifdef NC_NETCDF4
+    integer :: deflate_level  ! Compression: 0 (none) to 9 (most)
+    logical :: shuffle        ! Shuffle bytes before compression
+    integer :: chunksizes(NF90_MAX_VAR_DIMS)
+#endif
+    integer :: data_type
+    integer :: ndims
+    integer :: idimids_in(NF90_MAX_VAR_DIMS)
+    integer :: idimids_out(NF90_MAX_VAR_DIMS)
+    integer :: nattr
+    character(len=512) :: attr_name
+    character(len=512) :: dim_name
+
+    integer :: istatus
+    integer :: ivarid_in, ivarid_out
+    integer :: jattr, jdim
+
+    if (this%iverbose >= 4) then
+      write(nulout,'(a,a)') '  Copying definition of ', trim(var_name)
+    end if
+
+    ! Get variable ID from name
+    istatus = nf90_inq_varid(infile%ncid, var_name, ivarid_in) 
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,i0,a)') '*** Error inquiring about NetCDF variable "', &
+           & var_name, '": ', trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+    ! Get variable properties
+#ifdef NC_NETCDF4
+    istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type, ndims=ndims, &
+         &  dimids=idimids_in, chunksizes=chunksizes, deflate_level=deflate_level, &
+         &  shuffle=shuffle, natts=nattr)
+#else
+    istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type, ndims=ndims, &
+         &  dimids=idimids_in, natts=nattr)
+#endif
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a)') '*** Error reading NetCDF variable properties: ', &
+           trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+    ! Map dimension IDs
+    do jdim = 1,ndims
+      istatus = nf90_inquire_dimension(infile%ncid, idimids_in(jdim), name=dim_name)
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a)') '*** Error reading NetCDF dimension name: ', &
+             trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+
+      istatus = nf90_inq_dimid(this%ncid, trim(dim_name), idimids_out(jdim))
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a)') '*** Error reading NetCDF dimension ID: ', &
+             trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+    end do
+
+    ! Create variable
+#ifdef NC_NETCDF4
+    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids_out(1:ndims), &
+         & ivarid_out, deflate_level=deflate_level, shuffle=shuffle, chunksizes=chunksizes(1:ndims))
+#else
+    istatus = nf90_def_var(this%ncid, var_name, data_type, idimids_out(1:ndims), ivarid_out)
+#endif
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error defining variable "', var_name, &
+           &                    '": ', trim(nf90_strerror(istatus))
+      call my_abort('Error writing NetCDF file')
+    end if
+
+    ! Copy attributes
+    do jattr = 1,nattr
+      istatus = nf90_inq_attname(infile%ncid, ivarid_in, jattr, attr_name)
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a)') '*** Error reading attribute: ', &
+             &  trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+      istatus = nf90_copy_att(infile%ncid, ivarid_in, trim(attr_name), &
+           &                    this%ncid, ivarid_out)
+
+    end do
+
+  end subroutine copy_variable_definition
+
+
+  !---------------------------------------------------------------------
+  ! Copy variable from "infile" to "this"
+  subroutine copy_variable(this, infile, var_name)
+    class(netcdf_file)             :: this
+    class(netcdf_file), intent(in) :: infile
+    character(len=*),   intent(in) :: var_name
+
+    integer :: ivarid_in, ivarid_out
+    integer :: ndims
+    integer :: ndimlens(NF90_MAX_VAR_DIMS)
+    integer(kind=jpib) :: ntotal
+    integer :: data_type
+    integer :: istatus
+
+    ! We use the Fortran-77 functions because they don't check that
+    ! the rank of the arguments is correct
+    integer, external :: nf_get_var_double, nf_put_var_double
+    integer, external :: nf_get_var_int, nf_put_var_int
+
+    real(kind=jprd), allocatable :: data_real(:)
+    integer,         allocatable :: data_int(:)
+
+    ! If we are in define mode, exit define mode
+    call this%end_define_mode()
+
+    if (this%iverbose >= 4) then
+      write(nulout,'(a,a)') '  Copying ', trim(var_name)
+    end if
+
+    call infile%get_variable_id(var_name, ivarid_in)
+    call infile%get_array_dimensions(ivarid_in, ndims, ndimlens, ntotal)
+    istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type)
+    if (istatus /= NF90_NOERR) then
+      write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
+           &  trim(nf90_strerror(istatus))
+      call my_abort('Error reading NetCDF file')
+    end if
+
+    call infile%get_variable_id(var_name, ivarid_out)
+    if (data_type == NF90_DOUBLE .or. data_type == NF90_FLOAT) then
+      allocate(data_real(ntotal))
+      !istatus = nf90_get_var(infile%ncid, ivarid_in, data_real(1))
+      istatus = nf_get_var_double(infile%ncid, ivarid_in, data_real)
+      if (istatus /= NF90_NOERR) then
+        deallocate(data_real)
+        write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
+             &  trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+
+      !istatus = nf90_put_var(this%ncid, ivarid_out, data_real)
+      istatus = nf_put_var_double(this%ncid, ivarid_out, data_real)
+      deallocate(data_real)
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a,a,a)') '*** Error writing variable "', var_name, '": ', &
+             &  trim(nf90_strerror(istatus))
+        call my_abort('Error writing NetCDF file')
+      end if
+
+    else
+      allocate(data_int(ntotal))
+      !istatus = nf90_get_var(infile%ncid, ivarid_in, data_int)
+      istatus = nf_get_var_int(infile%ncid, ivarid_in, data_int)
+      if (istatus /= NF90_NOERR) then
+        deallocate(data_int)
+ 
+        write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
+             &  trim(nf90_strerror(istatus))
+        call my_abort('Error reading NetCDF file')
+      end if
+
+      !istatus = nf90_put_var(this%ncid, ivarid_out, data_int)
+      istatus = nf_put_var_int(this%ncid, ivarid_out, data_int)
+      deallocate(data_int)
+      if (istatus /= NF90_NOERR) then
+        write(nulerr,'(a,a,a,a)') '*** Error writing variable "', var_name, '": ', &
+             &  trim(nf90_strerror(istatus))
+        call my_abort('Error writing NetCDF file')
+      end if
+    end if
+
+  end subroutine copy_variable
+
 end module easy_netcdf