A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
MTV_int_tensor.f90
Go to the documentation of this file.
1 
2 ! MTV_int_tensor.f90
3 !
4 !> The MTV tensors used to integrate the MTV model
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !> See : Franzke, C., Majda, A. J., & Vanden-Eijnden, E. (2005). Low-order
14 !> stochastic mode reduction for a realistic barotropic model climate.
15 !> Journal of the atmospheric sciences, 62(6), 1722-1745.
16 !
17 !---------------------------------------------------------------------------!
18 
20 
21  !-----------------------------------------------------!
22  ! !
23  ! Preamble and variables declaration !
24  ! !
25  !-----------------------------------------------------!
26 
27  USE tensor
28  USE dec_tensor
29  USE corrmod
30  USE int_corr
31  USE params, only:ndim
32  USE stoch_params, only:mode
33  USE util, only: mat_trace, mat_contract
34 
35  IMPLICIT NONE
36 
37  PRIVATE
38 
39  PUBLIC :: init_mtv_int_tensor
40 
41  ! Constant vectors
42  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: h1 !< First constant vector
43  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: h2 !< Second constant vector
44  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: h3 !< Third constant vector
45  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: htot !< Total constant vector
46 
47  ! Tensors for the linear terms
48  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: l1 !< First linear tensor
49  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: l2 !< Second linear tensor
50  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: l3 !< Third linear tensor
51  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: ltot !< Total linear tensor
52 
53  ! Tensors for the quadratic terms
54  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: b1 !< First quadratic tensor
55  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: b2 !< Second quadratic tensor
56  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: btot !< Total quadratic tensor
57 
58  TYPE(coolist4), DIMENSION(:), ALLOCATABLE, PUBLIC :: mtot !< Tensor for the cubic terms
59 
60  ! Tensors for the noise covariance matrices
61  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: q1 !< Constant terms for the state-dependent noise covariance matrix
62  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: q2 !< Constant terms for the state-independent noise covariance matrix
63  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: utot !< Linear terms for the state-dependent noise covariance matrix
64  TYPE(coolist4), DIMENSION(:), ALLOCATABLE, PUBLIC :: vtot !< Quadratic terms for the state-dependent noise covariance matrix
65 
66  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dumb_vec !< Dummy vector
67  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat1 !< Dummy matrix
68  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat2 !< Dummy matrix
69  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat3 !< Dummy matrix
70  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat4 !< Dummy matrix
71 
72 
73  !-----------------------------------------------------!
74  ! !
75  ! End of preamble !
76  ! !
77  !-----------------------------------------------------!
78 
79 CONTAINS
80 
81  !-----------------------------------------------------!
82  ! !
83  ! Initialisation routine !
84  ! !
85  !-----------------------------------------------------!
86 
87  !> Subroutine to initialise the MTV tensor
88  SUBROUTINE init_mtv_int_tensor
89  INTEGER :: AllocStat,i,j,k,l
90 
91  print*, 'Initializing the decomposition tensors...'
92  CALL init_dec_tensor
93  print*, "Initializing the correlation matrices and tensors..."
94  CALL init_corrint
95  print*, "Computing the correlation integrated matrices and tensors..."
96  CALL comp_corrint
97 
98  !H part
99  print*, "Computing the H term..."
100 
101  ALLOCATE(h1(0:ndim), h2(0:ndim), h3(0:ndim), htot(0:ndim), stat=allocstat)
102  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
103 
104  ALLOCATE(dumb_mat1(ndim,ndim), dumb_mat2(ndim,ndim), stat=allocstat)
105  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
106 
107  ALLOCATE(dumb_mat3(ndim,ndim), dumb_mat4(ndim,ndim), stat=allocstat)
108  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
109 
110  !H1
112  dumb_mat2=matmul(dumb_mat1,corrint)
114 
115  ! H2
116  h2=0.d0
117  IF (mode.ne.'ures') THEN
120 
121  DO i=1,ndim
122  CALL coo_to_mat_i(i,bxyy,dumb_mat2)
127  ENDDO
128  ENDIF
129 
130  !H3
131  h3=0.d0
133 
134  !Htot
135  htot=0.d0
136  htot=h1+h2+h3
137 
138  print*, "Computing the L terms..."
139  ALLOCATE(l1(ndim), l2(ndim), l3(ndim), stat=allocstat)
140  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
141 
142  !L1
146  dumb_mat4=matmul(dumb_mat2,matmul(transpose(dumb_mat3),dumb_mat1))
147  CALL matc_to_coo(dumb_mat4,l1)
148 
149  !L2
150  dumb_mat4=0.d0
151  DO i=1,ndim
152  DO j=1,ndim
153  CALL coo_to_mat_i(i,bxyy,dumb_mat1)
155 
156  CALL coo_to_mat_j(j,byxy,dumb_mat1)
159  END DO
160  END DO
161  CALL matc_to_coo(dumb_mat4,l2)
162 
163  !L3
164  dumb_mat4=0.d0
165  DO i=1,ndim
166  DO j=1,ndim
167  CALL coo_to_mat_j(j,bxxy,dumb_mat1)
168  CALL coo_to_mat_i(i,bxxy,dumb_mat2)
169  dumb_mat4(i,j)=mat_trace(matmul(dumb_mat1,matmul(corrint,transpose(dumb_mat2))))
170  ENDDO
171  END DO
172  CALL matc_to_coo(dumb_mat4,l3)
173 
174  !Ltot
175 
176  ALLOCATE(ltot(ndim), stat=allocstat)
177  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
178 
179  CALL add_to_tensor(l1,ltot)
180  CALL add_to_tensor(l2,ltot)
181  CALL add_to_tensor(l3,ltot)
182 
183  print*, "Computing the B terms..."
184  ALLOCATE(b1(ndim), b2(ndim), btot(ndim), stat=allocstat)
185  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
186  ALLOCATE(dumb_vec(ndim), stat=allocstat)
187  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
188 
189  ! B1
192 
193  dumb_mat3=matmul(dumb_mat1,transpose(dumb_mat2))
194  DO j=1,ndim
195  DO k=1,ndim
196  CALL coo_to_vec_jk(j,k,byxx,dumb_vec)
197  dumb_vec=matmul(dumb_mat3,dumb_vec)
199  ENDDO
200  END DO
201 
202  ! B2
205 
206  dumb_mat4=matmul(transpose(dumb_mat2),dumb_mat3)
207  DO i=1,ndim
208  CALL coo_to_mat_i(i,bxxy,dumb_mat1)
211  ENDDO
212 
213  CALL add_to_tensor(b1,btot)
214  CALL add_to_tensor(b2,btot)
215 
216  !M
217 
218  print*, "Computing the M term..."
219 
220  ALLOCATE(mtot(ndim), stat=allocstat)
221  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
222 
224 
225  DO i=1,ndim
226  CALL coo_to_mat_i(i,bxxy,dumb_mat1)
227  dumb_mat3=matmul(dumb_mat1,transpose(dumb_mat2))
228  DO k=1,ndim
229  DO l=1,ndim
230  CALL coo_to_vec_jk(k,l,byxx,dumb_vec)
231  dumb_vec=matmul(dumb_mat3,dumb_vec)
233  ENDDO
234  END DO
235  END DO
236 
237  !Q
238 
239  print*, "Computing the Q terms..."
240  ALLOCATE(q1(ndim,ndim), q2(ndim,ndim), stat=allocstat)
241  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
242 
243  !Q1
244 
246  q1=matmul(dumb_mat1,matmul(corrint,transpose(dumb_mat1)))
247 
248  !Q2
249 
250  DO i=1,ndim
251  DO j=1,ndim
252  CALL coo_to_mat_i(i,bxyy,dumb_mat1)
253  CALL coo_to_mat_i(j,bxyy,dumb_mat2)
258  END DO
259  END DO
260 
261  !U
262 
263  ALLOCATE(utot(ndim), stat=allocstat)
264  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
265 
267  DO i=1,ndim
268  CALL coo_to_mat_i(i,bxxy,dumb_mat2)
269  dumb_mat3=matmul(dumb_mat1,matmul(corrint,transpose(dumb_mat2)))
271  ENDDO
272 
273  DO j=1,ndim
274  CALL coo_to_mat_i(j,bxxy,dumb_mat2)
275  dumb_mat3=matmul(dumb_mat1,matmul(corrint,transpose(dumb_mat2)))
276  DO k=1,ndim
277  CALL add_vec_jk_to_tensor(j,k,dumb_mat3(:,k),utot)
278  ENDDO
279  ENDDO
280 
281  !V
282 
283  ALLOCATE(vtot(ndim), stat=allocstat)
284  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
285 
286  DO i=1,ndim
287  DO j=1,ndim
288  CALL coo_to_mat_i(i,bxxy,dumb_mat1)
289  CALL coo_to_mat_i(j,bxxy,dumb_mat2)
290  dumb_mat3=matmul(dumb_mat1,matmul(corrint,transpose(dumb_mat2)))
292  ENDDO
293  ENDDO
294 
295  DEALLOCATE(dumb_mat1, dumb_mat2, stat=allocstat)
296  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
297 
298  DEALLOCATE(dumb_mat3, dumb_mat4, stat=allocstat)
299  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
300 
301  DEALLOCATE(dumb_vec, stat=allocstat)
302  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
303 
304 
305  END SUBROUTINE init_mtv_int_tensor
306 
307 
308 END MODULE mtv_int_tensor
309 
310 
311 
type(coolist), dimension(:), allocatable, public lxy
Tensor holding the linear part of the resolved tendencies involving the unresolved variables...
Definition: dec_tensor.f90:38
The stochastic models parameters module.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public add_matc_to_tensor4(i, j, src, dst)
Routine to add a matrix to a rank-4 tensor.
Definition: tensor.f90:537
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
subroutine, public comp_corrint
Routine that actually compute or load the integrals.
Definition: int_corr.f90:75
real(kind=8), dimension(:,:), allocatable, public corrint
Matrix holding the integral of the correlation matrix.
Definition: int_corr.f90:29
The MTV tensors used to integrate the MTV model.
type(coolist), dimension(:), allocatable, public l2
Second linear tensor.
Statistics accumulators.
Definition: stat.f90:14
subroutine, public add_vec_jk_to_tensor(j, k, src, dst)
Routine to add a vector to a rank-3 tensor.
Definition: tensor.f90:602
The resolved-unresolved components decomposition of the tensor.
Definition: dec_tensor.f90:16
real(kind=8), dimension(:), allocatable dumb_vec
Dummy vector.
subroutine, public sparse_mul4_with_mat_jl(coolist_ijkl, mat_jl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
Definition: tensor.f90:1169
type(coolist4), dimension(:), allocatable, public mtot
Tensor for the cubic terms.
type(coolist), dimension(:), allocatable, public l1
First linear tensor.
subroutine, public coo_to_mat_i(i, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1112
type(coolist), dimension(:), allocatable, public lyy
Tensor holding the linear part of the unresolved tendencies involving the unresolved variables...
Definition: dec_tensor.f90:45
subroutine, public coo_to_mat_j(j, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1148
Tensor utility module.
Definition: tensor.f90:18
real(kind=8) function, public mat_trace(A)
Definition: util.f90:150
real(kind=8) function, public mat_contract(A, B)
Definition: util.f90:162
type(coolist), dimension(:), allocatable, public btot
Total quadratic tensor.
type(coolist), dimension(:), allocatable, public b2
Second quadratic tensor.
real(kind=8), dimension(:), allocatable, public htot
Total constant vector.
real(kind=8), dimension(:), allocatable, public h1
First constant vector.
subroutine, public add_matc_to_tensor(i, src, dst)
Routine to add a matrix to a rank-3 tensor.
Definition: tensor.f90:474
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
type(coolist), dimension(:), allocatable, public ltot
Total linear tensor.
real(kind=8), dimension(:,:), allocatable dumb_mat3
Dummy matrix.
type(coolist4), dimension(:), allocatable, public corr2int
Tensor holding the integral of the correlation outer product with itself.
Definition: int_corr.f90:30
type(coolist), dimension(:), allocatable, public l3
Third linear tensor.
type(coolist), dimension(:), allocatable, public byxy
Tensor holding the quadratic part of the unresolved tendencies involving both resolved and unresolved...
Definition: dec_tensor.f90:47
real(kind=8), dimension(:,:), allocatable, public q2
Constant terms for the state-independent noise covariance matrix.
subroutine, public init_corrint
Subroutine to initialise the integrated matrices and tensors.
Definition: int_corr.f90:38
real(kind=8), dimension(:), allocatable, public h3
Third constant vector.
subroutine, public init_mtv_int_tensor
Subroutine to initialise the MTV tensor.
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
subroutine, public coo_to_mat_ik(src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix with i and k indices.
Definition: tensor.f90:1063
real(kind=8), dimension(:,:), allocatable, public q1
Constant terms for the state-dependent noise covariance matrix.
subroutine, public add_to_tensor(src, dst)
Routine to add a rank-3 tensor to another one.
Definition: tensor.f90:350
real(kind=8), dimension(:,:), allocatable dumb_mat4
Dummy matrix.
real(kind=8), dimension(:,:), allocatable dumb_mat1
Dummy matrix.
type(coolist), dimension(:), allocatable, public bxxy
Tensor holding the quadratic part of the resolved tendencies involving both resolved and unresolved v...
Definition: dec_tensor.f90:40
type(coolist), dimension(:), allocatable, public utot
Linear terms for the state-dependent noise covariance matrix.
4d coordinate list. Type used to represent the rank-4 sparse tensor.
Definition: tensor.f90:44
subroutine, public matc_to_coo(src, dst)
Routine to convert a matrix to a rank-3 tensor.
Definition: tensor.f90:1244
subroutine, public sparse_mul3_with_mat(coolist_ijk, mat_jk, res)
Sparse multiplication of a rank-3 tensor coolist with a matrix: .
Definition: tensor.f90:1220
subroutine, public add_vec_ikl_to_tensor4(i, k, l, src, dst)
Routine to add a vector to a rank-4 tensor.
Definition: tensor.f90:726
Module to initialize the correlation matrix of the unresolved variables.
Definition: corrmod.f90:16
type(coolist4), dimension(:), allocatable, public vtot
Quadratic terms for the state-dependent noise covariance matrix.
The model parameters module.
Definition: params.f90:18
type(coolist), dimension(:), allocatable, public bxyy
Tensor holding the quadratic part of the resolved tendencies involving unresolved variables...
Definition: dec_tensor.f90:41
real(kind=8), dimension(:), allocatable, public h2
Second constant vector.
type(coolist), dimension(:), allocatable, public lyx
Tensor holding the linear part of the unresolved tendencies involving the resolved variables...
Definition: dec_tensor.f90:44
type(coolist), dimension(:), allocatable, public byxx
Tensor holding the quadratic part of the unresolved tendencies involving resolved variables...
Definition: dec_tensor.f90:46
character(len=4) mode
Stochastic mode parameter.
real(kind=8), dimension(:,:), allocatable dumb_mat2
Dummy matrix.
subroutine, public init_dec_tensor
Subroutine that initialize and compute the decomposed tensors.
Definition: dec_tensor.f90:195
Module to compute or load the integrals of the correlation matrices.
Definition: int_corr.f90:16
subroutine, public coo_to_vec_jk(j, k, src, dst)
Routine to convert a rank-3 tensor coolist component into a vector.
Definition: tensor.f90:1129
type(coolist), dimension(:), allocatable, public b1
First quadratic tensor.