source: LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/math_lib.F90 @ 1329

Last change on this file since 1329 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

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