diff --git a/src/budget_interp_mod.f90 b/src/budget_interp_mod.f90
index 7af0e203..107e019f 100644
--- a/src/budget_interp_mod.f90
+++ b/src/budget_interp_mod.f90
@@ -126,9 +126,9 @@ SUBROUTINE interpolate_budget_scalar(IPOPT,grid_in,grid_out, &
 
     select type(grid_out)
     type is(ip_station_points_grid)
-       grid_desc_out2 = grid_out%descriptor
+       allocate(grid_desc_out2, source = grid_out%descriptor)
        grid_desc_out2%grid_num = 255 + grid_out%descriptor%grid_num
-       grid_out2 = init_grid(grid_desc_out2)
+       call init_grid(grid_out2, grid_desc_out2)
 
        CALL GDSWZD(grid_out2,-1,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
        IF(NO.EQ.0) then
@@ -460,10 +460,10 @@ SUBROUTINE interpolate_budget_vector(IPOPT,grid_in,grid_out, &
     ! The type of the subgrid is calculated by 255 + 
     select type(grid_out)
     type is(ip_station_points_grid)
-       desc_out_subgrid = grid_out%descriptor
+       allocate(desc_out_subgrid, source = grid_out%descriptor)
        desc_out_subgrid%grid_num = 255 + grid_out%descriptor%grid_num
 
-       grid_out2 = init_grid(desc_out_subgrid)
+       call init_grid(grid_out2, desc_out_subgrid)
        CALL GDSWZD(grid_out2,-1,MO,FILL,XPTS,YPTS, &
             RLON,RLAT,NO,CROT,SROT)
        IF(NO.EQ.0) IRET=3
diff --git a/src/gdswzd_mod.f90 b/src/gdswzd_mod.f90
index e0518204..4b14da18 100644
--- a/src/gdswzd_mod.f90
+++ b/src/gdswzd_mod.f90
@@ -679,7 +679,7 @@ SUBROUTINE GDSWZD_1D_ARRAY(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
     class(ip_grid), allocatable :: grid
 
     desc = init_descriptor(igdtnum, igdtlen, igdtmpl)
-    grid = init_grid(desc)
+    call init_grid(grid, desc)
     
     call gdswzd_grid(grid,IOPT,NPTS,FILL, &
          XPTS,YPTS,RLON,RLAT,NRET, &
@@ -771,7 +771,7 @@ SUBROUTINE GDSWZD_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, &
     class(ip_grid), allocatable :: grid
 
     desc = init_descriptor(kgds)
-    grid = init_grid(desc)
+    call init_grid(grid, desc)
     
     call gdswzd_grid(grid,IOPT,NPTS,FILL, &
          XPTS,YPTS,RLON,RLAT,NRET, &
@@ -866,7 +866,7 @@ SUBROUTINE GDSWZD_2d_array_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, &
     class(ip_grid), allocatable :: grid
 
     desc = init_descriptor(kgds)
-    grid = init_grid(desc)
+    call init_grid(grid, desc)
     
     call gdswzd_grid(grid,IOPT,NPTS,FILL, &
          XPTS,YPTS,RLON,RLAT,NRET, &
diff --git a/src/ip_grid_factory_mod.f90 b/src/ip_grid_factory_mod.f90
index 71ee6c14..1d30ba5d 100644
--- a/src/ip_grid_factory_mod.f90
+++ b/src/ip_grid_factory_mod.f90
@@ -25,34 +25,34 @@ module ip_grid_factory_mod
 
   !> Initializes a polymorphic ip_grid object from an ip_grid_descriptor.
   !!
+  !! @param[out] grid Grid to initialize
   !! @param[in] grid_desc Grid descriptor created from a grib1/grib2 template.
-  !! @return Initialized ip_grid.
   !!
   !! @author Kyle Gerheiser
   !! @date July 2021
-  function init_grid_generic(grid_desc) result(grid)
+  subroutine init_grid_generic(grid, grid_desc)
     class(ip_grid_descriptor), intent(in) :: grid_desc
-    class(ip_grid), allocatable :: grid
+    class(ip_grid), allocatable, intent(out) :: grid
 
     select type(grid_desc)
     type is(grib1_descriptor)
-       grid = init_grid_grib1(grid_desc)
+       call init_grid_grib1(grid, grid_desc)
     type is(grib2_descriptor)
-       grid = init_grid_grib2(grid_desc)
+       call init_grid_grib2(grid, grid_desc)
     end select
-  end function init_grid_generic
+  end subroutine init_grid_generic
 
   !> Initializes a polymorphic ip_grid from a grib1_descriptor.
   !! The concrete grid type is chosen based on the grid number in the descriptor.
   !!
+  !! @param[out] grid Grid to initialize
   !! @param[in] g1_desc
-  !! @return Initialized grid.
   !!
   !! @author Kyle Gerheiser
   !! @date July 2021
-  function init_grid_grib1(g1_desc) result(grid)
+  subroutine init_grid_grib1(grid, g1_desc)
     type(grib1_descriptor), intent(in) :: g1_desc
-    class(ip_grid), allocatable :: grid
+    class(ip_grid), allocatable, intent(out) :: grid
     
     select case(g1_desc%grid_num)
     case(:-1)
@@ -75,20 +75,20 @@ function init_grid_grib1(g1_desc) result(grid)
 
     call grid%init(g1_desc)
     allocate(grid%descriptor, source = g1_desc)
-  end function init_grid_grib1
+  end subroutine init_grid_grib1
 
   
   !> Initializes a polymorphic ip_grid from a grib2_descriptor.
   !! The concrete grid type is chosen based on the grid number in the descriptor.
   !!
-  !! @param[in] g2_desc
-  !! @return Initialized grid.
+  !! @param[out] grid Grid to initialize
+  !! @param[in] g2_desc Grib2 descriptor
   !!
   !! @author Kyle Gerheiser
   !! @date July 2021
-  function init_grid_grib2(g2_desc) result(grid)
+  subroutine init_grid_grib2(grid, g2_desc)
     type(grib2_descriptor), intent(in) :: g2_desc
-    class(ip_grid), allocatable :: grid
+    class(ip_grid), allocatable, intent(out) :: grid
 
     integer :: i_offset_odd, i_offset_even
 
@@ -120,6 +120,6 @@ function init_grid_grib2(g2_desc) result(grid)
 
     call grid%init(g2_desc)
     allocate(grid%descriptor, source = g2_desc)
-  end function init_grid_grib2
+  end subroutine init_grid_grib2
   
 end module ip_grid_factory_mod
diff --git a/src/ipolates.f90 b/src/ipolates.f90
index b990658b..b7e5ae78 100644
--- a/src/ipolates.f90
+++ b/src/ipolates.f90
@@ -9,7 +9,7 @@
 module ipolates_mod
   use ip_interpolators_mod
   use ip_grid_descriptor_mod
-  use ip_grid_factory_mod
+  use ip_grid_factory_mod, only: init_grid
   use ip_interpolators_mod
   use ip_grid_mod
   implicit none
@@ -207,8 +207,8 @@ subroutine ipolates_grib1(ip,ipopt,kgdsi,kgdso,mi,mo,km,ibi,li,gi, &
     desc_in = init_descriptor(kgdsi)
     desc_out = init_descriptor(kgdso)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
 
     call ipolates_grid(ip, ipopt, grid_in, grid_out, mi, mo, km, ibi, li, gi, no, rlat, rlon, ibo, lo, go, iret)
 
@@ -489,8 +489,8 @@ SUBROUTINE IPOLATES_grib2(IP,IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
     desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
     desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
 
     CALL ipolates_grid(ip,IPOPT,grid_in,grid_out,MI,MO,KM,IBI,LI,GI,NO,RLAT,RLON,IBO,LO,GO,IRET)
 
diff --git a/src/ipolatev.f90 b/src/ipolatev.f90
index d2cbaf1a..fe24159b 100644
--- a/src/ipolatev.f90
+++ b/src/ipolatev.f90
@@ -405,8 +405,8 @@ subroutine ipolatev_grib2(ip,ipopt,igdtnumi,igdtmpli,igdtleni, &
     desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
     desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
 
     CALL ipolatev_grid(ip,IPOPT,grid_in,grid_out, &
          MI,MO,KM,IBI,LI,UI,VI,&
@@ -577,8 +577,8 @@ subroutine ipolatev_grib1(ip,ipopt,kgdsi,kgdso,mi,mo,km,ibi,li,ui,vi, &
     desc_in = init_descriptor(kgdsi)
     desc_out = init_descriptor(kgdso)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
 
     CALL ipolatev_grid(ip,IPOPT,grid_in,grid_out, &
          MI,MO,KM,IBI,LI,UI,VI,&
diff --git a/src/spectral_interp_mod.f90 b/src/spectral_interp_mod.f90
index 1f34e964..94a29b2a 100644
--- a/src/spectral_interp_mod.f90
+++ b/src/spectral_interp_mod.f90
@@ -265,8 +265,8 @@ SUBROUTINE POLATES4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
     desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
     desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
     IRET=0
@@ -574,8 +574,9 @@ SUBROUTINE POLATES4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,GI, &
     desc_in = init_descriptor(kgdsi)
     desc_out = init_descriptor(kgdso)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
+    
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
     IRET=0
@@ -940,8 +941,8 @@ SUBROUTINE POLATEV4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
     desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
     desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
 
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
@@ -1277,8 +1278,8 @@ SUBROUTINE POLATEV4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,UI,VI, &
     desc_in = init_descriptor(kgdsi)
     desc_out = init_descriptor(kgdso)
 
-    grid_in = init_grid(desc_in)
-    grid_out = init_grid(desc_out)
+    call init_grid(grid_in, desc_in)
+    call init_grid(grid_out, desc_out)
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
     IRET=0