Skip to content

Commit

Permalink
added perturbations to root solve
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholaswogan committed Jun 19, 2024
1 parent 122916e commit 5d153ea
Showing 1 changed file with 38 additions and 19 deletions.
57 changes: 38 additions & 19 deletions src/adiabat/clima_adiabat_solve.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 5d153ea

Please sign in to comment.