polynomial_functions.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Polynomial functions
7
8implicit none ; private
9
11
12contains
13
14!> Pointwise evaluation of a polynomial in arbitrary units [A] at x
15!!
16!! The polynomial is defined by the coefficients contained in the
17!! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ...
18!! where C refers to the array 'coeff'.
19!! The number of coefficients is given by ncoef and x
20!! is the coordinate where the polynomial is to be evaluated.
21real function evaluation_polynomial( coeff, ncoef, x )
22 real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial, in units that
23 !! vary with the index k as [A H^(k-1)]
24 integer, intent(in) :: ncoef !< The number of polynomial coefficients
25 real, intent(in) :: x !< The position at which to evaluate the polynomial
26 !! in arbitrary thickness units [H]
27 ! Local variables
28 integer :: k
29 real :: f ! value of polynomial at x in arbitrary units [A]
30
31 f = 0.0
32 do k = 1,ncoef
33 f = f + coeff(k) * ( x**(k-1) )
34 enddo
35
37
38end function evaluation_polynomial
39
40!> Calculates the first derivative of a polynomial evaluated in arbitrary units of [A H-1]
41!! at a point x
42!!
43!! The polynomial is defined by the coefficients contained in the
44!! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ...
45!! where C refers to the array 'coeff'.
46!! The number of coefficients is given by ncoef and x
47!! is the coordinate where the polynomial's derivative is to be evaluated.
48real function first_derivative_polynomial( coeff, ncoef, x )
49 real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial, in units that
50 !! vary with the index k as [A H^(k-1)]
51 integer, intent(in) :: ncoef !< The number of polynomial coefficients
52 real, intent(in) :: x !< The position at which to evaluate the derivative
53 !! in arbitrary thickness units [H]
54 ! Local variables
55 integer :: k
56 real :: f ! value of the derivative at x in [A H-1]
57
58 f = 0.0
59 do k = 2,ncoef
60 f = f + real(k-1)*coeff(k) * ( x**(k-2) )
61 enddo
62
64
66
67!> Exact integration of polynomial of degree npoly in arbitrary units of [A H]
68!!
69!! The array of coefficients (Coeff) must be of size npoly+1.
70real function integration_polynomial( xi0, xi1, Coeff, npoly )
71 real, intent(in) :: xi0 !< The lower bound of the integral in arbitrary
72 !! thickness units [H]
73 real, intent(in) :: xi1 !< The upper bound of the integral in arbitrary
74 !! thickness units [H]
75 real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial, in units that
76 !! vary with the index k as [A H^(k-1)]
77 integer, intent(in) :: npoly !< The degree of the polynomial
78 ! Local variables
79 integer :: k
80 real :: integral ! The integral of the polynomial over the specified range in [A H]
81
82 integral = 0.0
83
84 do k = 1,npoly+1
85 integral = integral + coeff(k) * (xi1**k - xi0**k) / real(k)
86 enddo
87!
88!One non-answer-changing way of unrolling the above is:
89! k=1
90! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k)
91! if (npoly>=1) then
92! k=2
93! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k)
94! endif
95! if (npoly>=2) then
96! k=3
97! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k)
98! endif
99! if (npoly>=3) then
100! k=4
101! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k)
102! endif
103! if (npoly>=4) then
104! k=5
105! integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k)
106! endif
107!
108 integration_polynomial = integral
109
110end function integration_polynomial
111
112!> \namespace polynomial_functions
113!!
114!! Date of creation: 2008.06.12
115!! L. White
116!!
117!! This module contains routines that handle polynomials.
118
119end module polynomial_functions