add netcdf read and write support

Author: mjreno

autotest/test_netcdf_gwf_vsc03_sfr.py

@@ -227,7 +227,7 @@ def check_output(idx, test, export, gridded_input):
     aux = getattr(rch, "aux").array
     if export == "ugrid":
         rarr = xds["rcha-1_recharge_l1_p1"].data.flatten()
-        auxarr = xds["rcha-1_auxvar_l1_p1a1"].data.flatten()
+        auxarr = xds["rcha-1_temperature_l1_p1"].data.flatten()
     elif export == "structured":
         rarr = xds["rcha-1_recharge_p1"].data[0].flatten()
         auxarr = xds["rcha-1_temperature_p1"].data[0].flatten()

src/Utilities/Export/DisNCMesh.f90

@@ -26,7 +26,6 @@ module MeshModelModule
   public :: ncvar_deflate
   public :: ncvar_gridmap
   public :: ncvar_mf6attr
-  public :: export_varname
 
   !> @brief type for storing model export dimension ids
   !<
@@ -564,35 +563,4 @@ subroutine ncvar_mf6attr(ncid, varid, layer, iper, iaux, nc_tag, nc_fname)
     end if
   end subroutine ncvar_mf6attr
 
-  !> @brief build netcdf variable name
-  !<
-  function export_varname(varname, layer, iper, iaux) result(vname)
-    use InputOutputModule, only: lowcase
-    character(len=*), intent(in) :: varname
-    integer(I4B), optional, intent(in) :: layer
-    integer(I4B), optional, intent(in) :: iper
-    integer(I4B), optional, intent(in) :: iaux
-    character(len=LINELENGTH) :: vname
-    vname = ''
-    if (varname /= '') then
-      vname = varname
-      call lowcase(vname)
-      if (present(layer)) then
-        if (layer > 0) then
-          write (vname, '(a,i0)') trim(vname)//'_l', layer
-        end if
-      end if
-      if (present(iper)) then
-        if (iper > 0) then
-          write (vname, '(a,i0)') trim(vname)//'_p', iper
-        end if
-      end if
-      if (present(iaux)) then
-        if (iaux > 0) then
-          write (vname, '(a,i0)') trim(vname)//'a', iaux
-        end if
-      end if
-    end if
-  end function export_varname
-
 end module MeshModelModule

src/Utilities/Export/DisNCStructured.f90

@@ -24,7 +24,7 @@ module NCModelExportModule
   public :: NCExportAnnotation
   public :: ExportPackageType
   public :: NETCDF_UNDEF, NETCDF_STRUCTURED, NETCDF_MESH2D
-  public :: export_longname
+  public :: export_longname, export_varname
 
   !> @brief netcdf export types enumerator
   !<
@@ -410,14 +410,62 @@ function input_attribute(this, pkgname, idt) result(attr)
     end if
   end function input_attribute
 
+  !> @brief build netcdf variable name
+  !<
+  function export_varname(pkgname, tagname, mempath, layer, iper, iaux) &
+    result(varname)
+    use MemoryManagerModule, only: mem_setptr
+    use CharacterStringModule, only: CharacterStringType
+    use InputOutputModule, only: lowcase
+    character(len=*), intent(in) :: pkgname
+    character(len=*), intent(in) :: tagname
+    character(len=*), intent(in) :: mempath
+    integer(I4B), optional, intent(in) :: layer
+    integer(I4B), optional, intent(in) :: iper
+    integer(I4B), optional, intent(in) :: iaux
+    character(len=LINELENGTH) :: varname
+    type(CharacterStringType), dimension(:), pointer, &
+      contiguous :: auxnames
+    character(len=LINELENGTH) :: pname, vname
+    vname = tagname
+    pname = pkgname
+
+    if (present(iaux)) then
+      if (iaux > 0) then
+        if (tagname == 'AUX') then
+          ! reset vname to auxiliary variable name
+          call mem_setptr(auxnames, 'AUXILIARY', mempath)
+          vname = auxnames(iaux)
+        end if
+      end if
+    end if
+
+    call lowcase(vname)
+    call lowcase(pname)
+    varname = trim(pname)//'_'//trim(vname)
+
+    if (present(layer)) then
+      if (layer > 0) then
+        !write (varname, '(a,i0)') trim(varname)//'_L', layer
+        write (varname, '(a,i0)') trim(varname)//'_l', layer
+      end if
+    end if
+    if (present(iper)) then
+      if (iper > 0) then
+        !write (varname, '(a,i0)') trim(varname)//'_SP', iper
+        write (varname, '(a,i0)') trim(varname)//'_p', iper
+      end if
+    end if
+  end function export_varname
+
   !> @brief build netcdf variable longname
   !<
   function export_longname(longname, pkgname, tagname, layer, iper) result(lname)
     use InputOutputModule, only: lowcase
     character(len=*), intent(in) :: longname
     character(len=*), intent(in) :: pkgname
     character(len=*), intent(in) :: tagname
-    integer(I4B), intent(in) :: layer
+    integer(I4B), optional, intent(in) :: layer
     integer(I4B), optional, intent(in) :: iper
     character(len=LINELENGTH) :: lname
     character(len=LINELENGTH) :: pname, vname
@@ -430,8 +478,10 @@ function export_longname(longname, pkgname, tagname, layer, iper) result(lname)
     else
       lname = longname
     end if
-    if (layer > 0) then
-      write (lname, '(a,i0)') trim(lname)//' layer=', layer
+    if (present(layer)) then
+      if (layer > 0) then
+        write (lname, '(a,i0)') trim(lname)//' layer=', layer
+      end if
     end if
     if (present(iper)) then
       if (iper > 0) then

src/Utilities/Export/DisvNCMesh.f90

@@ -462,11 +462,11 @@ subroutine load_double1d_spd(dbl1d, mf6_input, mshape, idt, nc_vars, &
     character(len=*), intent(in) :: input_fname
     integer(I4B), optional, intent(in) :: iaux
     real(DP), dimension(:, :, :), contiguous, pointer :: dbl3d
-    integer(I4B) :: nvals, varid
+    integer(I4B) :: varid
     integer(I4B) :: n, i, j, k
 
     ! initialize
-    nvals = 0
+    n = 0
 
     ! set varid
     if (present(iaux)) then
@@ -475,20 +475,22 @@ subroutine load_double1d_spd(dbl1d, mf6_input, mshape, idt, nc_vars, &
       varid = nc_vars%varid(idt%mf6varname, period=iper)
     end if
 
-    if (idt%shape == 'NODES') then
-      ! TODO future support
-      write (errmsg, '(a)') &
-        'IDM NetCDF load_double1d_spd NODES var shape not supported => '// &
-        trim(idt%tagname)
-      call store_error(errmsg)
-      call store_error_filename(input_fname)
-    else if (idt%shape == 'NCPL' .or. idt%shape == 'NAUX NCPL') then
+    if (size(mshape) == 3) then
+      allocate (dbl3d(mshape(3), mshape(2), mshape(1)))
+      call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl3d), &
+                     nc_vars%nc_fname)
 
-      if (size(mshape) == 3) then
-        allocate (dbl3d(mshape(3), mshape(2), mshape(1)))
-        call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl3d), &
-                       nc_vars%nc_fname)
-        n = 0
+      if (idt%shape == 'NODES' .or. idt%shape == 'NAUX NODES') then
+        do k = 1, size(dbl3d, dim=3)
+          do i = 1, size(dbl3d, dim=2)
+            do j = 1, size(dbl3d, dim=1)
+              n = n + 1
+              dbl1d(n) = dbl3d(j, i, k)
+            end do
+          end do
+        end do
+
+      else if (idt%shape == 'NCPL' .or. idt%shape == 'NAUX NCPL') then
         do k = 1, size(dbl3d, dim=3)
           do i = 1, size(dbl3d, dim=2)
             do j = 1, size(dbl3d, dim=1)
@@ -503,15 +505,10 @@ subroutine load_double1d_spd(dbl1d, mf6_input, mshape, idt, nc_vars, &
             end do
           end do
         end do
-
-      else if (size(mshape) == 2) then
-        ! TODO
-        write (errmsg, '(a)') &
-          'IDM NetCDF load_double1d_spd DISV model not supported => '// &
-          trim(idt%tagname)
-        call store_error(errmsg)
-        call store_error_filename(input_fname)
       end if
+
+      ! clean up
+      deallocate (dbl3d)
     end if
   end subroutine load_double1d_spd
 
@@ -561,29 +558,33 @@ subroutine load_double1d_layered_spd(dbl1d, mf6_input, mshape, idt, nc_vars, &
     integer(I4B), optional, intent(in) :: iaux
     integer(I4B), dimension(:), allocatable :: layer_shape
     integer(I4B) :: nlay, varid
-    integer(I4B) :: k, n, ncpl
-    integer(I4B) :: index_start, index_stop
+    integer(I4B) :: k, n, ncpl, idx
     real(DP), dimension(:), contiguous, pointer :: dbl1d_ptr
 
     call get_layered_shape(mshape, nlay, layer_shape)
     ncpl = product(layer_shape)
     allocate (dbl1d_ptr(ncpl))
 
     do k = 1, nlay
-      index_start = 1
-      index_stop = index_start + ncpl - 1
       if (present(iaux)) then
         varid = nc_vars%varid(idt%mf6varname, layer=k, period=iper, iaux=iaux)
       else
         varid = nc_vars%varid(idt%mf6varname, layer=k, period=iper)
       end if
       call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl1d_ptr), &
                      nc_vars%nc_fname)
-      do n = 1, ncpl
-        if (dbl1d_ptr(n) /= DNODATA) then
-          dbl1d(n) = dbl1d_ptr(n)
-        end if
-      end do
+      if (idt%shape == 'NODES' .or. idt%shape == 'NAUX NODES') then
+        do n = 1, ncpl
+          idx = (k - 1) * ncpl + n
+          dbl1d(idx) = dbl1d_ptr(n)
+        end do
+      else if (idt%shape == 'NCPL' .or. idt%shape == 'NAUX NCPL') then
+        do n = 1, ncpl
+          if (dbl1d_ptr(n) /= DNODATA) then
+            dbl1d(n) = dbl1d_ptr(n)
+          end if
+        end do
+      end if
     end do
 
     ! cleanup

src/Utilities/Export/MeshNCModel.f90

@@ -116,7 +116,7 @@ function ncvars_varid(this, tagname, layer, period, iaux) result(varid)
         write (errmsg, '(a)') &
           'NetCDF variable not found, tagname="'//trim(tagname)//'"'
         if (present(layer)) then
-          write (errmsg, '(a,i0)') trim(errmsg)//', ilayer=', layer
+          write (errmsg, '(a,i0)') trim(errmsg)//', layer=', layer
         end if
         if (present(period)) then
           write (errmsg, '(a,i0)') trim(errmsg)//', period=', period