A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
Functions/Subroutines
int_comp Module Reference

Utility module containing the routines to perform the integration of functions. More...

Functions/Subroutines

subroutine, public integrate (func, ss)
 Routine to compute integrals of function from O to #maxint. More...
 
subroutine qromb (func, a, b, ss)
 Romberg integration routine. More...
 
subroutine qromo (func, a, b, ss, choose)
 Romberg integration routine on an open interval. More...
 
subroutine polint (xa, ya, n, x, y, dy)
 Polynomial interpolation routine. More...
 
subroutine trapzd (func, a, b, s, n)
 Trapezoidal rule integration routine. More...
 
subroutine midpnt (func, a, b, s, n)
 Midpoint rule integration routine. More...
 
subroutine midexp (funk, aa, bb, s, n)
 Midpoint routine for bb infinite with funk decreasing infinitely rapidly at infinity. More...
 

Detailed Description

Utility module containing the routines to perform the integration of functions.

Remarks
Most are taken from the Numerical Recipes

Function/Subroutine Documentation

subroutine, public int_comp::integrate ( external  func,
real(kind=8)  ss 
)

Routine to compute integrals of function from O to #maxint.

Parameters
funcfunction to integrate
ssresult of the integration

Definition at line 30 of file int_comp.f90.

30  REAL(KIND=8) :: ss,func,b
31  EXTERNAL func
32  b=maxint
33  ! CALL qromo(func,0.D0,1.D0,ss,midexp)
34  CALL qromb(func,0.d0,b,ss)
subroutine int_comp::midexp ( external  funk,
real(kind=8)  aa,
real(kind=8)  bb,
real(kind=8)  s,
integer  n 
)
private

Midpoint routine for bb infinite with funk decreasing infinitely rapidly at infinity.

Parameters
funkfunction to integrate
aalower limit of the integral
bbhigher limit of the integral
sresult of the integration
nhigher stage of the rule to be computed

Definition at line 200 of file int_comp.f90.

200  INTEGER :: n
201  REAL(KIND=8) :: aa,bb,s,funk
202  EXTERNAL funk
203  INTEGER :: it,j
204  REAL(KIND=8) :: ddel,del,sum,tnm,x,func,a,b
205  func(x)=funk(-log(x))/x
206  b=exp(-aa)
207  a=0.
208  if (n.eq.1) then
209  s=(b-a)*func(0.5*(a+b))
210  else
211  it=3**(n-2)
212  tnm=it
213  del=(b-a)/(3.*tnm)
214  ddel=del+del
215  x=a+0.5*del
216  sum=0.
217  do j=1,it
218  sum=sum+func(x)
219  x=x+ddel
220  sum=sum+func(x)
221  x=x+del
222  end do
223  s=(s+(b-a)*sum/tnm)/3.
224  endif
225  return
subroutine int_comp::midpnt ( external  func,
real(kind=8)  a,
real(kind=8)  b,
real(kind=8)  s,
integer  n 
)
private

Midpoint rule integration routine.

Parameters
funcfunction to integrate
alower limit of the integral
bhigher limit of the integral
sresult of the integration
nhigher stage of the rule to be computed

Definition at line 167 of file int_comp.f90.

167  INTEGER :: n
168  REAL(KIND=8) :: a,b,s,func
169  EXTERNAL func
170  INTEGER :: it,j
171  REAL(KIND=8) :: ddel,del,sum,tnm,x
172  if (n.eq.1) then
173  s=(b-a)*func(0.5*(a+b))
174  else
175  it=3**(n-2)
176  tnm=it
177  del=(b-a)/(3.*tnm)
178  ddel=del+del
179  x=a+0.5*del
180  sum=0.
181  do j=1,it
182  sum=sum+func(x)
183  x=x+ddel
184  sum=sum+func(x)
185  x=x+del
186  end do
187  s=(s+(b-a)*sum/tnm)/3.
188  endif
189  return
subroutine int_comp::polint ( real(kind=8), dimension(n)  xa,
real(kind=8), dimension(n)  ya,
integer  n,
real(kind=8)  x,
real(kind=8)  y,
real(kind=8)  dy 
)
private

Polynomial interpolation routine.

Definition at line 91 of file int_comp.f90.

91  INTEGER :: n,nmax
92  REAL(KIND=8) :: dy,x,y,xa(n),ya(n)
93  parameter(nmax=10)
94  INTEGER :: i,m,ns
95  REAL(KIND=8) :: den,dif,dift,ho,hp,w,c(nmax),d(nmax)
96  ns=1
97  dif=abs(x-xa(1))
98  do i=1,n
99  dift=abs(x-xa(i))
100  if (dift.lt.dif) then
101  ns=i
102  dif=dift
103  endif
104  c(i)=ya(i)
105  d(i)=ya(i)
106  end do
107  y=ya(ns)
108  ns=ns-1
109  do m=1,n-1
110  do i=1,n-m
111  ho=xa(i)-x
112  hp=xa(i+m)-x
113  w=c(i+1)-d(i)
114  den=ho-hp
115  if(den.eq.0.)stop 'failure in polint'
116  den=w/den
117  d(i)=hp*den
118  c(i)=ho*den
119  end do
120  if (2*ns.lt.n-m)then
121  dy=c(ns+1)
122  else
123  dy=d(ns)
124  ns=ns-1
125  endif
126  y=y+dy
127  end do
128  return
subroutine int_comp::qromb ( external  func,
real(kind=8)  a,
real(kind=8)  b,
real(kind=8)  ss 
)
private

Romberg integration routine.

Parameters
funcfunction to integrate
alower limit of the integral
bhigher limit of the integral
funcfunction to integrate
ssresult of the integration

Definition at line 44 of file int_comp.f90.

44  INTEGER :: jmax,jmaxp,k,km
45  REAL(KIND=8) :: a,b,func,ss,eps
46  EXTERNAL func
47  parameter(eps=1.d-6, jmax=20, jmaxp=jmax+1, k=5, km=k-1)
48  INTEGER j
49  REAL(KIND=8) :: dss,h(jmaxp),s(jmaxp)
50  h(1)=1.
51  DO j=1,jmax
52  CALL trapzd(func,a,b,s(j),j)
53  IF (j.ge.k) THEN
54  CALL polint(h(j-km),s(j-km),k,0.d0,ss,dss)
55  IF (abs(dss).le.eps*abs(ss)) RETURN
56  ENDIF
57  s(j+1)=s(j)
58  h(j+1)=0.25*h(j)
59  ENDDO
60  stop 'too many steps in qromb'
subroutine int_comp::qromo ( external  func,
real(kind=8)  a,
real(kind=8)  b,
real(kind=8)  ss,
external  choose 
)
private

Romberg integration routine on an open interval.

Parameters
alower limit of the integral
bhigher limit of the integral
funcfunction to integrate
ssresult of the integration
choseroutine to perform the integration

Definition at line 70 of file int_comp.f90.

70  INTEGER :: jmax,jmaxp,k,km
71  REAL(KIND=8) :: a,b,func,ss,eps
72  EXTERNAL func,choose
73  parameter(eps=1.e-6, jmax=14, jmaxp=jmax+1, k=5, km=k-1)
74  INTEGER :: j
75  REAL(KIND=8) :: dss,h(jmaxp),s(jmaxp)
76  h(1)=1.
77  DO j=1,jmax
78  CALL choose(func,a,b,s(j),j)
79  IF (j.ge.k) THEN
80  call polint(h(j-km),s(j-km),k,0.d0,ss,dss)
81  if (abs(dss).le.eps*abs(ss)) return
82  ENDIF
83  s(j+1)=s(j)
84  h(j+1)=h(j)/9.
85  ENDDO
86  stop 'too many steps in qromo'
subroutine int_comp::trapzd ( external  func,
real(kind=8)  a,
real(kind=8)  b,
real(kind=8)  s,
integer  n 
)
private

Trapezoidal rule integration routine.

Parameters
funcfunction to integrate
alower limit of the integral
bhigher limit of the integral
sresult of the integration
nhigher stage of the rule to be computed

Definition at line 138 of file int_comp.f90.

138  INTEGER :: n
139  REAL(KIND=8) :: a,b,s,func
140  EXTERNAL func
141  INTEGER :: it,j
142  REAL(KIND=8) :: del,sum,tnm,x
143  if (n.eq.1) then
144  s=0.5*(b-a)*(func(a)+func(b))
145  else
146  it=2**(n-2)
147  tnm=it
148  del=(b-a)/tnm
149  x=a+0.5*del
150  sum=0.
151  do j=1,it
152  sum=sum+func(x)
153  x=x+del
154  end do
155  s=0.5*(s+(b-a)*sum/tnm)
156  endif
157  return