diff --git a/components/homme/src/share/physical_constants.F90 b/components/homme/src/share/physical_constants.F90
index 17c4ce9f6ff0..9080902fae7d 100644
--- a/components/homme/src/share/physical_constants.F90
+++ b/components/homme/src/share/physical_constants.F90
@@ -58,6 +58,7 @@ module physical_constants
   real (kind=real_kind), public, parameter :: rearth0      = 6.376D6    ! m
   real (kind=real_kind), public            :: rearth       = rearth0    ! m
   real (kind=real_kind), public, parameter :: gravit            = 9.80616D0  ! m s^-2
+  real (kind=real_kind), public, parameter :: g            = 9.80616D0  ! m s^-2
   real (kind=real_kind), public            :: gravitinv         = 1.0_real_kind/gravit
   real (kind=real_kind), public, parameter :: omega0       = 7.292D-5   ! s^-1
   real (kind=real_kind), public            :: omega        = omega0
diff --git a/components/homme/src/test_src/dcmip12_wrapper.F90 b/components/homme/src/test_src/dcmip12_wrapper.F90
index c8152b62b638..5536cbd3b3e6 100644
--- a/components/homme/src/test_src/dcmip12_wrapper.F90
+++ b/components/homme/src/test_src/dcmip12_wrapper.F90
@@ -91,7 +91,7 @@ subroutine dcmip2012_test1_1(elem,hybrid,hvcoord,nets,nete,time,n0,n1)
 
       dp = pressure_thickness(ps,k,hvcoord)
       call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
-      if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
+      if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
 
   enddo; enddo; enddo; enddo
 
@@ -167,7 +167,7 @@ subroutine dcmip2012_test1_1_conv(elem,hybrid,hvcoord,nets,nete,time,n0,n1)
 
       dp = pressure_thickness(ps,k,hvcoord)
       call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
-      if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
+      if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
 
   enddo; enddo; enddo; enddo
 
@@ -235,7 +235,7 @@ subroutine dcmip2012_test1_2(elem,hybrid,hvcoord,nets,nete,time,n0,n1)
       call test1_advection_hadley(time,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q(1),q(2))
       dp = pressure_thickness(ps,k,hvcoord)
       call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
-      if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
+      if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
 
   enddo; enddo; enddo; enddo
 
@@ -306,7 +306,7 @@ subroutine dcmip2012_test1_3(elem,hybrid,hvcoord,nets,nete,time,n0,n1,deriv)
       call test1_advection_orography(lon,lat,p,z,zcoords,cfv,use_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q(1),q(2),q(3),q(4))
       dp = pressure_thickness(ps,k,hvcoord)
       call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
-      if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
+      if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
 
   enddo; enddo; enddo; enddo
 
@@ -373,7 +373,7 @@ subroutine dcmip2012_test2_0(elem,hybrid,hvcoord,nets,nete)
         !let's get an analytical \phi
         he = (T0 - T)/gamma
         call set_state(u,v,w,T,ps,phis,p,dp,he,g, i,j,k,elem(ie),1,nt)
-        call set_tracers(q,qsize,dp,he,i,j,k,lat,lon,elem(ie))
+        call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
      enddo; enddo; enddo; 
      do k=1,nlevp; do j=1,np; do i=1,np
         call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord)
@@ -432,7 +432,7 @@ subroutine dcmip2012_test2_x(elem,hybrid,hvcoord,nets,nete,shear)
         !    call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt)
         ! This test obtains analytical height and returns it, so, we use it for \phi ...
         call set_state(u,v,w,T,ps,phis,p,dp,z,g, i,j,k,elem(ie),1,nt)
-        call set_tracers(q,qsize,dp,z,i,j,k,lat,lon,elem(ie))
+        call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
         ! ... or we can use discrete hydro state to init \phi. 
         
      enddo; enddo; enddo; 
@@ -602,7 +602,7 @@ subroutine dcmip2012_test3(elem,hybrid,hvcoord,nets,nete)
         call test3_gravity_wave(lon,lat,p,z,zcoords,use_eta,hyam,hybm,u,v,w,T,T_mean,phis,ps,rho,rho_mean,q(1))
         dp = pressure_thickness(ps,k,hvcoord)
         call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt)
-        call set_tracers(q,qsize, dp,z,i,j,k,lat,lon,elem(ie))
+        call set_tracers(q,qsize, dp,i,j,k,lat,lon,elem(ie))
      enddo; enddo; enddo; 
      do k=1,nlevp; do j=1,np; do i=1,np
         call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord)
@@ -662,7 +662,7 @@ subroutine dcmip2012_test4_init(elem,hybrid,hvcoord,nets,nete)
 
         !init only <=qsize tracers
         qs = min(qsize,3)
-        call set_tracers(qarray(1:qs),qs, dp,z,i,j,k,lat,lon,elem(ie))
+        call set_tracers(qarray(1:qs),qs, dp,i,j,k,lat,lon,elem(ie))
       enddo; enddo; enddo; 
 
     do k=1,nlevp
@@ -785,19 +785,19 @@ real(rl) function pressure_thickness(ps,k,hv)
 
 
 !_____________________________________________________________________
-subroutine set_tracers(q,nq, dp,zm,i,j,k,lat,lon,elem)
+subroutine set_tracers(q,nq, dp,i,j,k,lat,lon,elem)
   use physical_constants, only: rearth
   use deep_atm_mod, only: r_hat_from_phi, phi_from_z
   ! set tracer values at node(i,j,k)
 
-  real(rl),         intent(in)    :: q(nq), dp,zm, lat, lon
+  real(rl),         intent(in)    :: q(nq), dp, lat, lon
   integer,          intent(in)    :: i,j,k,nq
   type(element_t),  intent(inout) :: elem
   real(rl) :: r_hat, pdensity
   real(rl), parameter :: wl = 1.0 ! checkerboard wavelength in dg
   integer :: qi
 
-  r_hat = r_hat_from_phi(phi_from_z(zm))
+  !r_hat = r_hat_from_phi(phi_from_z(zm))
   pdensity =  dp
 
   if (nq>qsize) call abortmp('qsize set too small for dcmip test case')
diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90
index ddefd9a8f7ff..64a82e05245c 100644
--- a/components/homme/src/test_src/dcmip16_wrapper.F90
+++ b/components/homme/src/test_src/dcmip16_wrapper.F90
@@ -83,7 +83,8 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete)
   integer,            intent(in)            :: nets,nete                ! start, end element index
 
   integer,  parameter :: use_zcoords  = 0                               ! use vertical pressure coordinates
-  logical,  parameter :: is_deep      = atm_is_deep                             ! use shallow atmosphere approximation
+  !logical,  parameter :: is_deep      = atm_is_deep                             ! use shallow atmosphere approximation
+  integer,  parameter :: is_deep      = 1                             ! use shallow atmosphere approximation
   integer,  parameter :: pertt        = 0                               ! use exponential perturbation type
   real(rl), parameter :: dcmip_X      = 1.0_rl                          ! full scale planet
   integer :: moist                                                      ! use moist version
@@ -115,172 +116,48 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete)
 
   if (qsize<6) call abortmp('ERROR: test requires qsize>=6')
   
+  ! set initial conditions
   do ie = nets,nete
+    do k=1,nlevp; do j=1,np; do i=1,np
 
-  do k=1,nlev
-    dp(:,:, k) =-( (hvcoord%etai(k)-hvcoord%etai(k+1)) *  p0)
-  end do
+        ! no surface topography
+        p_i(i,j,k)  = p0*hvcoord%etai(k)
 
-  do k=1,nlevp; do j=1,np; do i=1,np
         lon = elem(ie)%spherep(i,j)%lon
         lat = elem(ie)%spherep(i,j)%lat
 
-
-         p_i(i,j,k)  = p0*hvcoord%etai(k)
-         call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),&
+        w_i(i,j,k)   = 0.0d0
+        ! call this only to compute z_i, will ignore all other output
+        call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),&
             z_i(i,j,k),use_zcoords,u(i,j,1),v(i,j,1),T(i,j,1),thetav(i,j,1),phis(i,j),ps(i,j),rho(i,j,1),q(i,j,1,1))
-         if (k > 1) then
-           dphi_tmp(i,j,k-1) = (phi_from_z(z_i(i,j,k-1)) - phi_from_z(z_i(i,j,k)))/r_hat_from_phi(phi_from_z(z_i(i,j,k)))**2
-         end if 
-  end do; end do; end do
-  q(:,:,:,1:5) = 0.0d0
-  q(:,:,:,6) = 1
-  w(:,:,:)   = 0.0d0
-  w_i(:,:,:)   = 0.0d0
-    
-
-  dphi_tmp = 10
-  
-  phi_prev(:,:,nlevp) = phis
-  z_i(:, :, nlevp) = z_from_phi(phis)
-  pnh_prev(:,:,nlevp) = ps
-  do k=nlev,1,-1; do j=1,np; do i=1,np
-      
-       lon = elem(ie)%spherep(i,j)%lon
-       lat = elem(ie)%spherep(i,j)%lat
-
-
-! ===================================
-      if (k.eq. nlev ) then
-         dmass_i(i,j,k) = -dp(i,j,k)/2
-       else
-         dmass_i(i,j,k) = -(dp(i,j,k)+dp(i,j,k+1))/2
-       end if
-       dmass(i,j,k) = -dp(i,j,k)
- 
-    phi_guess(i,j,k) = phi_prev(i,j,k+1) + dphi_tmp(i,j,k)  
-
-   if (mu_instead_of_eos_rootfinding) then
-    
-    do newton_idx=1, max_newton 
-       z_guess(i,j,k) = z_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k))/2)
-       z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1) +phi_guess(i,j,k)+newton_eps)/2)
-       call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
-                                 pnh_next(i,j,k),z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
-       call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
-                                 pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u_pert(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
-
-
-       r_hat(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k))/2)
-       r_hat_pert(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2)
-
-       pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1)))
-       pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1)))
-       !end if
-       residual(i,j,k) = -u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + &
-                          g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 )
-        if (abs(residual(i,j,k)) < 1e-10) then
-         exit
-        end if
-       deriv(i,j,k) =   ( (-u_pert(i,j,k)**2/(z_guess_pert(i,j,k)+rearth) - 2 * omega * cos(lat) * u_pert(i,j,k) + &
-                          g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next_pert(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 ))-&
-                          (-u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + &
-                          g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 )))/newton_eps
-                          
-
-  
-        if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then ! .or. pnh_next(i,j,k) > pnh_prev(i,j,k+1)) then
-           write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k
-           write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1)
-           write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1)
-           call abort
-        end if
-
-        phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k)
-     
-
-    end do
-
-    else
-    do newton_idx=1, max_newton
-       z_guess(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)
-       z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1)+phi_guess(i,j,k)+newton_eps)/2)
-       call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
-                                 p_tmp,z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
-       call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
-                                 pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
-
-       r_hat(i,j,k) = r_hat_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)
-       r_hat_pert(i,j,k) = r_hat_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2)
+    enddo; enddo; enddo
+    do k=1,nlev; do j=1,np; do i=1,np
 
+        ! no surface topography
+        p(i,j,k)  = p0*hvcoord%etam(k)
+        dp(i,j,k) = (hvcoord%etai(k+1)-hvcoord%etai(k))*p0
 
+        lon = elem(ie)%spherep(i,j)%lon
+        lat = elem(ie)%spherep(i,j)%lat
 
-       pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1)))
-       pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1)))
-       residual(i,j,k) = pnh_next(i,j,k) - p_tmp
-       if (abs(residual(i,j,k)) < 1e-10) then
-         exit
-        end if
-       deriv(i,j,k) = (pnh_next_pert(i,j,k) &
-                       - pnh_next(i,j,k))/newton_eps
-
-        if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then 
-           write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k
-           write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1)
-           write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1)
-           call abort
-        end if
-
-        phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k)
-
-
-    end do
+        q(i,j,k,1:5) = 0.0d0
+        q(i,j,k,6) = 1
+        w(i,j,k)   = 0.0d0
 
+        call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p(i,j,k),&
+            z(i,j,k),use_zcoords,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
 
+        ! initialize tracer chemistry
+        call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) )
+        call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),i,j,k,lat,lon,elem(ie))
 
-    end if
-    if (abs(residual(i,j,k)) > 1e-9) then
-      write(*,*)'WARNING: Deep atm geopotential rootfinding failed at level ', k, ' residual: ', abs(residual(i,j,k))
-      call abort
-    end if
-    pnh_prev(i,j,k) = pnh_next(i,j,k)
-    p(i,j,k) = pnh_next(i,j,k)
-    z_i(i, j, k) = z_from_phi(phi_guess(i, j,k)) 
-    z(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)
- 
-    phi_prev(i,j,k) = phi_guess(i,j,k)
-! ===================================
+        min_thetav =  min( min_thetav,   thetav(i,j,k) )
+        max_thetav =  max( max_thetav,   thetav(i,j,k) )
 
-  end do; end do; end do
-  !do k=1,nlev
-  ! write(*,*) "k: ", k,  " max pnh: ", maxval(p(:,:,k)), "min pnh: ", minval(p(:,:,k))
-  !end do
+    enddo; enddo; enddo
 
-        
-!  do k=1,nlev; do j=1,np; do i=1,np
-!        ! initialize tracer chemistry
-!        call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) )
-!        call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
-!
-!        min_thetav =  min( min_thetav,   thetav(i,j,k) )
-!        max_thetav =  max( max_thetav,   thetav(i,j,k) )
-!
-!    enddo; enddo; enddo
- 
     call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1)
     call tests_finalize(elem(ie),hvcoord)
-     
-    call pnh_and_exner_from_eos(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,nt),&
-       elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%phinh_i(:,:,:,nt),pnh,exner,dpnh_dp_i) 
-!    call pnh_and_exner_from_eos(hvcoord,T*dp*((p/p0)**(-kappa)),&
-!       dp,phi_from_z(z_i, nlevp),pnh,exner,dpnh_dp_i) 
-
-  do k=1,nlevp
-    if (maxval(abs(dpnh_dp_i(:,:,k)-1)) > 1e-9) then
-         write(*,*) "dpnh_dp_i error at lev ", k, ": ", maxval(dpnh_dp_i(:,:,k)-1), "minval: ", minval(dpnh_dp_i(:,:,k)-1)
-         !call abort
-    end if 
-  end do
 
   enddo
   sample_period = 1800.0 ! sec
@@ -392,7 +269,7 @@ subroutine dcmip2016_test2(elem,hybrid,hvcoord,nets,nete)
         q(2)=0
         q(3)=0
 
-        call set_tracers(q(:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
+        call set_tracers(q(:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie))
      enddo; enddo; enddo;
 
     call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1)
@@ -498,7 +375,7 @@ subroutine dcmip2016_test3(elem,hybrid,hvcoord,nets,nete)
         q   (i,j,k,2)= 0 ! no initial clouds
         q   (i,j,k,3)= 0 ! no initial rain
 
-        call set_tracers(q(i,j,k,:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
+        call set_tracers(q(i,j,k,:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie))
 
       enddo; enddo
     enddo
diff --git a/components/homme/src/test_src/dcmip2016-baroclinic.F90 b/components/homme/src/test_src/dcmip2016-baroclinic.F90
index db99434bc7f9..db02714c0042 100644
--- a/components/homme/src/test_src/dcmip2016-baroclinic.F90
+++ b/components/homme/src/test_src/dcmip2016-baroclinic.F90
@@ -40,7 +40,7 @@ MODULE baroclinic_wave
 !
 !=======================================================================
 
-use physical_constants, only: g0=>gravit,kappa0=>kappa,Rgas,Cp0=>Cp,Rwater_vapor,rearth0,omega0, dd_pi
+use physical_constants, only: g0=>g,kappa0=>kappa,Rgas,Cp0=>Cp,Rwater_vapor,rearth0,omega0, dd_pi
 
   IMPLICIT NONE
 
@@ -88,9 +88,9 @@ MODULE baroclinic_wave
        lapse      = 0.005d0             ! lapse rate parameter
 
   REAL(8), PARAMETER ::               &
-       pertu0     = 0.0d0      ,      & ! SF Perturbation wind velocity (m/s)
+       pertu0     = 0.5d0      ,      & ! SF Perturbation wind velocity (m/s)
        pertr      = 1.d0/6.d0  ,      & ! SF Perturbation radius (Earth radii)
-       pertup     = 0.0d0      ,      & ! Exp. perturbation wind velocity (m/s)
+       pertup     = 1.0d0      ,      & ! Exp. perturbation wind velocity (m/s)
        pertexpr   = 0.1d0      ,      & ! Exp. perturbation radius (Earth radii)
        pertlon    = pi/9.d0    ,      & ! Perturbation longitude
        pertlat    = 2.d0*pi/9.d0,     & ! Perturbation latitude
@@ -122,9 +122,9 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the
 !     input/output params parameters at given location
 !-----------------------------------------------------------------------
     INTEGER, INTENT(IN)  :: &
+                deep,       & ! Deep (1) or Shallow (0) test case
                 moist,      & ! Moist (1) or Dry (0) test case
                 pertt         ! Perturbation type
-    LOGICAL, INTENT(IN) :: deep
 
     REAL(8), INTENT(IN)  :: &
                 lon,        & ! Longitude (radians)
@@ -182,7 +182,7 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the
     inttau2 = constC * z * exp(- scaledZ**2)
 
     ! radius ratio
-    if (deep) then
+    if (deep .eq. 0) then
       rratio = 1.d0
     else
       rratio = (z + aref) / aref;
@@ -198,7 +198,7 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the
     !-----------------------------------------------------
     inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0)
     bigU = g / aref * K * inttau2 * inttermU * t
-    if (deep) then
+    if (deep .eq. 0) then
       rcoslat = aref * cos(lat)
     else
       rcoslat = (z + aref) * cos(lat)
@@ -270,7 +270,7 @@ END SUBROUTINE baroclinic_wave_test
 !-----------------------------------------------------------------------
   SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
 
-    LOGICAL, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
+    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
 
     REAL(8), INTENT(IN)  :: &
                 X,          & ! Earth scaling ratio
@@ -315,7 +315,7 @@ SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
     !--------------------------------------------
     !    radius ratio
     !--------------------------------------------
-    if (deep) then
+    if (deep .eq. 0) then
       rratio = 1.d0
     else
       rratio = (z + aref) / aref;
@@ -344,7 +344,7 @@ END SUBROUTINE evaluate_pressure_temperature
 !-----------------------------------------------------------------------
   SUBROUTINE evaluate_z_temperature(deep, X, lon, lat, p, z, t)
     
-    LOGICAL, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
+    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
 
     REAL(8), INTENT(IN)  :: &
                 X,          & ! Earth scaling ratio
diff --git a/components/homme/src/test_src/dry_planar.F90 b/components/homme/src/test_src/dry_planar.F90
index b3150ca3395d..3fa5df80105f 100644
--- a/components/homme/src/test_src/dry_planar.F90
+++ b/components/homme/src/test_src/dry_planar.F90
@@ -104,7 +104,7 @@ subroutine general_gravity_wave_init(elem,hybrid,hvcoord,nets,nete,d,f)
         call gravity_wave(x,y,p,z,zcoords,use_eta,hyam,hybm,d,u,v,w,T,T_mean,phis,ps,rho,rho_mean,q(1))
         dp = pressure_thickness(ps,k,hvcoord)
         call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt)
-        call set_tracers(q,qsize, dp,zm(k),i,j,k,y,x,elem(ie))
+        call set_tracers(q,qsize, dp,i,j,k,y,x,elem(ie))
      enddo; enddo; enddo;
      do k=1,nlevp; do j=1,np; do i=1,np
         call get_xycoordinates(x,y,hyai,hybi, i,j,k,elem(ie),hvcoord)
diff --git a/components/homme/src/theta-l/share/element_ops.F90 b/components/homme/src/theta-l/share/element_ops.F90
index d1bce092cf8e..8c12578adcaa 100644
--- a/components/homme/src/theta-l/share/element_ops.F90
+++ b/components/homme/src/theta-l/share/element_ops.F90
@@ -157,7 +157,7 @@ subroutine get_field_i(elem,name,field,hvcoord,nt)
       endif
     case('geo_i');
       if(theta_hydrostatic_mode) then
-         call phi_from_eos(hvcoord,elem%state%phis,elem%state%ps_v(:,:,nt),elem%state%vtheta_dp(:,:,:,nt),elem%state%dp3d(:,:,:,nt),field)
+         call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,nt),elem%state%dp3d(:,:,:,nt),field)
       else
           field = elem%state%phinh_i(:,:,1:nlevp,nt)
       endif
@@ -366,7 +366,7 @@ subroutine get_phi(elem,phi,phi_i,hvcoord,nt)
 
     if(theta_hydrostatic_mode) then
        dp=elem%state%dp3d(:,:,:,nt)
-       call phi_from_eos(hvcoord,elem%state%phis,elem%state%ps_v(:,:,nt),elem%state%vtheta_dp(:,:,:,nt),dp,phi_i)
+       call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,nt),dp,phi_i)
     else
        phi_i = elem%state%phinh_i(:,:,:,nt)
     endif
@@ -726,9 +726,11 @@ subroutine tests_finalize(elem,hvcoord,ie)
   logical, parameter :: do_rootfinding = .false., pnh_instead_of_mu=.false.
 
   tl=1
-  
-  !call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,tl),&
-  !     elem%state%dp3d(:,:,:,tl),elem%state%phinh_i(:,:,:,tl))
+ 
+!       call phi_from_eos(hvcoord,elem%state%phis,elem%state%ps_v(:,:,nt),elem%state%vtheta_dp(:,:,:,nt),dp,phi_i)
+ 
+  call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,tl),&
+       elem%state%dp3d(:,:,:,tl),elem%state%phinh_i(:,:,:,tl))
 
   ! Disable the following check in CUDA bfb builds,
   ! since the calls to pow are inexact
@@ -742,6 +744,7 @@ subroutine tests_finalize(elem,hvcoord,ie)
  
   ! =====================
    
+#if 0
   do k=1,nlev
      pi(:,:,k) = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem%state%ps_v(:,:,tl)
      if (maxval(abs(1-dpnh_dp_i(:,:,k))) > 1e-9) then
@@ -750,6 +753,7 @@ subroutine tests_finalize(elem,hvcoord,ie)
         write(iulog,*) 'pnh',pi(1,1,k),pnh(1,1,k)
      endif
   enddo
+#endif
 #endif
 
   do tl = 2,timelevels
@@ -796,7 +800,7 @@ subroutine initialize_reference_states(hvcoord, phis, &
 
    ! compute phi_ref
    temp = theta_ref*dp_ref
-   call phi_from_eos(hvcoord, phis, ps_ref, temp, dp_ref, phi_ref)
+   call phi_from_eos(hvcoord, phis, temp, dp_ref, phi_ref)
 
    ! keep profiles, based on the value of hv_ref_profiles
    if (hv_ref_profiles == 0) then
diff --git a/components/homme/src/theta-l/share/eos.F90 b/components/homme/src/theta-l/share/eos.F90
index 45b6f6dbe677..4a013819af92 100644
--- a/components/homme/src/theta-l/share/eos.F90
+++ b/components/homme/src/theta-l/share/eos.F90
@@ -78,7 +78,11 @@ subroutine pnh_and_exner_from_eos(hvcoord,vtheta_dp,dp3d,phi_i,pnh,exner,&
   end subroutine pnh_and_exner_from_eos
 
 
+#undef OGEOS
+!#define OGEOS
 
+#ifndef OGEOS
+!Owen's version
 subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
      dpnh_dp_i,caller,pnh_i_out)
 use deep_atm_mod, only: r_hat_from_phi, g_from_phi
@@ -98,8 +102,8 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
 ! instead of hvcoord%ps0, which is set by CAM to ~1021mb
 !  
   type (hvcoord_t),     intent(in)  :: hvcoord             ! hybrid vertical coordinate struct
-  real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev)   
-  real (kind=real_kind), intent(in) :: dp3d(np,np,nlev)   
+  real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev)
+  real (kind=real_kind), intent(in) :: dp3d(np,np,nlev)
   real (kind=real_kind), intent(in) :: dphi(np,np,nlev)
   real (kind=real_kind), intent(in) :: phi_i(np,np,nlevp) ! Needed for deep atmosphere
   real (kind=real_kind), intent(out) :: pnh(np,np,nlev)        ! nh nonhyrdo pressure
@@ -112,10 +116,10 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
   !   local
   real (kind=real_kind) :: p_over_exner(np,np,nlev)
   real (kind=real_kind) :: pi(np,np,nlev)
-  real (kind=real_kind) :: exner_i(np,np,nlevp) 
-  real (kind=real_kind) :: pnh_i(np,np,nlevp)  
+  real (kind=real_kind) :: exner_i(np,np,nlevp)
+  real (kind=real_kind) :: pnh_i(np,np,nlevp)
   real (kind=real_kind) :: dp3d_i(np,np,nlevp)
-  real (kind=real_kind) :: pi_i(np,np,nlevp) 
+  real (kind=real_kind) :: pi_i(np,np,nlevp)
   real (kind=real_kind) :: r_hat(np,np)
   integer :: i,j,k,k2
   logical :: ierr
@@ -169,8 +173,8 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
 
      pnh = pi ! copy hydrostatic pressure into output variable
      dpnh_dp_i = 1
-     if (present(pnh_i_out)) then  
-       pnh_i_out=pi_i 
+     if (present(pnh_i_out)) then
+       pnh_i_out=pi_i
      endif
   else
 
@@ -199,7 +203,13 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
    ! an approximation (hydrostatic) so that dpnh/dpi = 1
    ! DO NOT CHANGE this approximation.  it is required by 
    ! compute_andor_apply_rhs()
+
+!OWEN's
    pnh_i(:,:,1) = pnh(:,:,1) - dp3d(:, :,1)/(2.0 * r_hat**2) !* (1 + u_top/g_from_phi(phi_i(:, :, 1)))
+
+!OG pressure top
+!   pnh_i(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0/r_hat/r_hat
+
    pnh_i(:,:,nlevp) = pnh(:,:,nlev) + dp3d(:,:,nlev)/2
 
 
@@ -216,7 +226,7 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
    dpnh_dp_i(:,:,nlevp)  = 2*(pnh_i(:,:,nlevp)-pnh(:,:,nlev))/dp3d_i(:,:,nlevp) * r_hat**2
    do k=2,nlev
       r_hat = r_hat_from_phi(phi_i(:, :, k)) !DA_CHANGE
-      dpnh_dp_i(:,:,k) = (pnh(:,:,k) -pnh(:,:,k-1) )/dp3d_i(:,:,k) * r_hat**2 
+      dpnh_dp_i(:,:,k) = (pnh(:,:,k) -pnh(:,:,k-1) )/dp3d_i(:,:,k) * r_hat**2
    end do
    if (present(pnh_i_out)) then
       ! boundary values already computed. interpolate interior
@@ -226,14 +236,202 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
          pnh_i(:,:,k)=(dp3d(:,:,k-1)*pnh(:,:,k)+dp3d(:,:,k)*pnh(:,:,k-1))/&
               (dp3d(:,:,k-1)+dp3d(:,:,k))
       enddo
-      pnh_i_out=pnh_i    
+      pnh_i_out=pnh_i
    endif
-   
+
   endif ! hydrostatic/nonhydrostatic version
+  end subroutine
+#endif
+
+
+#ifdef OGEOS
+!OG version but with Owen's interface, phi_i instead of phis
+subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,&
+     dpnh_dp_i,caller,pnh_i_out)
+use deep_atm_mod, only: r_hat_from_phi, g_from_phi
+  use physical_constants, only : p0, kappa, gravit, Rgas, rearth
+implicit none
+  type (hvcoord_t),     intent(in)  :: hvcoord             ! hybrid vertical coordinate struct
+  real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev)   
+  real (kind=real_kind), intent(in) :: dp3d(np,np,nlev)   
+  real (kind=real_kind), intent(in) :: dphi(np,np,nlev)
+  real (kind=real_kind), intent(in) :: phi_i(np,np,nlevp) ! Needed for deep atmosphere
+  real (kind=real_kind), intent(out) :: pnh(np,np,nlev)        ! nh nonhyrdo pressure
+  real (kind=real_kind), intent(out) :: dpnh_dp_i(np,np,nlevp) ! d(pnh) / d(pi)
+  real (kind=real_kind), intent(out) :: exner(np,np,nlev)      ! exner nh pressure
+  character(len=*),      intent(in), optional  :: caller       ! name for error
+  real (kind=real_kind), intent(out), optional :: pnh_i_out(np,np,nlevp)  ! pnh on interfaces
+
+  !   local
+  real (kind=real_kind) :: p_over_exner(np,np,nlev)
+  real (kind=real_kind) :: pi(np,np,nlev)
+  real (kind=real_kind) :: exner_i(np,np,nlevp)
+  real (kind=real_kind) :: pnh_i(np,np,nlevp)
+  real (kind=real_kind) :: dp3d_i(np,np,nlevp)
+  real (kind=real_kind) :: pi_i(np,np,nlevp), g
+  integer :: i,j,k,k2
+  logical :: ierr
+
+  real (kind=real_kind) ::  rheighti(np,np,nlevp), rheightm(np,np,nlev), rhatm(np,np,nlev), r0
+  real (kind=real_kind) ::  rhati(np,np,nlevp), invrhatm(np,np,nlev), invrhati(np,np,nlevp), &
+                            newrhatsquared(np,np,nlev)
+
+  !construct phi_i here
+  !phi_i(:,:,nlevp) = phis(:,:)
+  !do k=nlev,1,-1
+  !  phi_i(:,:,k) = phi_i(:,:,k+1) - dphi(:,:,k)
+  !enddo
+
+  g = gravit
+  r0=rearth
+
+  rheighti = phi_i/g + r0
+  rheightm(:,:,1:nlev) = (rheighti(:,:,1:nlev) + rheighti(:,:,2:nlevp))/2.0
+  rhati = rheighti/r0 ! r/r0
+  rhatm = rheightm/r0
+  invrhatm = 1.0/rhatm
+  invrhati = 1.0/rhati
+
+  newrhatsquared = (rhati(:,:,1:nlev)*rhati(:,:,1:nlev)   + &
+                    rhati(:,:,2:nlevp)*rhati(:,:,2:nlevp) + &
+                    rhati(:,:,1:nlev)*rhati(:,:,2:nlevp))/3.0
+
+  ! check for bad state that will crash exponential function below
+  if (theta_hydrostatic_mode) then
+    ierr= any(dp3d(:,:,:) < 0 )
+  else
+    ierr= any(vtheta_dp(:,:,:) < 0 )  .or. &
+          any(dp3d(:,:,:) < 0 ) .or. &
+          any(dphi(:,:,:) > 0 )
+  endif
+
+  if (ierr) then
+     print *,'bad state in EOS, called from: ',caller
+     do j=1,np
+     do i=1,np
+     do k=1,nlev
+        if ( (vtheta_dp(i,j,k) < 0) .or. (dp3d(i,j,k)<0)  .or. &
+             (dphi(i,j,k)>0)  ) then
+           print *,'bad i,j,k=',i,j,k
+           print *,'vertical column: phi_i,dphi,dp3d,vtheta_dp'
+           do k2=1,nlev
+              write(*,'(i3,5f14.4)') k2,phi_i(i,j,k),dphi(i,j,k2),dp3d(i,j,k2),vtheta_dp(i,j,k2)
+           enddo
+print *, 'phi_i', phi_i(1,1,:)
+           call abortmp('EOS bad state: d(phi), dp3d or vtheta_dp < 0')
+        endif
+     enddo
+     enddo
+     enddo
+  endif
+
+  if (theta_hydrostatic_mode) then
+     ! hydrostatic pressure
+     pi_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0
+     do k=1,nlev
+        pi_i(:,:,k+1)=pi_i(:,:,k) + dp3d(:,:,k)
+     enddo
+#ifdef HOMMEXX_BFB_TESTING
+     do k=1,nlev
+        pi(:,:,k) = (pi_i(:,:,k+1)+pi_i(:,:,k))/2
+     enddo
+     exner  = bfb_pow(pi/p0,kappa)
+#else
+     do k=1,nlev
+        pi(:,:,k)=pi_i(:,:,k) + dp3d(:,:,k)/2
+     enddo
+     exner  = (pi/p0)**kappa
+#endif
+
+     pnh = pi ! copy hydrostatic pressure into output variable
+     dpnh_dp_i = 1
+     if (present(pnh_i_out)) then
+       pnh_i_out=pi_i
+     endif
+  else
+
+!==============================================================
+!  non-hydrostatic EOS
+!==============================================================
+  do k=1,nlev
+     p_over_exner(:,:,k) = Rgas*vtheta_dp(:,:,k)/(-dphi(:,:,k))
+
+     !da
+#ifdef DA
+     !p_over_exner(:,:,k) = p_over_exner(:,:,k)*invrhatm(:,:,k)*invrhatm(:,:,k)
+     p_over_exner(:,:,k) = p_over_exner(:,:,k)/newrhatsquared(:,:,k)
+
+!print *, 'newrhatsq', newrhatsquared(1,1,k)
+#endif
+
+#ifndef HOMMEXX_BFB_TESTING
+     pnh(:,:,k) = p0 * (p_over_exner(:,:,k)/p0)**(1/(1-kappa))
+#else
+     pnh(:,:,k) = p0 * bfb_pow(p_over_exner(:,:,k)/p0,1/(1-kappa))
+#endif
+     exner(:,:,k) =  pnh(:,:,k)/ p_over_exner(:,:,k)
+  enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! boundary terms
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
+
+!!!!!
+   pnh_i(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0*invrhati(:,:,1)*invrhati(:,:,1)  ! hydrostatic ptop/rhat^2
+
+   ! surface boundary condition pnh_i determined by w equation to enforce
+   ! w b.c.  This is computed in the RHS calculation.  Here, we use
+   ! an approximation (hydrostatic) so that dpnh/dpi = 1
+   ! DO NOT CHANGE this approximation.  it is required by 
+   ! compute_andor_apply_rhs()
+
+   pnh_i(:,:,nlevp) = pnh(:,:,nlev) + dp3d(:,:,nlev)/2
+!#ifdef DA
+!   pnh_i(:,:,nlevp) = pnh(:,:,nlev) + dp3d(:,:,nlev)/2*invrhatm(:,:,nlev)*invrhatm(:,:,nlev)
+!#endif
+
+   ! compute d(pnh)/d(pi) at interfaces
+   ! use one-sided differences at boundaries
+   dp3d_i(:,:,1) = dp3d(:,:,1)
+   dp3d_i(:,:,nlevp) = dp3d(:,:,nlev)
+   do k=2,nlev
+      dp3d_i(:,:,k)=(dp3d(:,:,k)+dp3d(:,:,k-1))/2
+   end do
+
+!original
+   dpnh_dp_i(:,:,1)  = 2*(pnh(:,:,1)-pnh_i(:,:,1))/dp3d_i(:,:,1)
+
+   dpnh_dp_i(:,:,nlevp)  = 2*(pnh_i(:,:,nlevp)-pnh(:,:,nlev))/dp3d_i(:,:,nlevp)
+   do k=2,nlev
+      dpnh_dp_i(:,:,k) = (pnh(:,:,k)-pnh(:,:,k-1))/dp3d_i(:,:,k)
+   end do
+
+#ifdef DA
+   !da
+   !keep the bottom val unchanged and set to 1
+   dpnh_dp_i(:,:,1:nlev) = dpnh_dp_i(:,:,1:nlev)*rhati(:,:,1:nlev)*rhati(:,:,1:nlev)
+#endif   
+
+   if (present(pnh_i_out)) then
+      ! boundary values already computed. interpolate interior
+      ! use linear interpolation in hydrostatic pressure coordinate
+      ! if pnh=pi, then pnh_i will recover pi_i
+      do k=2,nlev
+         pnh_i(:,:,k)=(dp3d(:,:,k-1)*pnh(:,:,k)+dp3d(:,:,k)*pnh(:,:,k-1))/&
+              (dp3d(:,:,k-1)+dp3d(:,:,k))
+      enddo
+      pnh_i_out=pnh_i
+   endif
+
+  endif ! hydrostatic/nonhydrostatic version
+
   end subroutine 
+#endif
+
 
 
 
+#if 0
+!!!! Owen's version with ps_v
   !_____________________________________________________________________
   subroutine phi_from_eos(hvcoord,phis,ps,vtheta_dp,dp,phi_i)
   use deep_atm_mod, only: r_hat_from_phi
@@ -318,5 +516,90 @@ subroutine phi_from_eos(hvcoord,phis,ps,vtheta_dp,dp,phi_i)
   enddo
   end subroutine
 
+#endif
+
+
+
+  subroutine phi_from_eos(hvcoord,phis,vtheta_dp,dp,phi_i)
+!
+! Use Equation of State to compute HYDROSTATIC geopotential
+!
+! input:  dp, phis, vtheta_dp  
+! output:  phi
+!
+! used to initialize phi for dry and wet test cases
+! used to compute background phi for reference state
+!
+! NOTE1: dp is pressure layer thickness.  If pnh is used to compute thickness, this
+! routine should be the discrete inverse of pnh_and_exner_from_eos().
+! This routine is usually called with hydrostatic layer thickness (dp3d), 
+! in which case it returns a hydrostatic PHI
+!
+! NOTE2: Exner pressure is defined in terms of p0=1000mb.  Be sure to use global constant p0,
+! instead of hvcoord%ps0, which is set by CAM to ~1021mb
+!  
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  implicit none
+
+  type (hvcoord_t),      intent(in)  :: hvcoord                      ! hybrid vertical coordinate struct
+  real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev)
+  real (kind=real_kind), intent(in) :: dp(np,np,nlev)
+  real (kind=real_kind), intent(in) :: phis(np,np)
+  real (kind=real_kind), intent(out) :: phi_i(np,np,nlevp)
+
+  !   local
+  real (kind=real_kind) :: p(np,np,nlev) ! pressure at cell centers 
+  real (kind=real_kind) :: p_i(np,np,nlevp)  ! pressure on interfaces
+
+  integer :: k
+
+#ifndef NDEBUG
+  logical :: ierr
+  integer :: i,j,k2
+
+  ierr= any(vtheta_dp(:,:,:) < 0 )  .or. &
+          any(dp(:,:,:) < 0 )
+
+  if (ierr) then
+     print *,'bad state in phi_from_eos:'
+     do j=1,np
+     do i=1,np
+     do k=1,nlev
+        if ( (vtheta_dp(i,j,k) < 0) .or. (dp(i,j,k)<0) ) then
+           print *,'bad i,j,k=',i,j,k
+           print *,'vertical column: dp,vtheta_dp'
+           do k2=1,nlev
+              write(*,'(i3,4f14.4)') k2,dp(i,j,k2),vtheta_dp(i,j,k2)
+           enddo
+           call abortmp('EOS bad state: dp or vtheta_dp < 0')
+        endif
+     enddo
+     enddo
+     enddo
+  endif
+#endif
+  ! compute pressure on interfaces                                                                                   
+  p_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0
+  do k=1,nlev
+     p_i(:,:,k+1)=p_i(:,:,k) + dp(:,:,k)
+  enddo
+  do k=1,nlev
+     p(:,:,k) = (p_i(:,:,k+1)+p_i(:,:,k))/2
+  enddo
+
+  phi_i(:,:,nlevp) = phis(:,:)
+
+  do k=nlev,1,-1
+#ifdef HOMMEXX_BFB_TESTING
+     phi_i(:,:,k) = phi_i(:,:,k+1)+ (Rgas*vtheta_dp(:,:,k)*bfb_pow(p(:,:,k)/p0,(kappa-1)))/p0
+#else
+     phi_i(:,:,k) = phi_i(:,:,k+1)+(Rgas*vtheta_dp(:,:,k)*(p(:,:,k)/p0)**(kappa-1))/p0
+#endif
+  enddo
+
+  end subroutine phi_from_eos
+
+
+
 end module
 
diff --git a/components/homme/src/theta-l/share/imex_mod.F90 b/components/homme/src/theta-l/share/imex_mod.F90
index b6898be7254e..a2b0930391b7 100644
--- a/components/homme/src/theta-l/share/imex_mod.F90
+++ b/components/homme/src/theta-l/share/imex_mod.F90
@@ -145,7 +145,7 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h
        itercount=0
        itererr=0
        do ie=nets,nete
-          call phi_from_eos(hvcoord,elem(ie)%state%phis, elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),&
+          call phi_from_eos(hvcoord,elem(ie)%state%phis, elem(ie)%state%vtheta_dp(:,:,:,np1),&
                elem(ie)%state%dp3d(:,:,:,np1),elem(ie)%state%phinh_i(:,:,:,np1))
           elem(ie)%state%w_i(:,:,:,np1)=0
        enddo
@@ -227,7 +227,7 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h
 #endif
 #if 1
        ! use hydrostatic for initial guess
-       call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),&
+       call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),&
             elem(ie)%state%dp3d(:,:,:,np1),phi_np1)
        w_np1(:,:,1:nlev) = (z_from_phi(phi_np1(:,:,1:nlev), nlev) -  z_from_phi(phi_n0(:,:,1:nlev), nlev ))/(dt2) ! DA_CHANGE
 #endif
@@ -361,12 +361,13 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h
 #endif
        itererr=min(max_reserr,max_deltaerr) ! passed back to ARKODE
 
-
+#if 0
        if (itercount >= maxiter) then
           if (verbosity > 0) then
           write(iulog,*) 'WARNING:IMEX solver failed b/c max iteration count was met',deltaerr,reserr
           end if
        end if
+#endif
     end do ! end do for the ie=nets,nete loop
 #ifdef NEWTONCOND
     if (hybrid%masterthread) print *,'max J condition number (mpi task0): ',1/min_rcond
@@ -401,7 +402,7 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
     real (kind=real_kind), intent(out) :: JacD(np,np,nlev)
     real (kind=real_kind), intent(out) :: JacL(np,np,nlev-1),JacU(np,np,nlev-1)
     real (kind=real_kind), intent(in)  :: dp3d(np,np,nlev)
-    real (kind=real_kind), intent(inout) :: pnh(np,np,nlev)
+    real (kind=real_kind), intent(inout)  :: pnh(np,np,nlev)
     real (kind=real_kind), intent(in)  :: dphi(np,np,nlev)
     real (kind=real_kind), intent(in)  :: phi_i(np,np,nlevp)
     real (kind=real_kind), intent(in)  :: dt2
@@ -420,11 +421,41 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
     real (kind=real_kind) :: ds(np,np,nlev),delta_mu(np,np,nlevp),r_hat(np,np,nlevp)
     real (kind=real_kind) :: a(np,np),b(np,np),ck(np,np),ckm1(np,np)
     !
+!original SA jacobian
+!#define ORIGINALI
+
+!new approx DA jacobian from Owen
+#undef ORIGINALI
+
     integer :: k,l,k2
     if (exact.eq.1) then ! use exact JacobianA
 
        ! this code will need to change when the equation of state is changed.
        ! add special cases for k==1 and k==nlev+1
+#ifdef ORIGINALI
+       a  = (dt2*gravit)**2/(1-kappa)
+       k  = 1 ! Jacobian row 1
+       b  = a/dp3d(:,:,k)
+       ck = pnh(:,:,k)/dphi(:,:,k)
+       JacU(:,:,k) = 2*b*ck
+       JacD(:,:,k) = 1 - JacU(:,:,k)
+       ckm1 = ck
+       do k = 2,nlev-1 ! Jacobian row k
+          b  = 2*a/(dp3d(:,:,k-1) + dp3d(:,:,k))
+          ck = pnh(:,:,k)/dphi(:,:,k)
+          JacL(:,:,k-1) = b*ckm1
+          JacU(:,:,k  ) = b*ck
+          JacD(:,:,k  ) = 1 - JacL(:,:,k-1) - JacU(:,:,k)
+          ckm1 = ck
+       end do
+       k  = nlev ! Jacobian row nlev
+       b  = 2*a/(dp3d(:,:,k) + dp3d(:,:,k-1))
+       ck = pnh(:,:,k)/dphi(:,:,k)
+       JacL(:,:,k-1) = b*ckm1
+       JacD(:,:,k  ) = 1 - JacL(:,:,k-1) - b*ck
+#else
+       ! this code will need to change when the equation of state is changed.
+       ! add special cases for k==1 and k==nlev+1
        phi_i_tmp = phi_i
        k  = 1 ! Jacobian row 1
        r_hat(:,:,k) = r_hat_from_phi(phi_i_tmp(:,:,k)) ! DA_CHANGE
@@ -453,6 +484,8 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
        ck = pnh(:,:,k)/dphi(:,:,k)
        JacL(:,:,k-1) = r_hat(:,:,k)**2 * b*ckm1
        JacD(:,:,k  ) = 1 - JacL(:,:,k-1) - r_hat(:,:,k)**2 * b*ck
+#endif
+
 
     else ! use finite difference approximation to Jacobian with differencing size espie
       ! compute Jacobian of F(dphi) = phi +const + (dt*g)^2 *(1-dp/dpi) column wise
@@ -461,6 +494,36 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
       ! we only form the tridagonal entries and this code can easily be modified to
       ! accomodate sparse non-tridigonal and dense Jacobians, however, testing only
       ! the tridiagonal of a Jacobian is probably sufficient for testing purpose
+#ifdef ORIGINALI
+      do k=1,nlev
+        e=0
+        e(:,:,k)=1
+        dphi_temp(:,:,:) = dphi(:,:,:)
+        ! dphi_tmp(k)= ( phi(k+1)-(phi(k)+eps ) = dphi(k)-eps
+        ! dphi_tmp(k-1)= ( phi(k)+eps-(phi(k-1) ) = dphi(k-1)+eps
+        dphi_temp(:,:,k) = dphi(:,:,k) - epsie*e(:,:,k)
+        if (k>1) dphi_temp(:,:,k-1) = dphi(:,:,k-1) + epsie*e(:,:,k)
+
+        if (theta_hydrostatic_mode) then
+           !dpnh_dp_i_epsie(:,:,:)=1.d0
+           delta_mu=0
+        else
+!hvcoord,vtheta_dp,dp3d,dphi_temp,phi_i_tmp, pnh,exner,dpnh_dp_i_epsie,'get_dirk_jacobian'
+           call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi_temp,phi_i,pnh,exner,dpnh_dp_i_epsie,'get_dirk_jacobian')
+           delta_mu(:,:,:)=(gravit*dt2)**2*(dpnh_dp_i(:,:,:)-dpnh_dp_i_epsie(:,:,:))/epsie
+        end if
+
+        JacD(:,:,k) = 1 +  delta_mu(:,:,k)
+        if (k.eq.1) then
+           JacL(:,:,k) =   delta_mu(:,:,k+1)
+        elseif (k.eq.nlev) then
+           JacU(:,:,k-1) = delta_mu(:,:,k-1)
+        else
+           JacL(:,:,k)   = delta_mu(:,:,k+1)
+           JacU(:,:,k-1) = delta_mu(:,:,k-1)
+        end if
+      end do
+#else
       phi_i_tmp = phi_i
       do k=1,nlev
         e=0
@@ -478,7 +541,7 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
            do k2=nlev,1,-1  ! scan
               phi_i_tmp(:,:,k2) = phi_i_tmp(:,:,k2+1)-dphi_temp(:,:,k2)
            enddo
- 
+
            call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi_temp,phi_i_tmp, pnh,exner,dpnh_dp_i_epsie,'get_dirk_jacobian')
            r_hat = r_hat_from_phi(phi_i_tmp, nlevp) !DA_CHANGE
            delta_mu(:,:,:)= (g_from_phi(phi_i_tmp,nlevp)*dt2)**2*(dpnh_dp_i(:,:,:)-dpnh_dp_i_epsie(:,:,:))/epsie ! DA_CHANGE
@@ -495,6 +558,8 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,&
         end if
 
       end do
+
+#endif
     end if
   end subroutine get_dirk_jacobian
 
diff --git a/components/homme/src/theta-l/share/vertremap_mod.F90 b/components/homme/src/theta-l/share/vertremap_mod.F90
index cd98b29d5794..0ad1dd4185bc 100644
--- a/components/homme/src/theta-l/share/vertremap_mod.F90
+++ b/components/homme/src/theta-l/share/vertremap_mod.F90
@@ -107,7 +107,7 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete)
 
      if (rsplit>0) then
         ! remove hydrostatic phi befor remap
-        call phi_from_eos(hvcoord,elem(ie)%state%phis, elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),dp_star,phi_ref)
+        call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),dp_star,phi_ref)
         elem(ie)%state%phinh_i(:,:,:,np1)=&
              elem(ie)%state%phinh_i(:,:,:,np1) -phi_ref(:,:,:)
  
@@ -141,7 +141,7 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete)
         enddo
 
         ! depends on theta, so do this after updating theta:
-        call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),dp,phi_ref)
+        call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),dp,phi_ref)
         elem(ie)%state%phinh_i(:,:,:,np1)=&
              elem(ie)%state%phinh_i(:,:,:,np1)+phi_ref(:,:,:)