source: LMDZ6/trunk/libf/phylmd/cospv2/math_lib.f90 @ 5451

Last change on this file since 5451 was 5268, checked in by abarral, 2 months ago

.f90 <-> .F90 depending on cpp key use
File size: 12.5 KB
RevLine 
[3491]1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2015, Regents of the University of Colorado
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History:
30! July 2006: John Haynes      - Initial version
31! May 2015:  Dustin Swales    - Modified for COSPv2.0
32!
33! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34module math_lib
35  USE COSP_KINDS,     ONLY: wp
36  use mod_cosp_error, ONLY: errorMessage
37  implicit none
38
39contains
40  ! ########################################################################## 
41  !                           function PATH_INTEGRAL
42  ! ########################################################################## 
43  function path_integral(f,s,i1,i2)
44    use m_mrgrnk
45    use array_lib
46    implicit none
47    ! ########################################################################
48    ! Purpose:
49    !   evalues the integral (f ds) between f(index=i1) and f(index=i2)
50    !   using the AVINT procedure
51    !
52    ! Inputs:
53    !   [f]    functional values
54    !   [s]    abscissa values
55    !   [i1]   index of lower limit
56    !   [i2]   index of upper limit
57    !
58    ! Returns:
59    !   result of path integral
60    !
61    ! Notes:
62    !   [s] may be in forward or reverse numerical order
63    !
64    ! Requires:
65    !   mrgrnk package
66    !
67    ! Created:
68    !   02/02/06  John Haynes (haynes@atmos.colostate.edu)
69    ! ########################################################################
70
71    ! INPUTS
72    real(wp),intent(in), dimension(:) :: &
73         f,  & ! Functional values
74         s     ! Abscissa values 
75    integer, intent(in) :: &
76         i1, & ! Index of lower limit
77         i2    ! Index of upper limit
78
79    ! OUTPUTS
80    real(wp) :: path_integral 
81 
82    ! Internal variables
83    real(wp) :: sumo, deltah, val
84    integer :: nelm, j
85    integer, dimension(i2-i1+1) :: idx
86    real(wp), dimension(i2-i1+1) :: f_rev, s_rev
87
88    nelm = i2-i1+1
89    if (nelm > 3) then
90       call mrgrnk(s(i1:i2),idx)
91       s_rev = s(idx)
92       f_rev = f(idx)
93       call avint(f_rev(i1:i2),s_rev(i1:i2),(i2-i1+1), &
94            s_rev(i1),s_rev(i2), val)
95       path_integral = val
96    else
97       sumo = 0._wp
98       do j=i1,i2
99          deltah = abs(s(i1+1)-s(i1))
100          sumo = sumo + f(j)*deltah
101       enddo
102       path_integral = sumo
103    endif
104   
105    return
106  end function path_integral
107 
108  ! ##########################################################################
109  !                            subroutine AVINT
110  ! ##########################################################################
111  subroutine avint ( ftab, xtab, ntab, a_in, b_in, result )
112    implicit none
113    ! ######################################################################## 
114    ! Purpose:
115    !   estimate the integral of unevenly spaced data
116    !
117    ! Inputs:
118    !   [ftab]     functional values
119    !   [xtab]     abscissa values
120    !   [ntab]     number of elements of [ftab] and [xtab]
121    !   [a]        lower limit of integration
122    !   [b]        upper limit of integration
123    !
124    ! Outputs:
125    !   [result]   approximate value of integral
126    !
127    ! Reference:
128    !   From SLATEC libraries, in public domain
129    !
130    !***********************************************************************
131    !
132    !  AVINT estimates the integral of unevenly spaced data.
133    !
134    !  Discussion:
135    !
136    !    The method uses overlapping parabolas and smoothing.
137    !
138    !  Modified:
139    !
140    !    30 October 2000
141    !    4 January 2008, A. Bodas-Salcedo. Error control for XTAB taken out of
142    !                    loop to allow vectorization.
143    !
144    !  Reference:
145    !
146    !    Philip Davis and Philip Rabinowitz,
147    !    Methods of Numerical Integration,
148    !    Blaisdell Publishing, 1967.
149    !
150    !    P E Hennion,
151    !    Algorithm 77,
152    !    Interpolation, Differentiation and Integration,
153    !    Communications of the Association for Computing Machinery,
154    !    Volume 5, page 96, 1962.
155    !
156    !  Parameters:
157    !
158    !    Input, real ( kind = 8 ) FTAB(NTAB), the function values,
159    !    FTAB(I) = F(XTAB(I)).
160    !
161    !    Input, real ( kind = 8 ) XTAB(NTAB), the abscissas at which the
162    !    function values are given.  The XTAB's must be distinct
163    !    and in ascending order.
164    !
165    !    Input, integer NTAB, the number of entries in FTAB and
166    !    XTAB.  NTAB must be at least 3.
167    !
168    !    Input, real ( kind = 8 ) A, the lower limit of integration.  A should
169    !    be, but need not be, near one endpoint of the interval
170    !    (X(1), X(NTAB)).
171    !
172    !    Input, real ( kind = 8 ) B, the upper limit of integration.  B should
173    !    be, but need not be, near one endpoint of the interval
174    !    (X(1), X(NTAB)).
175    !
176    !    Output, real ( kind = 8 ) RESULT, the approximate value of the integral.
177    ! ########################################################################## 
178
179    ! INPUTS
180    integer,intent(in) :: &
181         ntab    ! Number of elements of [ftab] and [xtab]
182    real(wp),intent(in) :: &
183         a_in, & ! Lower limit of integration
184         b_in    ! Upper limit of integration
185    real(wp),intent(in),dimension(ntab) :: &
186         ftab, & ! Functional values
187         xtab    ! Abscissa value
188   
189    ! OUTPUTS
190    real(wp),intent(out) :: result  ! Approximate value of integral
191
192    ! Internal varaibles
193    real(wp) :: a, atemp, b, btemp,ca,cb,cc,ctemp,sum1,syl,term1,term2,term3,x1,x2,x3
194    integer  :: i,ihi,ilo,ind
195    logical  :: lerror
196 
197    lerror = .false.
198    a = a_in
199    b = b_in 
200 
201    if ( ntab < 3 ) then
202       call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Ntab is less than 3.')
203       return
204    end if
205   
206    do i = 2, ntab
207       if ( xtab(i) <= xtab(i-1) ) then
208          lerror = .true.
209          exit
210       end if
211    end do
212   
213    if (lerror) then
214       call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Xtab(i) is not greater than Xtab(i-1).')
215       return
216    end if
217   
218!ds    result = 0.0D+00
219    result = 0._wp
220   
221    if ( a == b ) then
222       call errorMessage('WARNING(optics/math_lib.f90:AVINT): A=B => integral=0')
223       return
224    end if
225   
226    !  If B < A, temporarily switch A and B, and store sign.
227    if ( b < a ) then
228       syl = b
229       b = a
230       a = syl
231       ind = -1
232    else
233       syl = a
234       ind = 1
235    end if
236   
237    !  Bracket A and B between XTAB(ILO) and XTAB(IHI).
238    ilo = 1
239    ihi = ntab
240   
241    do i = 1, ntab
242       if ( a <= xtab(i) ) then
243          exit
244       end if
245       ilo = ilo + 1
246    end do
247   
248    ilo = max ( 2, ilo )
249    ilo = min ( ilo, ntab - 1 )
250   
251    do i = 1, ntab
252       if ( xtab(i) <= b ) then
253          exit
254       end if
255       ihi = ihi - 1
256    end do
257   
258    ihi = min ( ihi, ntab - 1 )
259    ihi = max ( ilo, ihi - 1 )
260   
261    !  Carry out approximate integration from XTAB(ILO) to XTAB(IHI).
262    sum1 = 0._wp
263!ds    sum1 = 0.0D+00
264   
265    do i = ilo, ihi
266       
267       x1 = xtab(i-1)
268       x2 = xtab(i)
269       x3 = xtab(i+1)
270       
271       term1 = ftab(i-1) / ( ( x1 - x2 ) * ( x1 - x3 ) )
272       term2 = ftab(i)   / ( ( x2 - x1 ) * ( x2 - x3 ) )
273       term3 = ftab(i+1) / ( ( x3 - x1 ) * ( x3 - x2 ) )
274       
275       atemp = term1 + term2 + term3
276       
277       btemp = - ( x2 + x3 ) * term1 &
278            - ( x1 + x3 ) * term2 &
279            - ( x1 + x2 ) * term3
280       
281       ctemp = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3
282       
283       if ( i <= ilo ) then
284          ca = atemp
285          cb = btemp
286          cc = ctemp
287       else
288          ca = 0.5_wp * ( atemp + ca )
289          cb = 0.5_wp * ( btemp + cb )
290          cc = 0.5_wp * ( ctemp + cc )
291!ds          ca = 0.5D+00 * ( atemp + ca )
292!ds          cb = 0.5D+00 * ( btemp + cb )
293!ds          cc = 0.5D+00 * ( ctemp + cc )
294       end if
295       
296       sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3._wp &
297            + cb * 0.5_wp * ( x2**2 - syl**2 ) + cc * ( x2 - syl )
298!ds       sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3.0D+00 &
299!ds            + cb * 0.5D+00 * ( x2**2 - syl**2 ) + cc * ( x2 - syl )
300       
301       ca = atemp
302       cb = btemp
303       cc = ctemp
304       
305       syl = x2
306       
307    end do
308
309    result = sum1 + ca * ( b**3 - syl**3 ) / 3._wp &
310         + cb * 0.5_wp * ( b**2 - syl**2 ) + cc * ( b - syl )
311!ds    result = sum1 + ca * ( b**3 - syl**3 ) / 3.0D+00 &
312!ds         + cb * 0.5D+00 * ( b**2 - syl**2 ) + cc * ( b - syl )
313
314    !  Restore original values of A and B, reverse sign of integral
315    !  because of earlier switch.
316    if ( ind /= 1 ) then
317       ind = 1
318       syl = b
319       b = a
320       a = syl
321       result = -result
322    end if
323   
324    return
325  end subroutine avint
326  ! ######################################################################################
327  ! SUBROUTINE gamma
328  ! Purpose:
329  !   Returns the gamma function
330  !
331  ! Input:
332  !   [x]   value to compute gamma function of
333  !
334  ! Returns:
335  !   gamma(x)
336  !
337  ! Coded:
338  !   02/02/06  John Haynes (haynes@atmos.colostate.edu)
339  !   (original code of unknown origin)
340  ! ######################################################################################
341  function gamma(x)
342    ! Inputs
343    real(wp), intent(in) :: x
344
345    ! Outputs
346    real(wp) :: gamma
347   
348    ! Local variables
349    real(wp) :: pi,ga,z,r,gr
350    integer  :: k,m1,m 
351   
352    ! Parameters
353    real(wp),dimension(26),parameter :: &
354         g = (/1.0,0.5772156649015329, -0.6558780715202538, -0.420026350340952e-1,     &
355               0.1665386113822915,-0.421977345555443e-1,-0.96219715278770e-2,              &
356               0.72189432466630e-2,-0.11651675918591e-2, -0.2152416741149e-3,                &
357               0.1280502823882e-3, -0.201348547807e-4, -0.12504934821e-5, 0.11330272320e-5,  &
358               -0.2056338417e-6, 0.61160950e-8,0.50020075e-8, -0.11812746e-8, 0.1043427e-9,   &
359               0.77823e-11, -0.36968e-11, 0.51e-12, -0.206e-13, -0.54e-14, 0.14e-14, 0.1e-15/) 
360!ds    real(wp),dimension(26),parameter :: &
361!ds         g = (/1.0d0,0.5772156649015329d0, -0.6558780715202538d0, -0.420026350340952d-1,     &
362!ds               0.1665386113822915d0,-0.421977345555443d-1,-0.96219715278770d-2,              &
363!ds               0.72189432466630d-2,-0.11651675918591d-2, -0.2152416741149d-3,                &
364!ds               0.1280502823882d-3, -0.201348547807d-4, -0.12504934821d-5, 0.11330272320d-5,  &
365!ds               -0.2056338417d-6, 0.61160950d-8,0.50020075d-8, -0.11812746d-8, 0.1043427d-9,   &
366!ds               0.77823d-11, -0.36968d-11, 0.51d-12, -0.206d-13, -0.54d-14, 0.14d-14, 0.1d-15/) 
367   
368    pi = acos(-1._wp)   
369    if (x ==int(x)) then
370       if (x > 0.0) then
371          ga=1._wp
372          m1=x-1
373          do k=2,m1
374             ga=ga*k
375          enddo
376       else
377          ga=1._wp+300
378       endif
379    else
380       if (abs(x) > 1.0) then
381          z=abs(x)
382          m=int(z)
383          r=1._wp
384          do k=1,m
385             r=r*(z-k)
386          enddo
387          z=z-m
388       else
389          z=x
390       endif
391       gr=g(26)
392       do k=25,1,-1
393          gr=gr*z+g(k)
394       enddo
395       ga=1._wp/(gr*z)
396       if (abs(x) > 1.0) then
397          ga=ga*r
398          if (x < 0.0) ga=-pi/(x*ga*sin(pi*x))
399       endif
400    endif
401    gamma = ga
402    return
403  end function gamma
404end module math_lib
Note: See TracBrowser for help on using the repository browser.