Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Newton's method returns error if finds a stationary point. Change roo… #66

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Newton's method returns error if finds a stationary point. Change roo… #66

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 15 additions & 3 deletions src/shared/cvmix_kpp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ module cvmix_kpp
use cvmix_math, only : CVMIX_MATH_INTERP_LINEAR, &
CVMIX_MATH_INTERP_QUAD, &
CVMIX_MATH_INTERP_CUBE_SPLINE, &
CVMIX_MATH_ROOT_NOT_FOUND, &
cvmix_math_poly_interp, &
cvmix_math_cubic_root_find, &
cvmix_math_evaluate_cubic
Expand Down Expand Up @@ -1302,7 +1303,8 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
real(kind=cvmix_r8), dimension(:), pointer :: depth
real(kind=cvmix_r8), dimension(4) :: coeffs
real(kind=cvmix_r8) :: Ekman, MoninObukhov, OBL_Limit
integer :: nlev, k
real(kind=cvmix_r8), dimension(3) :: root_guess_mul
integer :: nlev, k, r, errval
logical :: lstable

type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in
Expand Down Expand Up @@ -1400,8 +1402,18 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
end if
coeffs(1) = coeffs(1)-CVmix_kpp_params_in%ri_crit

OBL_depth = -cvmix_math_cubic_root_find(coeffs, 0.5_cvmix_r8 * &
(depth(k)+depth(k+1)))
root_guess_mul = (/ 0.5_cvmix_r8, 0.25_cvmix_r8, 0.75_cvmix_r8 /)
do r=1,size(root_guess_mul)
OBL_depth = -cvmix_math_cubic_root_find(coeffs, root_guess_mul(r) * &
(depth(k)+depth(k+1)), errval)
if (errval /= CVMIX_MATH_ROOT_NOT_FOUND) then
exit
endif
enddo
if (errval == CVMIX_MATH_ROOT_NOT_FOUND) then
print*, "ERROR: cvmix_math_cubic_root_find could not find root."
stop 1
endif

! OBL_depth needs to be at or below the center of the top level
! Note: OBL_depth can only be computed to be above this point if k=1,
Expand Down
11 changes: 10 additions & 1 deletion src/shared/cvmix_math.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,14 @@ module cvmix_math
real(cvmix_r8), parameter :: CVMIX_MATH_NEWTON_TOL = 1.0e-12_cvmix_r8
integer, parameter :: CVMIX_MATH_MAX_NEWTON_ITERS = 100

integer, parameter :: CVMIX_MATH_ROOT_NOT_FOUND = 1

! !PUBLIC MEMBER FUNCTIONS:

public :: cvmix_math_poly_interp
public :: cvmix_math_cubic_root_find
public :: cvmix_math_evaluate_cubic
public :: CVMIX_MATH_ROOT_NOT_FOUND

!EOP

Expand Down Expand Up @@ -180,21 +183,27 @@ subroutine cvmix_math_poly_interp(coeffs, interp_type, x, y, x0, y0)

end subroutine cvmix_math_poly_interp

function cvmix_math_cubic_root_find(coeffs, x0)
function cvmix_math_cubic_root_find(coeffs, x0, errval)

real(cvmix_r8), dimension(4), intent(in) :: coeffs
real(cvmix_r8), intent(in) :: x0
integer, intent(out) :: errval

real(cvmix_r8) :: cvmix_math_cubic_root_find
real(cvmix_r8) :: fun_val, root, slope
integer :: it_cnt

errval = 0
root = x0
fun_val = coeffs(4)*(root**3)+coeffs(3)*(root**2)+coeffs(2)*root+coeffs(1)
do it_cnt = 1, CVMIX_MATH_MAX_NEWTON_ITERS
if (abs(fun_val).lt.CVMIX_MATH_NEWTON_TOL) &
exit
slope = 3.0_cvmix_r8*coeffs(4)*(root**2)+2.0_cvmix_r8*coeffs(3)*root+coeffs(2)
if (slope == 0) then
errval = CVMIX_MATH_ROOT_NOT_FOUND
exit
endif
root = root - fun_val/slope
fun_val = coeffs(4)*(root**3)+coeffs(3)*(root**2)+coeffs(2)*root+coeffs(1)
end do
Expand Down