From 1949c64092e9eb653399031d33f83ecd6395754d Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 15 Jun 2016 16:20:33 +1000 Subject: [PATCH] Newton's method returns error if finds a stationary point. Change root guess if this is the case. Closes #65. --- src/shared/cvmix_kpp.F90 | 18 +++++++++++++++--- src/shared/cvmix_math.F90 | 11 ++++++++++- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/shared/cvmix_kpp.F90 b/src/shared/cvmix_kpp.F90 index e5b675521..f8fa8d4b2 100644 --- a/src/shared/cvmix_kpp.F90 +++ b/src/shared/cvmix_kpp.F90 @@ -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 @@ -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 @@ -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, diff --git a/src/shared/cvmix_math.F90 b/src/shared/cvmix_math.F90 index cdf6fcdca..d89870440 100644 --- a/src/shared/cvmix_math.F90 +++ b/src/shared/cvmix_math.F90 @@ -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 @@ -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