From 5d153ea28941ea0e5980b3a853b0bb6805fdd6ec Mon Sep 17 00:00:00 2001
From: Nick Wogan <nicholaswogan@gmail.com>
Date: Wed, 19 Jun 2024 16:56:33 -0700
Subject: [PATCH] added perturbations to root solve

---
 src/adiabat/clima_adiabat_solve.f90 | 57 +++++++++++++++++++----------
 1 file changed, 38 insertions(+), 19 deletions(-)

diff --git a/src/adiabat/clima_adiabat_solve.f90 b/src/adiabat/clima_adiabat_solve.f90
index 45f52ab..5d200c7 100644
--- a/src/adiabat/clima_adiabat_solve.f90
+++ b/src/adiabat/clima_adiabat_solve.f90
@@ -96,7 +96,8 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec
     type(MinpackHybrj) :: mv
     logical, allocatable :: convecting_with_below_save(:,:)
     real(dp), allocatable :: difference(:)
-    real(dp), allocatable :: T_in(:)
+    real(dp), allocatable :: T_in(:), x_init(:)
+    real(dp) :: perturbation
     integer :: i, j, k
 
     if (.not.self%double_radiative_grid) then
@@ -132,30 +133,48 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec
     do i = 1,self%max_rc_iters
       j = i
 
-      mv = MinpackHybrj(fcn, size(self%inds_Tx))
-      mv%x(1) = self%T_surf
+      if (self%verbose) then
+        print"(1x,'Iteration =',i3)", i
+      endif
+
+      if (allocated(x_init)) deallocate(x_init)
+      allocate(x_init(size(self%inds_Tx)))
+      x_init(1) = self%T_surf
       do k = 2,size(self%inds_Tx)
-        mv%x(k) = self%T(self%inds_Tx(k)-1)
+        x_init(k) = self%T(self%inds_Tx(k)-1)
       enddo
-      mv%nprint = 1 ! enable printing
+
+      mv = MinpackHybrj(fcn, size(self%inds_Tx))
       mv%xtol = self%xtol_rc
+      mv%nprint = 1
 
-      if (self%verbose) then
-        print"(1x,'Iteration =',i3)", i
-      endif
+      k = 0
+      do
+        if (mod(k,2) == 0) then
+          perturbation = real(k,dp)*1.0_dp
+        else
+          perturbation = -real(k,dp)*1.0_dp
+        endif
 
-      call mv%hybrj()
-      if (mv%info == 0) then
-        err = 'hybrj root solve failed in surface_temperature.'
-        return
-      elseif (any(mv%info == [2, 3, 4, 5])) then
-        if (self%verbose) then
-          print'(3x,A)','Minpack Warning: '//mv%code_to_message(mv%info)
+        if (self%verbose .and. k > 0) then
+          print'(3x,"Perturbation = ",f7.1)',perturbation
         endif
-      elseif (mv%info < 0) then
-        err = 'hybrj root solve failed in surface_temperature: '//err
-        return
-      endif
+
+        mv%x = x_init + perturbation
+        call mv%hybrj()
+        if (mv%info == 1) then
+          exit
+        else
+          if (mv%info < 0) deallocate(err)
+        endif
+
+        if (k > 6) then
+          err = 'hybrj root solve failed in RCE.'
+          return
+        endif
+
+        k = k + 1
+      enddo
 
       ! Update all variables to the current root
       call AdiabatClimate_objective(self, P_i_surf, mv%x, .false., mv%fvec, err)