A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
int_comp.f90
Go to the documentation of this file.
1 
2 ! int_comp.f90
3 !
4 !> Utility module containing the routines to perform the integration of functions
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !> Most are taken from the Numerical Recipes
14 !
15 !---------------------------------------------------------------------------!
16 
17 MODULE int_comp
18  USE stoch_params, only: maxint
19 
20  PRIVATE
21 
22  PUBLIC :: integrate
23 
24 
25 CONTAINS
26  !> Routine to compute integrals of function from O to #maxint
27  !> @param func function to integrate
28  !> @param ss result of the integration
29  SUBROUTINE integrate(func,ss)
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)
35  END SUBROUTINE integrate
36 
37  !> Romberg integration routine
38  !> @param func function to integrate
39  !> @param a lower limit of the integral
40  !> @param b higher limit of the integral
41  !> @param func function to integrate
42  !> @param ss result of the integration
43  SUBROUTINE qromb(func,a,b,ss)
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'
61  END SUBROUTINE qromb
62 
63  !> Romberg integration routine on an open interval
64  !> @param a lower limit of the integral
65  !> @param b higher limit of the integral
66  !> @param func function to integrate
67  !> @param ss result of the integration
68  !> @param chose routine to perform the integration
69  SUBROUTINE qromo(func,a,b,ss,choose)
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'
87  END SUBROUTINE qromo
88 
89  !> Polynomial interpolation routine
90  SUBROUTINE polint(xa,ya,n,x,y,dy)
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
129  END SUBROUTINE polint
130 
131  !> Trapezoidal rule integration routine
132  !> @param func function to integrate
133  !> @param a lower limit of the integral
134  !> @param b higher limit of the integral
135  !> @param s result of the integration
136  !> @param n higher stage of the rule to be computed
137  SUBROUTINE trapzd(func,a,b,s,n)
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
158  END SUBROUTINE trapzd
159 
160  !> Midpoint rule integration routine
161  !> @param func function to integrate
162  !> @param a lower limit of the integral
163  !> @param b higher limit of the integral
164  !> @param s result of the integration
165  !> @param n higher stage of the rule to be computed
166  SUBROUTINE midpnt(func,a,b,s,n)
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
190  END SUBROUTINE midpnt
191 
192  !> Midpoint routine for bb infinite with funk decreasing infinitely rapidly at
193  !> infinity
194  !> @param funk function to integrate
195  !> @param aa lower limit of the integral
196  !> @param bb higher limit of the integral
197  !> @param s result of the integration
198  !> @param n higher stage of the rule to be computed
199  SUBROUTINE midexp(funk,aa,bb,s,n)
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
226  END SUBROUTINE midexp
227 END MODULE int_comp
The stochastic models parameters module.
subroutine qromo(func, a, b, ss, choose)
Romberg integration routine on an open interval.
Definition: int_comp.f90:70
Utility module containing the routines to perform the integration of functions.
Definition: int_comp.f90:17
subroutine midpnt(func, a, b, s, n)
Midpoint rule integration routine.
Definition: int_comp.f90:167
subroutine polint(xa, ya, n, x, y, dy)
Polynomial interpolation routine.
Definition: int_comp.f90:91
subroutine trapzd(func, a, b, s, n)
Trapezoidal rule integration routine.
Definition: int_comp.f90:138
subroutine, public integrate(func, ss)
Routine to compute integrals of function from O to #maxint.
Definition: int_comp.f90:30
real(kind=8) maxint
Upper integration limit of the correlations.
subroutine qromb(func, a, b, ss)
Romberg integration routine.
Definition: int_comp.f90:44
subroutine midexp(funk, aa, bb, s, n)
Midpoint routine for bb infinite with funk decreasing infinitely rapidly at infinity.
Definition: int_comp.f90:200