A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
corrmod.f90
Go to the documentation of this file.
1 
2 ! corrmod.f90
3 !
4 !> Module to initialize the correlation matrix of the unresolved variables
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !
14 !---------------------------------------------------------------------------!
15 
16 MODULE corrmod
17  USE params, only: ndim
19  USE sf_def, only: n_unres,ind
20  USE util, only: invmat
21 
22  PRIVATE
23 
24  PUBLIC :: init_corr
25 
26  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: mean !< Vector holding the mean of the unresolved dynamics (reduced version)
27  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: mean_full !< Vector holding the mean of the unresolved dynamics (full version)
28  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: corr_i_full !< Covariance matrix of the unresolved variables (full version)
29  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: inv_corr_i_full !< Inverse of the covariance matrix of the unresolved variables (full version)
30  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: corr_i !< Covariance matrix of the unresolved variables (reduced version)
31  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: inv_corr_i !< Inverse of the covariance matrix of the unresolved variables (reduced version)
32  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: corr_ij !< Matrix holding the correlation matrix at a given time
33  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: y2 !< Vector holding coefficient of the spline and exponential correlation representation
34  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: ya !< Vector holding coefficient of the spline and exponential correlation representation
35  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: xa !< Vector holding coefficient of the spline and exponential correlation representation
36 
37  !> Integers needed by the spline representation of the correlation
38  INTEGER :: nspl,klo,khi
39 
40  !> Pointer to the correlation computation routine
41  PROCEDURE(corrcomp_from_spline), POINTER, PUBLIC :: corrcomp
42 
43 CONTAINS
44  !> Subroutine to initialise the computation of the correlation
45  SUBROUTINE init_corr
46  INTEGER :: AllocStat,i,j,k,nf
47  REAL(KIND=8), DIMENSION(5) :: dumb
48  LOGICAL :: ex
49 
50  ! Selection of the loading mode
51  SELECT CASE (load_mode)
52  CASE ('defi')
54  CASE ('spli')
55  INQUIRE(file='corrspline.def',exist=ex)
56  IF (.not.ex) stop "*** File corrspline.def not found ! ***"
57  OPEN(20,file='corrspline.def',status='old')
58  READ(20,*) nf,nspl
59  IF (nf /= n_unres) stop "*** Dimension in files corrspline.def and sf.nml do not correspond ! ***"
60  ALLOCATE(xa(nspl), ya(n_unres,n_unres,nspl), y2(n_unres,n_unres,nspl), stat=allocstat)
61  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
62  READ(20,*) xa
63  maxint=xa(nspl)/2
64  DO k=1,n_unres*n_unres
65  READ(20,*) i,j
66  READ(20,*) ya(i,j,:)
67  READ(20,*) y2(i,j,:)
68  ENDDO
69  CLOSE(20)
71  klo=1
72  khi=nspl
73  CASE ('expo')
74  INQUIRE(file='correxpo.def',exist=ex)
75  IF (.not.ex) stop "*** File correxpo.def not found ! ***"
76  OPEN(20,file='correxpo.def',status='old')
77  READ(20,*) nf,maxint
78  IF (nf /= n_unres) stop "*** Dimension in files correxpo.def and sf.nml do not correspond ! ***"
79  ALLOCATE(ya(n_unres,n_unres,5), stat=allocstat)
80  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
81  DO k=1,n_unres*n_unres
82  READ(20,*) i,j,dumb
83  ya(i,j,:)=dumb
84  ENDDO
85  CLOSE(20)
87  CASE DEFAULT
88  stop '*** LOAD_MODE variable not properly defined in corrmod.nml ***'
89  END SELECT
90 
91  ALLOCATE(mean(n_unres),mean_full(0:ndim), stat=allocstat)
92  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
93 
94  ALLOCATE(inv_corr_i(n_unres,n_unres), stat=allocstat)
95  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
96 
97  ALLOCATE(corr_i(n_unres,n_unres), stat=allocstat)
98  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
99 
100  ALLOCATE(corr_ij(n_unres,n_unres), stat=allocstat)
101  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
102 
103  ALLOCATE(corr_i_full(ndim,ndim), stat=allocstat)
104  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
105 
106  ALLOCATE(inv_corr_i_full(ndim,ndim), stat=allocstat)
107  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
108 
109  corr_ij=0.d0
110 
111  CALL corrcomp(0.d0)
114 
115  corr_i_full=0.d0
116  DO i=1,n_unres
117  DO j=1,n_unres
118  corr_i_full(ind(i),ind(j))=corr_i(i,j)
119  ENDDO
120  ENDDO
121 
122  inv_corr_i_full=0.d0
123  DO i=1,n_unres
124  DO j=1,n_unres
125  inv_corr_i_full(ind(i),ind(j))=inv_corr_i(i,j)
126  ENDDO
127  ENDDO
128 
129  mean=0.d0
130  INQUIRE(file='mean.def',exist=ex)
131  IF (ex) THEN
132  OPEN(20,file='mean.def',status='old')
133  READ(20,*) mean
134  CLOSE(20)
135  ENDIF
136 
137  mean_full=0.d0
138  DO i=1,n_unres
139  mean_full(ind(i))=mean(i)
140  ENDDO
141 
142  END SUBROUTINE init_corr
143 
144  !> Subroutine to compute the correlation of the unresolved variables \f$\langle Y \otimes Y^s \rangle\f$
145  !> at time \f$s\f$ from the definition given inside the module
146  !> @param s time \f$s\f$ at which the correlation is computed
147  SUBROUTINE corrcomp_from_def(s)
148  REAL(KIND=8), INTENT(IN) :: s
149  REAL(KIND=8) :: y
150  INTEGER :: i,j
151 
152  ! Definition of the corr_ij matrix as a function of time
153 
154  corr_ij(10,10)=((7.66977252307669*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (1.0240906173830213*cos(&
155  &0.07283568782600224*s))/exp(0.017262015588746404*s) - (0.6434985372062336*sin(0.039597160512071454*s&
156  &))/exp(0.06567483898489704*s) + (0.6434985372062335*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
157  corr_ij(10,9)=((0.6434985372062321*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) - (0.6434985372062324*co&
158  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) + (7.669772523076694*sin(0.039597160512071454*&
159  &s))/exp(0.06567483898489704*s) + (1.024090617383021*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
160  corr_ij(10,8)=0
161  corr_ij(10,7)=0
162  corr_ij(10,6)=0
163  corr_ij(10,5)=((-2.236364132659011*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (6.952804148086198*cos&
164  &(0.07283568782600224*s))/exp(0.017262015588746404*s) - (1.4494534432272481*sin(0.039597160512071454*&
165  &s))/exp(0.06567483898489704*s) - (0.6818177416446283*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
166  corr_ij(10,4)=((1.4494534432272483*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (0.6818177416446293*co&
167  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) - (2.2363641326590127*sin(0.039597160512071454&
168  &*s))/exp(0.06567483898489704*s) + (6.952804148086195*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
169  corr_ij(10,3)=0
170  corr_ij(10,2)=0
171  corr_ij(10,1)=0
172  corr_ij(9,10)=((-0.6434985372062344*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (0.643498537206234*co&
173  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) - (7.669772523076689*sin(0.039597160512071454*&
174  &s))/exp(0.06567483898489704*s) - (1.0240906173830204*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
175  corr_ij(9,9)=((7.66977252307669*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (1.0240906173830204*cos(&
176  &0.07283568782600224*s))/exp(0.017262015588746404*s) - (0.643498537206233*sin(0.039597160512071454*s)&
177  &)/exp(0.06567483898489704*s) + (0.6434985372062327*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
178  corr_ij(9,8)=0
179  corr_ij(9,7)=0
180  corr_ij(9,6)=0
181  corr_ij(9,5)=((-1.4494534432272477*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) - (0.6818177416446249*c&
182  &os(0.07283568782600224*s))/exp(0.017262015588746404*s) + (2.2363641326590105*sin(0.03959716051207145&
183  &4*s))/exp(0.06567483898489704*s) - (6.952804148086195*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
184  corr_ij(9,4)=((-2.2363641326590127*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (6.952804148086194*co&
185  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) - (1.4494534432272486*sin(0.039597160512071454&
186  &*s))/exp(0.06567483898489704*s) - (0.6818177416446249*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
187  corr_ij(9,3)=0
188  corr_ij(9,2)=0
189  corr_ij(9,1)=0
190  corr_ij(8,10)=0
191 
192  corr_ij(8,9)=0
193  corr_ij(8,8)=(9.135647293470983/exp(0.05076718239027029*s) + 2.2233889637758932/exp(0.016285467000648854*s))
194  corr_ij(8,7)=0
195  corr_ij(8,6)=0
196  corr_ij(8,5)=0
197  corr_ij(8,4)=0
198  corr_ij(8,3)=(-5.938566084855411/exp(0.05076718239027029*s) + 11.97129741027622/exp(0.016285467000648854*s))
199  corr_ij(8,2)=0
200  corr_ij(8,1)=0
201  corr_ij(7,10)=0
202 
203  corr_ij(7,9)=0
204  corr_ij(7,8)=0
205  corr_ij(7,7)=((11.518026982819887*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.05351107793747755*c&
206  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.14054811601869432*sin(0.0293414097268719&
207  &26*s))/exp(0.04435489221745234*s) + (0.14054811601869702*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
208  corr_ij(7,6)=((0.14054811601869532*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) - (0.14054811601869702*&
209  &cos(0.11425932545092894*s))/exp(0.019700737327669783*s) + (11.518026982819887*sin(0.0293414097268719&
210  &26*s))/exp(0.04435489221745234*s) + (0.0535110779374777*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
211  corr_ij(7,5)=0
212  corr_ij(7,4)=0
213  corr_ij(7,3)=0
214  corr_ij(7,2)=((-0.732907009016115*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (2.728845031386875*cos&
215  &(0.11425932545092894*s))/exp(0.019700737327669783*s) - (2.4717920234033532*sin(0.029341409726871926*&
216  &s))/exp(0.04435489221745234*s) - (0.24003801347124257*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
217  corr_ij(7,1)=((2.4717920234033532*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.2400380134712426*co&
218  &s(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.7329070090161153*sin(0.029341409726871926&
219  &*s))/exp(0.04435489221745234*s) + (2.728845031386876*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
220  corr_ij(6,10)=0
221 
222  corr_ij(6,9)=0
223  corr_ij(6,8)=0
224  corr_ij(6,7)=((-0.1405481160186977*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.14054811601869713*&
225  &cos(0.11425932545092894*s))/exp(0.019700737327669783*s) - (11.518026982819885*sin(0.0293414097268719&
226  &26*s))/exp(0.04435489221745234*s) - (0.05351107793747755*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
227  corr_ij(6,6)=((11.518026982819885*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.05351107793747768*c&
228  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.14054811601869832*sin(0.0293414097268719&
229  &26*s))/exp(0.04435489221745234*s) + (0.14054811601869707*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
230  corr_ij(6,5)=0
231  corr_ij(6,4)=0
232  corr_ij(6,3)=0
233  corr_ij(6,2)=((-2.471792023403353*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) - (0.2400380134712425*co&
234  &s(0.11425932545092894*s))/exp(0.019700737327669783*s) + (0.7329070090161155*sin(0.029341409726871926&
235  &*s))/exp(0.04435489221745234*s) - (2.7288450313868755*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
236  corr_ij(6,1)=((-0.7329070090161154*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (2.728845031386876*co&
237  &s(0.11425932545092894*s))/exp(0.019700737327669783*s) - (2.4717920234033524*sin(0.029341409726871926&
238  &*s))/exp(0.04435489221745234*s) - (0.24003801347124343*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
239  corr_ij(5,10)=((0.5794534449999711*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (4.136986570427212*cos&
240  &(0.07283568782600224*s))/exp(0.017262015588746404*s) - (1.0360597341248128*sin(0.039597160512071454*&
241  &s))/exp(0.06567483898489704*s) + (3.167330918996692*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
242  corr_ij(5,9)=((1.0360597341248134*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) - (3.1673309189966856*co&
243  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) + (0.5794534449999746*sin(0.039597160512071454&
244  &*s))/exp(0.06567483898489704*s) + (4.1369865704272115*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
245  corr_ij(5,8)=0
246  corr_ij(5,7)=0
247  corr_ij(5,6)=0
248  corr_ij(5,5)=((-0.37825091063447547*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (30.094690926061638*&
249  &cos(0.07283568782600224*s))/exp(0.017262015588746404*s) + (0.16085380971100194*sin(0.039597160512071&
250  &454*s))/exp(0.06567483898489704*s) - (0.1608538097109995*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
251  corr_ij(5,4)=((-0.16085380971100238*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (0.16085380971100127&
252  &*cos(0.07283568782600224*s))/exp(0.017262015588746404*s) - (0.37825091063447586*sin(0.03959716051207&
253  &1454*s))/exp(0.06567483898489704*s) + (30.09469092606163*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
254  corr_ij(5,3)=0
255  corr_ij(5,2)=0
256  corr_ij(5,1)=0
257  corr_ij(4,10)=((-1.0360597341248106*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (3.167330918996689*co&
258  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) - (0.5794534449999716*sin(0.039597160512071454&
259  &*s))/exp(0.06567483898489704*s) - (4.1369865704272115*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
260  corr_ij(4,9)=((0.5794534449999711*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (4.1369865704272115*co&
261  &s(0.07283568782600224*s))/exp(0.017262015588746404*s) - (1.0360597341248114*sin(0.039597160512071454&
262  &*s))/exp(0.06567483898489704*s) + (3.1673309189966843*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
263  corr_ij(4,8)=0
264  corr_ij(4,7)=0
265  corr_ij(4,6)=0
266  corr_ij(4,5)=((0.16085380971100194*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) - (0.16085380971100371*&
267  &cos(0.07283568782600224*s))/exp(0.017262015588746404*s) + (0.37825091063447497*sin(0.039597160512071&
268  &454*s))/exp(0.06567483898489704*s) - (30.094690926061617*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
269  corr_ij(4,4)=((-0.37825091063447536*cos(0.039597160512071454*s))/exp(0.06567483898489704*s) + (30.094690926061617*&
270  &cos(0.07283568782600224*s))/exp(0.017262015588746404*s) + (0.16085380971100172*sin(0.039597160512071&
271  &454*s))/exp(0.06567483898489704*s) - (0.16085380971100616*sin(0.07283568782600224*s))/exp(0.017262015588746404*s))
272  corr_ij(4,3)=0
273  corr_ij(4,2)=0
274  corr_ij(4,1)=0
275  corr_ij(3,10)=0
276 
277  corr_ij(3,9)=0
278  corr_ij(3,8)=(0.24013456462471527/exp(0.05076718239027029*s) + 5.792596760796093/exp(0.016285467000648854*s))
279  corr_ij(3,7)=0
280  corr_ij(3,6)=0
281  corr_ij(3,5)=0
282  corr_ij(3,4)=0
283  corr_ij(3,3)=(-0.15609785880208227/exp(0.05076718239027029*s) + 31.18882918422289/exp(0.016285467000648854*s))
284  corr_ij(3,2)=0
285  corr_ij(3,1)=0
286  corr_ij(2,10)=0
287 
288  corr_ij(2,9)=0
289  corr_ij(2,8)=0
290  corr_ij(2,7)=((1.6172201305728584*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.37871789179790255*c&
291  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) + (1.2889451151208258*sin(0.02934140972687192&
292  &6*s))/exp(0.04435489221745234*s) + (1.4228849217537705*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
293  corr_ij(2,6)=((-1.2889451151208255*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) - (1.4228849217537702*c&
294  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) + (1.6172201305728586*sin(0.02934140972687192&
295  &6*s))/exp(0.04435489221745234*s) + (0.3787178917979035*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
296  corr_ij(2,5)=0
297  corr_ij(2,4)=0
298  corr_ij(2,3)=0
299  corr_ij(2,2)=((0.1789135645266575*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (26.817024457844113*co&
300  &s(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.4268927977731004*sin(0.029341409726871926&
301  &*s))/exp(0.04435489221745234*s) + (0.4268927977730982*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
302  corr_ij(2,1)=((0.4268927977731007*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) - (0.42689279777309963*c&
303  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) + (0.17891356452665746*sin(0.0293414097268719&
304  &26*s))/exp(0.04435489221745234*s) + (26.81702445784412*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
305  corr_ij(1,10)=0
306 
307  corr_ij(1,9)=0
308  corr_ij(1,8)=0
309  corr_ij(1,7)=((1.288945115120824*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (1.4228849217537711*cos&
310  &(0.11425932545092894*s))/exp(0.019700737327669783*s) - (1.617220130572856*sin(0.029341409726871926*s&
311  &))/exp(0.04435489221745234*s) - (0.3787178917979028*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
312  corr_ij(1,6)=((1.6172201305728564*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.37871789179790377*c&
313  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) + (1.2889451151208242*sin(0.02934140972687192&
314  &6*s))/exp(0.04435489221745234*s) + (1.4228849217537711*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
315  corr_ij(1,5)=0
316  corr_ij(1,4)=0
317  corr_ij(1,3)=0
318  corr_ij(1,2)=((-0.4268927977731002*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (0.4268927977730981*c&
319  &os(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.1789135645266573*sin(0.02934140972687192&
320  &6*s))/exp(0.04435489221745234*s) - (26.81702445784412*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
321  corr_ij(1,1)=((0.1789135645266574*cos(0.029341409726871926*s))/exp(0.04435489221745234*s) + (26.817024457844113*co&
322  &s(0.11425932545092894*s))/exp(0.019700737327669783*s) - (0.42689279777310024*sin(0.02934140972687192&
323  &6*s))/exp(0.04435489221745234*s) + (0.4268927977730997*sin(0.11425932545092894*s))/exp(0.019700737327669783*s))
324 
325  corr_ij=q_au**2*corr_ij
326 
327  END SUBROUTINE corrcomp_from_def
328 
329  !> Subroutine to compute the correlation of the unresolved variables \f$\langle Y \otimes Y^s \rangle\f$
330  !> at time \f$s\f$ from the spline representation
331  !> @param s time \f$s\f$ at which the correlation is computed
332  SUBROUTINE corrcomp_from_spline(s)
333  REAL(KIND=8), INTENT(IN) :: s
334  REAL(KIND=8) :: y
335  INTEGER :: i,j
336  corr_ij=0.d0
337  DO i=1,n_unres
338  DO j=1,n_unres
339  CALL splint(xa,ya(i,j,:),y2(i,j,:),nspl,s,y)
340  corr_ij(i,j)=y
341  END DO
342  END DO
343  END SUBROUTINE corrcomp_from_spline
344 
345  !> Routine to compute the spline representation parameters
346  SUBROUTINE splint(xa,ya,y2a,n,x,y)
347  INTEGER, INTENT(IN) :: n
348  REAL(KIND=8), INTENT(IN), DIMENSION(n) :: xa,y2a,ya
349  REAL(KIND=8), INTENT(IN) :: x
350  REAL(KIND=8), INTENT(OUT) :: y
351  INTEGER :: k
352  REAL(KIND=8) :: a,b,h
353  if ((khi-klo.gt.1).or.(xa(klo).gt.x).or.(xa(khi).lt.x)) then
354  if ((khi-klo.eq.1).and.(xa(klo).lt.x)) then
355  khi=klo
356  DO WHILE (xa(khi).lt.x)
357  khi=khi+1
358  END DO
359  klo=khi-1
360  else
361  khi=n
362  klo=1
363  DO WHILE (khi-klo.gt.1)
364  k=(khi+klo)/2
365  if(xa(k).gt.x)then
366  khi=k
367  else
368  klo=k
369  endif
370  END DO
371  end if
372  ! print*, "search",x,khi-klo,xa(klo),xa(khi)
373  ! else
374  ! print*, "ok",x,khi-klo,xa(klo),xa(khi)
375  endif
376  h=xa(khi)-xa(klo)
377  if (h.eq.0.) stop 'bad xa input in splint'
378  a=(xa(khi)-x)/h
379  b=(x-xa(klo))/h
380  y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
381  return
382  END SUBROUTINE splint
383 
384  !> Exponential fit function
385  !> @param s time \f$s\f$ at which the function is evaluated
386  !> @param p vector holding the coefficients of the fit function
387  FUNCTION fs(s,p)
388  REAL(KIND=8), INTENT(IN) :: s
389  REAL(KIND=8), DIMENSION(5), INTENT(IN) :: p
390  REAL(KIND=8) :: fs
391  fs=p(1)*exp(-s/p(2))*cos(p(3)*s+p(4))
392  RETURN
393  END FUNCTION fs
394 
395  !> Subroutine to compute the correlation of the unresolved variables \f$\langle Y \otimes Y^s \rangle\f$
396  !> at time \f$s\f$ from the exponential representation
397  !> @param s time \f$s\f$ at which the correlation is computed
398  SUBROUTINE corrcomp_from_fit(s)
399  REAL(KIND=8), INTENT(IN) :: s
400  REAL(KIND=8) :: y
401  INTEGER :: i,j
402 
403  corr_ij=0.d0
404  DO i=1,n_unres
405  DO j=1,n_unres
406  corr_ij(i,j)=fs(s,ya(i,j,:))
407  END DO
408  END DO
409  END SUBROUTINE corrcomp_from_fit
410 
411 
412 END MODULE corrmod
The stochastic models parameters module.
real(kind=8) q_au
Atmospheric unresolved component noise amplitude.
integer khi
Definition: corrmod.f90:38
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
real(kind=8), dimension(:,:), allocatable, public corr_ij
Matrix holding the correlation matrix at a given time.
Definition: corrmod.f90:32
character(len=4) load_mode
Loading mode for the correlations.
real(kind=8), dimension(:,:), allocatable, public corr_i_full
Covariance matrix of the unresolved variables (full version)
Definition: corrmod.f90:28
Utility module.
Definition: util.f90:12
Statistics accumulators.
Definition: stat.f90:14
real(kind=8), dimension(:,:,:), allocatable y2
Vector holding coefficient of the spline and exponential correlation representation.
Definition: corrmod.f90:33
real(kind=8) q_or
Oceanic resolved component noise amplitude.
real(kind=8), dimension(:), allocatable, public mean_full
Vector holding the mean of the unresolved dynamics (full version)
Definition: corrmod.f90:27
subroutine corrcomp_from_spline(s)
Subroutine to compute the correlation of the unresolved variables at time from the spline represent...
Definition: corrmod.f90:333
subroutine corrcomp_from_fit(s)
Subroutine to compute the correlation of the unresolved variables at time from the exponential repr...
Definition: corrmod.f90:399
subroutine, public init_corr
Subroutine to initialise the computation of the correlation.
Definition: corrmod.f90:46
real(kind=8), dimension(:), allocatable, public mean
Vector holding the mean of the unresolved dynamics (reduced version)
Definition: corrmod.f90:26
character(len=4) int_corr_mode
Correlation integration mode.
real(kind=8), dimension(:), allocatable xa
Vector holding coefficient of the spline and exponential correlation representation.
Definition: corrmod.f90:35
real(kind=8) q_ou
Oceanic unresolved component noise amplitude.
real(kind=8), dimension(:,:), allocatable, public inv_corr_i_full
Inverse of the covariance matrix of the unresolved variables (full version)
Definition: corrmod.f90:29
real(kind=8), dimension(:,:,:), allocatable ya
Vector holding coefficient of the spline and exponential correlation representation.
Definition: corrmod.f90:34
subroutine corrcomp_from_def(s)
Subroutine to compute the correlation of the unresolved variables at time from the definition given...
Definition: corrmod.f90:148
integer, public n_unres
Number of unresolved variables.
Definition: sf_def.f90:26
real(kind=8) function, dimension(size(a, 1), size(a, 2)), public invmat(A)
Definition: util.f90:218
Module to initialize the correlation matrix of the unresolved variables.
Definition: corrmod.f90:16
real(kind=8), dimension(:,:), allocatable, public corr_i
Covariance matrix of the unresolved variables (reduced version)
Definition: corrmod.f90:30
integer, dimension(:), allocatable, public ind
Definition: sf_def.f90:24
procedure(corrcomp_from_spline), pointer, public corrcomp
Pointer to the correlation computation routine.
Definition: corrmod.f90:41
The model parameters module.
Definition: params.f90:18
subroutine splint(xa, ya, y2a, n, x, y)
Routine to compute the spline representation parameters.
Definition: corrmod.f90:347
real(kind=8) maxint
Upper integration limit of the correlations.
integer nspl
Integers needed by the spline representation of the correlation.
Definition: corrmod.f90:38
Module to select the resolved-unresolved components.
Definition: sf_def.f90:12
real(kind=8), dimension(:,:), allocatable, public inv_corr_i
Inverse of the covariance matrix of the unresolved variables (reduced version)
Definition: corrmod.f90:31
real(kind=8) q_ar
Atmospheric resolved component noise amplitude.
real(kind=8) function fs(s, p)
Exponential fit function.
Definition: corrmod.f90:388
integer klo
Definition: corrmod.f90:38