source: LMDZ5/branches/IPSLCM6.0.8/libf/phylmd/cosp/math_lib.F90 @ 2791

Last change on this file since 2791 was 2435, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2396:2434 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 9.1 KB
Line 
1! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
2! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/math_lib.f90 $
3! MATH_LIB: Mathematics procedures for F90
4! Compiled/Modified:
5!   07/01/06  John Haynes (haynes@atmos.colostate.edu)
6!
7! gamma (function)
8! path_integral (function)
9! avint (subroutine)
10 
11  module math_lib
12  implicit none
13
14  contains
15
16! ----------------------------------------------------------------------------
17! function GAMMA
18! ----------------------------------------------------------------------------
19  function gamma(x)
20  implicit none
21!
22! Purpose:
23!   Returns the gamma function
24!
25! Input:
26!   [x]   value to compute gamma function of
27!
28! Returns:
29!   gamma(x)
30!
31! Coded:
32!   02/02/06  John Haynes (haynes@atmos.colostate.edu)
33!   (original code of unknown origin)
34
35! ----- INPUTS -----
36  real*8, intent(in) :: x
37 
38! ----- OUTPUTS -----
39  real*8 :: gamma
40
41! ----- INTERNAL ----- 
42  real*8 :: pi,ga,z,r,gr
43  real*8 :: g(26)
44  integer :: k,m1,m
45       
46  pi = acos(-1.)   
47  if (x ==int(x)) then
48    if (x > 0.0) then
49      ga=1.0
50      m1=x-1
51      do k=2,m1
52        ga=ga*k
53      enddo
54    else
55      ga=1.0+300
56    endif
57  else
58    if (abs(x) > 1.0) then
59      z=abs(x)
60      m=int(z)
61      r=1.0
62      do k=1,m
63        r=r*(z-k)
64      enddo
65      z=z-m
66    else
67      z=x
68    endif
69    data g/1.0,0.5772156649015329, &
70           -0.6558780715202538, -0.420026350340952d-1, &
71           0.1665386113822915,-.421977345555443d-1, &
72           -.96219715278770d-2, .72189432466630d-2, &
73           -.11651675918591d-2, -.2152416741149d-3, &
74           .1280502823882d-3, -.201348547807d-4, &
75           -.12504934821d-5, .11330272320d-5, &
76           -.2056338417d-6, .61160950d-8, &
77           .50020075d-8, -.11812746d-8, &
78           .1043427d-9, .77823d-11, &
79          -.36968d-11, .51d-12, &
80          -.206d-13, -.54d-14, .14d-14, .1d-15/
81    gr=g(26)
82    do k=25,1,-1
83      gr=gr*z+g(k)
84    enddo
85    ga=1.0/(gr*z)
86    if (abs(x) > 1.0) then
87      ga=ga*r
88      if (x < 0.0) ga=-pi/(x*ga*sin(pi*x))
89    endif
90  endif
91  gamma = ga
92  return
93  end function gamma
94 
95! ----------------------------------------------------------------------------
96! function PATH_INTEGRAL
97! ----------------------------------------------------------------------------
98  function path_integral(f,s,i1,i2)
99  use m_mrgrnk
100  use array_lib
101  implicit none
102!
103! Purpose:
104!   evalues the integral (f ds) between f(index=i1) and f(index=i2)
105!   using the AVINT procedure
106!
107! Inputs:
108!   [f]    functional values
109!   [s]    abscissa values
110!   [i1]   index of lower limit
111!   [i2]   index of upper limit
112!
113! Returns:
114!   result of path integral
115!
116! Notes:
117!   [s] may be in forward or reverse numerical order
118!
119! Requires:
120!   mrgrnk package
121!
122! Created:
123!   02/02/06  John Haynes (haynes@atmos.colostate.edu)
124
125! ----- INPUTS ----- 
126  real*8, intent(in), dimension(:) :: f,s 
127  integer, intent(in) :: i1, i2
128
129! ---- OUTPUTS -----
130  real*8 :: path_integral 
131 
132! ----- INTERNAL -----   
133  real*8 :: sumo, deltah, val
134  integer*4 :: nelm, j
135  integer*4, dimension(i2-i1+1) :: idx
136  real*8, dimension(i2-i1+1) :: f_rev, s_rev
137
138  nelm = i2-i1+1
139  if (nelm > 3) then
140    call mrgrnk(s(i1:i2),idx)
141    s_rev = s(idx)
142    f_rev = f(idx)
143    call avint(f_rev(i1:i2),s_rev(i1:i2),(i2-i1+1), &
144      s_rev(i1),s_rev(i2), val)
145    path_integral = val
146  else
147     sumo = 0.
148     do j=i1,i2
149       deltah = abs(s(i1+1)-s(i1))
150       sumo = sumo + f(j)*deltah
151    enddo
152    path_integral = sumo
153  endif
154  ! print *, sumo
155
156  return
157  end function path_integral
158 
159! ----------------------------------------------------------------------------
160! subroutine AVINT
161! ----------------------------------------------------------------------------
162  subroutine avint ( ftab, xtab, ntab, a_in, b_in, result )
163  implicit none
164!
165! Purpose:
166!   estimate the integral of unevenly spaced data
167!
168! Inputs:
169!   [ftab]     functional values
170!   [xtab]     abscissa values
171!   [ntab]     number of elements of [ftab] and [xtab]
172!   [a]        lower limit of integration
173!   [b]        upper limit of integration
174!
175! Outputs:
176!   [result]   approximate value of integral
177!
178! Reference:
179!   From SLATEC libraries, in public domain
180!
181!***********************************************************************
182!
183!  AVINT estimates the integral of unevenly spaced data.
184!
185!  Discussion:
186!
187!    The method uses overlapping parabolas and smoothing.
188!
189!  Modified:
190!
191!    30 October 2000
192!    4 January 2008, A. Bodas-Salcedo. Error control for XTAB taken out of
193!                    loop to allow vectorization.
194!
195!  Reference:
196!
197!    Philip Davis and Philip Rabinowitz,
198!    Methods of Numerical Integration,
199!    Blaisdell Publishing, 1967.
200!
201!    P E Hennion,
202!    Algorithm 77,
203!    Interpolation, Differentiation and Integration,
204!    Communications of the Association for Computing Machinery,
205!    Volume 5, page 96, 1962.
206!
207!  Parameters:
208!
209!    Input, real ( kind = 8 ) FTAB(NTAB), the function values,
210!    FTAB(I) = F(XTAB(I)).
211!
212!    Input, real ( kind = 8 ) XTAB(NTAB), the abscissas at which the
213!    function values are given.  The XTAB's must be distinct
214!    and in ascending order.
215!
216!    Input, integer NTAB, the number of entries in FTAB and
217!    XTAB.  NTAB must be at least 3.
218!
219!    Input, real ( kind = 8 ) A, the lower limit of integration.  A should
220!    be, but need not be, near one endpoint of the interval
221!    (X(1), X(NTAB)).
222!
223!    Input, real ( kind = 8 ) B, the upper limit of integration.  B should
224!    be, but need not be, near one endpoint of the interval
225!    (X(1), X(NTAB)).
226!
227!    Output, real ( kind = 8 ) RESULT, the approximate value of the integral.
228
229  integer, intent(in) :: ntab
230
231  integer,parameter :: KR8 = selected_real_kind(15,300)
232  real ( kind = KR8 ), intent(in) :: a_in
233  real ( kind = KR8 ) a
234  real ( kind = KR8 ) atemp
235  real ( kind = KR8 ), intent(in) :: b_in
236  real ( kind = KR8 ) b
237  real ( kind = KR8 ) btemp
238  real ( kind = KR8 ) ca
239  real ( kind = KR8 ) cb
240  real ( kind = KR8 ) cc
241  real ( kind = KR8 ) ctemp
242  real ( kind = KR8 ), intent(in) :: ftab(ntab)
243  integer i
244  integer ihi
245  integer ilo
246  integer ind
247  real ( kind = KR8 ), intent(out) :: result
248  real ( kind = KR8 ) sum1
249  real ( kind = KR8 ) syl
250  real ( kind = KR8 ) term1
251  real ( kind = KR8 ) term2
252  real ( kind = KR8 ) term3
253  real ( kind = KR8 ) x1
254  real ( kind = KR8 ) x2
255  real ( kind = KR8 ) x3
256  real ( kind = KR8 ), intent(in) :: xtab(ntab)
257  logical lerror
258 
259  lerror = .false.
260  a = a_in
261  b = b_in 
262 
263  if ( ntab < 3 ) then
264    write ( *, '(a)' ) ' '
265    write ( *, '(a)' ) 'AVINT - Fatal error!'
266    write ( *, '(a,i6)' ) '  NTAB is less than 3.  NTAB = ', ntab
267    stop
268  end if
269 
270  do i = 2, ntab
271    if ( xtab(i) <= xtab(i-1) ) then
272       lerror = .true.
273       exit
274    end if
275  end do
276 
277  if (lerror) then
278      write ( *, '(a)' ) ' '
279      write ( *, '(a)' ) 'AVINT - Fatal error!'
280      write ( *, '(a)' ) '  XTAB(I) is not greater than XTAB(I-1).'
281      write ( *, '(a,i6)' ) '  Here, I = ', i
282      write ( *, '(a,g14.6)' ) '  XTAB(I-1) = ', xtab(i-1)
283      write ( *, '(a,g14.6)' ) '  XTAB(I) =   ', xtab(i)
284      stop 
285  end if
286 
287  result = 0.0D+00
288 
289  if ( a == b ) then
290    write ( *, '(a)' ) ' '
291    write ( *, '(a)' ) 'AVINT - Warning!'
292    write ( *, '(a)' ) '  A = B, integral=0.'
293    return
294  end if
295!
296!  If B < A, temporarily switch A and B, and store sign.
297!
298  if ( b < a ) then
299    syl = b
300    b = a
301    a = syl
302    ind = -1
303  else
304    syl = a
305    ind = 1
306  end if
307!
308!  Bracket A and B between XTAB(ILO) and XTAB(IHI).
309!
310  ilo = 1
311  ihi = ntab
312
313  do i = 1, ntab
314    if ( a <= xtab(i) ) then
315      exit
316    end if
317    ilo = ilo + 1
318  end do
319
320  ilo = max ( 2, ilo )
321  ilo = min ( ilo, ntab - 1 )
322
323  do i = 1, ntab
324    if ( xtab(i) <= b ) then
325      exit
326    end if
327    ihi = ihi - 1
328  end do
329 
330  ihi = min ( ihi, ntab - 1 )
331  ihi = max ( ilo, ihi - 1 )
332!
333!  Carry out approximate integration from XTAB(ILO) to XTAB(IHI).
334!
335  sum1 = 0.0D+00
336 
337  do i = ilo, ihi
338 
339    x1 = xtab(i-1)
340    x2 = xtab(i)
341    x3 = xtab(i+1)
342   
343    term1 = ftab(i-1) / ( ( x1 - x2 ) * ( x1 - x3 ) )
344    term2 = ftab(i)   / ( ( x2 - x1 ) * ( x2 - x3 ) )
345    term3 = ftab(i+1) / ( ( x3 - x1 ) * ( x3 - x2 ) )
346 
347    atemp = term1 + term2 + term3
348
349    btemp = - ( x2 + x3 ) * term1 &
350            - ( x1 + x3 ) * term2 &
351            - ( x1 + x2 ) * term3
352
353    ctemp = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3
354 
355    if ( i <= ilo ) then
356      ca = atemp
357      cb = btemp
358      cc = ctemp
359    else
360      ca = 0.5D+00 * ( atemp + ca )
361      cb = 0.5D+00 * ( btemp + cb )
362      cc = 0.5D+00 * ( ctemp + cc )
363    end if
364 
365    sum1 = sum1 &
366          + ca * ( x2**3 - syl**3 ) / 3.0D+00 &
367          + cb * 0.5D+00 * ( x2**2 - syl**2 ) &
368          + cc * ( x2 - syl )
369 
370    ca = atemp
371    cb = btemp
372    cc = ctemp
373 
374    syl = x2
375 
376  end do
377 
378  result = sum1 &
379        + ca * ( b**3 - syl**3 ) / 3.0D+00 &
380        + cb * 0.5D+00 * ( b**2 - syl**2 ) &
381        + cc * ( b - syl )
382!
383!  Restore original values of A and B, reverse sign of integral
384!  because of earlier switch.
385!
386  if ( ind /= 1 ) then
387    ind = 1
388    syl = b
389    b = a
390    a = syl
391    result = -result
392  end if
393 
394  return
395  end subroutine avint
396 
397  end module math_lib
Note: See TracBrowser for help on using the repository browser.