diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90
index e217196c..fbf64a67 100644
--- a/libglide/glide_types.F90
+++ b/libglide/glide_types.F90
@@ -2293,6 +2293,10 @@ subroutine glide_allocarr(model)
        call coordsystem_allocate(model%general%ice_grid, model%inversion%powerlaw_c_prescribed)
        call coordsystem_allocate(model%general%ice_grid, model%inversion%usrf_inversion)
        call coordsystem_allocate(model%general%ice_grid, model%inversion%dthck_dt_inversion)
+    else
+       ! Allocate powerlaw_c_inversion with size 1, since it is passed into subroutine calcbeta
+       !  by the velocity solver.
+       allocate(model%inversion%powerlaw_c_inversion(1,1))
     endif
 
     ! climate arrays
diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90
index d7aeab16..0e5f0279 100644
--- a/libglissade/glissade_velo_higher.F90
+++ b/libglissade/glissade_velo_higher.F90
@@ -1979,7 +1979,11 @@ subroutine glissade_velo_higher_solve(model,                &
     do j = 1+nhalo, ny-nhalo+1
        do i = 1+nhalo, nx-nhalo+1
           ! gn = exponent in Glen's flow law (= 3 by default)
-          flwafact(:,i,j) = 0.5d0 * flwa(:,i,j)**(-1.d0/real(gn,dp))  
+          do k = 1, nz-1
+             if (flwa(k,i,j) > 0.0d0) then
+                flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/real(gn,dp))  
+             endif
+          enddo
        enddo
     enddo