diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ad86f09..95d1bf5 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -25,6 +25,7 @@ set(RADIATE_SOURCES
 set(ADIABAT_SOURCES
   adiabat/clima_adiabat_general.f90
   adiabat/clima_adiabat_rc.f90
+  adiabat/clima_adiabat_solve.f90
   adiabat/clima_adiabat.f90
 )
 
diff --git a/src/adiabat/clima_adiabat.f90 b/src/adiabat/clima_adiabat.f90
index 34d3d4d..62d83eb 100644
--- a/src/adiabat/clima_adiabat.f90
+++ b/src/adiabat/clima_adiabat.f90
@@ -10,6 +10,15 @@ module clima_adiabat
 
   public :: AdiabatClimate
 
+  type :: RCEquilibriumWrk
+    integer :: n_convecting_zones
+    integer, allocatable :: ind_conv_lower(:)
+    integer, allocatable :: ind_conv_upper(:)
+
+    real(dp), allocatable :: lapse_rate(:)
+    real(dp), allocatable :: lapse_rate_obs(:)
+  end type
+
   type :: AdiabatClimate
 
     ! settings and free parameters
@@ -68,6 +77,9 @@ module clima_adiabat
     real(dp), allocatable :: densities_r(:,:)
     real(dp), allocatable :: dz_r(:)
 
+    ! work
+    type(RCEquilibriumWrk), allocatable :: rcwrk
+
     ! tolerances
     !> Relative tolerance of integration
     real(dp) :: rtol = 1.0e-12_dp
@@ -75,6 +87,8 @@ module clima_adiabat
     real(dp) :: atol = 1.0e-12_dp
     !> Tolerance for nonlinear solve in make_column
     real(dp) :: tol_make_column = 1.0e-8_dp
+    !> Perturbation for the jacobian
+    real(dp) :: epsj = 1.0e-2_dp
     
     ! State of the atmosphere
     real(dp) :: P_surf !! Surface pressure (dynes/cm^2)
@@ -105,8 +119,9 @@ module clima_adiabat
     procedure :: surface_temperature => AdiabatClimate_surface_temperature
     procedure :: surface_temperature_column => AdiabatClimate_surface_temperature_column
     procedure :: surface_temperature_bg_gas => AdiabatClimate_surface_temperature_bg_gas
-    ! Test
-    procedure :: make_profile_rc => AdiabatClimate_make_profile_rc
+    ! Routines for full radiative convective equilibrium
+    procedure, private :: make_profile_rc => AdiabatClimate_make_profile_rc
+    procedure :: RCE => AdiabatClimate_RCE
     ! Utilities
     procedure :: set_ocean_solubility_fcn => AdiabatClimate_set_ocean_solubility_fcn
     procedure :: to_regular_grid => AdiabatClimate_to_regular_grid
@@ -121,6 +136,21 @@ module clima_adiabat
   interface AdiabatClimate
     module procedure :: create_AdiabatClimate
   end interface
+
+  interface
+    module subroutine AdiabatClimate_make_profile_rc(self, P_i_surf, T_surf, T, err)
+      class(AdiabatClimate), intent(inout) :: self
+      real(dp), intent(in) :: P_i_surf(:) !! dynes/cm^2
+      real(dp), intent(in) :: T_surf, T(:)
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+    module subroutine AdiabatClimate_RCE(self, P_i_surf, T_guess, err)
+      class(AdiabatClimate), intent(inout) :: self
+      real(dp), intent(in) :: P_i_surf(:)
+      real(dp), intent(in) :: T_guess
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+  end interface
   
 contains
   
@@ -204,6 +234,9 @@ function create_AdiabatClimate(species_f, settings_f, star_f, data_dir, double_r
     c%rad = Radtran(c%species_names, particle_names, s, star_f, s%number_of_zenith_angles, s%surface_albedo, c%nz_r, data_dir, err)
     if (allocated(err)) return
 
+    allocate(c%rcwrk)
+    allocate(c%rcwrk%lapse_rate(c%nz),c%rcwrk%lapse_rate_obs(c%nz))
+
     ! allocate ocean functions
     allocate(c%ocean_fcns(c%sp%ng))
 
@@ -767,64 +800,6 @@ subroutine fcn(n_, x_, fvec_, iflag_)
     end subroutine
   end function
 
-  subroutine AdiabatClimate_make_profile_rc(self, P_i_surf, T_surf, T, err)
-    use clima_const, only: k_boltz, N_avo
-    use clima_adiabat_rc, only: make_profile_rc
-    class(AdiabatClimate), intent(inout) :: self
-    real(dp), intent(in) :: P_i_surf(:) !! dynes/cm^2
-    real(dp), intent(in) :: T_surf, T(:)
-    character(:), allocatable, intent(out) :: err
-
-    real(dp), allocatable :: P_e(:), z_e(:), f_i_e(:,:), lapse_rate_e(:)
-    real(dp), allocatable :: density(:)
-    integer :: i, j
-
-    allocate(P_e(2*self%nz+1),z_e(2*self%nz+1),f_i_e(2*self%nz+1,self%sp%ng),lapse_rate_e(2*self%nz+1))
-    allocate(density(self%nz))
-
-    if (size(P_i_surf) /= self%sp%ng) then
-      err = "P_i_surf has the wrong dimension"
-      return
-    endif
-    if (size(T) /= self%nz) then
-      err = "T has the wrong dimension"
-    endif
-
-    self%T_surf = T_surf
-    self%T = T
-    call make_profile_rc(self%T_surf, self%T, P_i_surf, &
-                         self%sp, self%nz, self%planet_mass, &
-                         self%planet_radius, self%P_top, self%RH, &
-                         self%rtol, self%atol, &
-                         self%ocean_fcns, self%ocean_args_p, &
-                         P_e, z_e, f_i_e, lapse_rate_e, &
-                         self%N_surface, self%N_ocean, &
-                         err)
-    if (allocated(err)) return
-
-    self%T_surf = self%T_surf
-    self%P_surf = P_e(1)
-    do i = 1,self%nz
-      self%P(i) = P_e(2*i)
-      self%z(i) = z_e(2*i)
-      self%dz(i) = z_e(2*i+1) - z_e(2*i-1)
-      do j =1,self%sp%ng
-        self%f_i(i,j) = f_i_e(2*i,j)
-      enddo
-    enddo
-
-    density = self%P/(k_boltz*self%T)
-    do j =1,self%sp%ng
-      self%densities(:,j) = self%f_i(:,j)*density(:)
-    enddo
-
-    do i = 1,self%sp%ng
-      ! mol/cm^2 in atmosphere
-      self%N_atmos(i) = sum(density*self%f_i(:,i)*self%dz)/N_avo
-    enddo
-
-  end subroutine
-
   !> Sets a function for describing how gases dissolve in a liquid ocean.
   subroutine AdiabatClimate_set_ocean_solubility_fcn(self, species, fcn, err)
     class(AdiabatClimate), intent(inout) :: self
diff --git a/src/adiabat/clima_adiabat_solve.f90 b/src/adiabat/clima_adiabat_solve.f90
new file mode 100644
index 0000000..8363e65
--- /dev/null
+++ b/src/adiabat/clima_adiabat_solve.f90
@@ -0,0 +1,292 @@
+
+submodule(clima_adiabat) clima_adiabat_solve
+  implicit none
+
+contains
+
+  module subroutine AdiabatClimate_make_profile_rc(self, P_i_surf, T_surf, T, err)
+    use clima_const, only: k_boltz, N_avo
+    use clima_adiabat_rc, only: make_profile_rc
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(in) :: P_i_surf(:) !! dynes/cm^2
+    real(dp), intent(in) :: T_surf, T(:)
+    character(:), allocatable, intent(out) :: err
+
+    real(dp), allocatable :: P_e(:), z_e(:), f_i_e(:,:), lapse_rate_e(:)
+    real(dp), allocatable :: density(:)
+    integer :: i, j
+
+    allocate(P_e(2*self%nz+1),z_e(2*self%nz+1),f_i_e(2*self%nz+1,self%sp%ng),lapse_rate_e(2*self%nz+1))
+    allocate(density(self%nz))
+
+    if (size(P_i_surf) /= self%sp%ng) then
+      err = "P_i_surf has the wrong dimension"
+      return
+    endif
+    if (size(T) /= self%nz) then
+      err = "T has the wrong dimension"
+    endif
+
+    self%T_surf = T_surf
+    self%T = T
+    call make_profile_rc(self%T_surf, self%T, P_i_surf, &
+                         self%sp, self%nz, self%planet_mass, &
+                         self%planet_radius, self%P_top, self%RH, &
+                         self%rtol, self%atol, &
+                         self%ocean_fcns, self%ocean_args_p, &
+                         P_e, z_e, f_i_e, lapse_rate_e, &
+                         self%N_surface, self%N_ocean, &
+                         err)
+    if (allocated(err)) return
+
+    self%P_surf = P_e(1)
+    do i = 1,self%nz
+      self%P(i) = P_e(2*i)
+      self%z(i) = z_e(2*i)
+      self%dz(i) = z_e(2*i+1) - z_e(2*i-1)
+      do j =1,self%sp%ng
+        self%f_i(i,j) = f_i_e(2*i,j)
+      enddo
+    enddo
+
+    density = self%P/(k_boltz*self%T)
+    do j =1,self%sp%ng
+      self%densities(:,j) = self%f_i(:,j)*density(:)
+    enddo
+
+    do i = 1,self%sp%ng
+      ! mol/cm^2 in atmosphere
+      self%N_atmos(i) = sum(density*self%f_i(:,i)*self%dz)/N_avo
+    enddo
+
+    ! lapse rates
+    self%rcwrk%lapse_rate(1) = lapse_rate_e(1)
+    do i = 2,self%nz
+      self%rcwrk%lapse_rate(i) = lapse_rate_e(2*i-1)
+    enddo
+
+    self%rcwrk%lapse_rate_obs(1) = (log(self%T(1)) - log(self%T_surf))/(log(self%P(1)) - log(self%P_surf))
+    do i = 2,self%nz
+      self%rcwrk%lapse_rate_obs(i) = (log(self%T(i)) - log(self%T(i-1)))/(log(self%P(i)) - log(self%P(i-1)))
+    enddo
+
+  end subroutine
+
+  module subroutine AdiabatClimate_RCE(self, P_i_surf, T_guess, err)
+    use minpack_module, only: hybrd1
+    use clima_useful, only: MinpackHybrj
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(in) :: P_i_surf(:)
+    real(dp), intent(in) :: T_guess
+    character(:), allocatable, intent(out) :: err
+
+    real(dp) :: T_surf
+    type(MinpackHybrj) :: mv
+
+    ! Create the initial temperature profile
+    T_surf = self%surface_temperature(P_i_surf, T_guess, err)
+    if (allocated(err)) return
+
+    ! call self%make_profile(T_guess, P_i_surf, err)
+    ! if (allocated(err)) return
+
+    ! Establish the convective zone
+    if (self%P_trop > 0.0_dp) then
+      self%rcwrk%n_convecting_zones = 1
+      allocate(self%rcwrk%ind_conv_lower(1))
+      allocate(self%rcwrk%ind_conv_upper(1))
+      self%rcwrk%ind_conv_lower(1) = 1
+      self%rcwrk%ind_conv_upper(1) = minloc(abs(self%P_trop - [self%P_surf, self%P]),dim=1)
+    else
+      err = 'error!'
+      return
+    endif
+
+    open(unit=2,file='test.dat',form='unformatted',status='replace')
+    write(unit=2) self%nz
+    write(unit=2) self%sp%ng
+  
+    ! Non-linear solve
+    mv = MinpackHybrj(fcn,self%nz+1)
+    mv%x(1) = self%T_surf
+    mv%x(2:) = self%T
+    mv%nprint = 1 ! enable printing
+    mv%xtol = 1.0e-8_dp
+    call mv%hybrj()
+    if (mv%info == 0 .or. mv%info > 1) then
+      print*,mv%info
+      err = 'hybrj root solve failed in surface_temperature.'
+      return
+    elseif (mv%info < 0) then
+      err = 'hybrj root solve failed in surface_temperature: '//err
+      return
+    endif
+
+    close(2)
+
+  contains
+
+    subroutine fcn(n_, x_, fvec_, fjac_, ldfjac_, iflag_)
+      implicit none
+      integer, intent(in) :: n_
+      real(dp), dimension(n_), intent(in) :: x_
+      integer, intent(in) :: ldfjac_
+      real(dp), dimension(n_), intent(inout) :: fvec_
+      real(dp), dimension(ldfjac_, n_), intent(inout) :: fjac_
+      integer, intent(inout) :: iflag_
+
+      if (allocated(err)) return
+
+      if (iflag_ == 1) then
+        ! Compute right-hand-side
+        call AdiabatClimate_objective(self, P_i_surf, x_, fvec_, err)
+        if (allocated(err)) then
+          iflag_ = -1
+          return
+        endif
+      elseif (iflag_ == 2) then
+        ! Compute jacobian
+        call AdiabatClimate_jacobian(self, P_i_surf, x_, fjac_, err)
+        if (allocated(err)) then
+          iflag_ = -1
+          return
+        endif
+      endif
+
+      if (iflag_ == 0) then
+        print*,maxval(x_),minval(x_),sum(fvec_**2.0_dp),mv%nfev, mv%njev
+        write(2) [self%P_surf,self%P]
+        write(2) x_
+        write(2) fvec_
+        write(2) sum(fvec_**2.0_dp)
+        write(2) self%f_i
+      endif
+
+    end subroutine
+
+  end subroutine
+
+  subroutine AdiabatClimate_objective(self, P_i_surf, T_in, res, err)
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(in) :: P_i_surf(:)
+    real(dp), intent(in) :: T_in(:)
+    real(dp), intent(out) :: res(:)
+    character(:), allocatable, intent(out) :: err
+
+    ! Sets, self%T_surf, self%T, self%N_surface, self%N_ocean, self%P_surf, self%P,
+    ! self%z, self%dz, self%f_i, self%densities, self%N_atmos, self%rcwrk%lapse_rate
+    ! self%rcwrk%lapse_rate_obs
+    self%T_surf = T_in(1)
+    self%T = T_in(2:)
+    call self%make_profile_rc(P_i_surf, self%T_surf, self%T, err)
+    if (allocated(err)) return
+
+    ! resets self%T_surf, self%T, self%densities, self%rcwrk%lapse_rate_obs
+    ! also does radiative transfer and computes res
+    call AdiabatClimate_objective_(self, P_i_surf, T_in, res, err)
+    if (allocated(err)) return
+
+  end subroutine
+
+  subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, res, err)
+    use clima_const, only: k_boltz
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(in) :: P_i_surf(:)
+    real(dp), intent(in) :: T_in(:)
+    real(dp), intent(out) :: res(:)
+    character(:), allocatable, intent(out) :: err
+
+    real(dp), allocatable :: fluxes(:)
+    real(dp), allocatable :: density(:)
+    integer :: i, j
+
+    ! work storage
+    allocate(fluxes(self%nz+1))
+    allocate(density(self%nz))
+
+    ! Stuff set in make_profile_rc that does not change
+    ! - self%N_surface, self%N_ocean, self%P_surf, self%P, self%z, 
+    ! - self%dz, self%f_i, self%N_atmos, self%rcwrk%lapse_rate
+    ! Stuff that does change
+    ! - self%T_surf, self%T, self%densities, self%rcwrk%lapse_rate_obs
+    self%T_surf = T_in(1)
+    self%T = T_in(2:)
+
+    density = self%P/(k_boltz*self%T)
+    do j =1,self%sp%ng
+      self%densities(:,j) = self%f_i(:,j)*density(:)
+    enddo
+
+    self%rcwrk%lapse_rate_obs(1) = (log(self%T(1)) - log(self%T_surf))/(log(self%P(1)) - log(self%P_surf))
+    do i = 2,self%nz
+      self%rcwrk%lapse_rate_obs(i) = (log(self%T(i)) - log(self%T(i-1)))/(log(self%P(i)) - log(self%P(i-1)))
+    enddo
+
+    ! radiative transfer
+    call self%copy_atm_to_radiative_grid()
+    call self%rad%radiate(self%T_surf, self%T_r, self%P_r/1.0e6_dp, self%densities_r, self%dz_r, err=err)
+    if (allocated(err)) return
+
+    ! Energy going into each layer (ergs/(cm^2*s))
+    fluxes(1) = self%rad%f_total(1)
+    do i = 1,self%nz
+      fluxes(i+1) = (self%rad%f_total(2*(i-1)+3) - self%rad%f_total(2*(i-1)+1))
+    enddo
+
+    ! Radiative layers
+    res(j+1:) = fluxes(j+1:)
+
+    ! Convective layers
+    j = self%rcwrk%ind_conv_upper(1)
+    i = j - 1 ! i is the atmospheric layer
+    res(j) = self%rad%f_total(2*(i-1)+3)
+    do i = 1,j-1
+      res(i) = self%rcwrk%lapse_rate_obs(i) - self%rcwrk%lapse_rate(i)
+    enddo
+
+  end subroutine
+
+  subroutine AdiabatClimate_jacobian(self, P_i_surf, T_in, jac, err)
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(in) :: P_i_surf(:)
+    real(dp), intent(in) :: T_in(:)
+    real(dp), intent(out) :: jac(:,:)
+    character(:), allocatable, intent(out) :: err
+
+    real(dp), allocatable :: res(:)
+    real(dp), allocatable :: T_perturb(:), res_perturb(:)
+    real(dp) :: deltaT
+
+    integer :: i
+
+    ! allocate work
+    allocate(res(self%nz+1))
+    allocate(T_perturb(self%nz+1),res_perturb(self%nz+1))
+
+    ! First evaluate res at T.
+    call AdiabatClimate_objective(self, P_i_surf, T_in, res, err)
+    if (allocated(err)) return
+
+    ! Now compute Jacobian
+    T_perturb = T_in
+    do i = 1,size(T_perturb)
+
+      ! Perturb temperature
+      deltaT = self%epsj*abs(T_in(i))
+      T_perturb(i) = T_in(i) + deltaT
+
+      ! Compute rhs_perturb
+      call AdiabatClimate_objective_(self, P_i_surf, T_perturb, res_perturb, err)
+      if (allocated(err)) return
+
+      ! Compute jacobian
+      jac(:,i) = (res_perturb(:) - res(:))/deltaT
+      
+      ! unperturb T
+      T_perturb(i) = T_in(i)
+    enddo
+
+  end subroutine
+
+
+end submodule
\ No newline at end of file
diff --git a/src/clima_useful.f90 b/src/clima_useful.f90
index d1add53..78ce164 100644
--- a/src/clima_useful.f90
+++ b/src/clima_useful.f90
@@ -1,9 +1,10 @@
 module clima_useful
   use clima_const, only: dp
+  use minpack_module, only: fcn_hybrj
   implicit none
   private
 
-  public :: MinpackHybrd1Vars
+  public :: MinpackHybrd1Vars, MinpackHybrj
 
   type :: MinpackHybrd1Vars
     integer :: n
@@ -18,6 +19,34 @@ module clima_useful
     procedure :: create_MinpackHybrd1Vars
   end interface
 
+  type :: MinpackHybrj
+    procedure(fcn_hybrj), nopass, pointer :: fcn => null()
+    integer :: n
+    real(dp), allocatable :: x(:)
+    real(dp), allocatable :: fvec(:)
+    real(dp), allocatable :: fjac(:,:)
+    integer :: Ldfjac
+    real(dp) :: xtol
+    integer :: maxfev
+    real(dp), allocatable :: diag(:)
+    integer :: mode
+    real(dp) :: factor
+    integer :: nprint
+    integer :: info
+    integer :: nfev
+    integer :: njev
+    real(dp), allocatable :: r(:)
+    integer :: Lr
+    real(dp), allocatable :: qtf(:)
+    real(dp), allocatable :: wa1(:), wa2(:), wa3(:), wa4(:)
+    integer :: lwa
+  contains
+    procedure :: hybrj => MinpackHybrj_hybrj
+  end type
+  interface MinpackHybrj
+    procedure :: create_MinpackHybrj
+  end interface
+
 contains
 
   function create_MinpackHybrd1Vars(n, tol) result(v)
@@ -34,4 +63,44 @@ function create_MinpackHybrd1Vars(n, tol) result(v)
     allocate(v%wa(v%lwa))
   end function
 
+  function create_MinpackHybrj(fcn, n) result(v)
+    procedure(fcn_hybrj) :: fcn
+    integer, intent(in) :: n
+    type(MinpackHybrj) :: v
+
+    v%fcn => fcn
+    v%n = n
+    allocate(v%x(v%n))
+    v%x = 1.0_dp
+    allocate(v%fvec(v%n))
+    allocate(v%fjac(v%n,v%n))
+    v%ldfjac = v%n
+    v%xtol = 1.0e-8_dp ! not default in hybrdj1
+    v%maxfev = 100*(v%n + 1)
+    allocate(v%diag(v%n))
+    v%diag(:) = 1.0_dp
+    v%mode = 2
+    v%factor = 100.0_dp
+    v%nprint = 0
+    v%lr = (v%n*(v%n+1))/2
+    allocate(v%r(v%lr))
+    allocate(v%qtf(v%n))
+    allocate(v%wa1(v%n), v%wa2(v%n), v%wa3(v%n), v%wa4(v%n))
+
+  end function
+
+  subroutine MinpackHybrj_hybrj(self)
+    use minpack_module, only: hybrj
+    class(MinpackHybrj), intent(inout) :: self
+
+    if (.not.associated(self%fcn)) return
+
+    call hybrj(self%fcn, self%n, self%x, self%Fvec, self%Fjac, &
+               self%Ldfjac, self%Xtol, self%Maxfev, self%Diag, self%Mode, &
+               self%Factor, self%Nprint, self%Info, self%Nfev, self%Njev, &
+               self%r, self%Lr, self%Qtf, self%Wa1, self%Wa2, &
+               self%Wa3, self%Wa4)
+    
+  end subroutine
+
 end module
\ No newline at end of file