A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
MTV_sigma_tensor.f90
Go to the documentation of this file.
1 
2 ! MTV_sigma_tensor.f90
3 !
4 !> The MTV noise sigma matrices 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 
19 MODULE sigma
20  USE params, only:ndim
21  USE mtv_int_tensor, only: q1,q2,utot,vtot
22  USE tensor
25  IMPLICIT NONE
26 
27  PRIVATE
28 
30 
31  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: sig1 !< \f$\sigma_1(X)\f$ state-dependent noise matrix
32  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: sig2 !< \f$\sigma_2\f$ state-independent noise matrix
33  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: sig1r !< Reduced \f$\sigma_1(X)\f$ state-dependent noise matrix
34 
35  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat1 !< Dummy matrix
36  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat2 !< Dummy matrix
37  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat3 !< Dummy matrix
38  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat4 !< Dummy matrix
39  INTEGER, DIMENSION(:), ALLOCATABLE :: ind1,rind1,ind2,rind2 !< Reduction indices
40 
41  INTEGER :: n1,n2
42 
43 CONTAINS
44 
45 
46  !> Subroutine to initialize the sigma matices
47  SUBROUTINE init_sigma(mult,Q1fill)
48  LOGICAL, INTENT(OUT) :: mult,Q1fill
49  INTEGER :: AllocStat,info1,info2
50 
51  CALL init_sqrt
52 
53  ALLOCATE(sig1(ndim,ndim), sig2(ndim,ndim), sig1r(ndim,ndim),stat=allocstat)
54  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
55 
56  ALLOCATE(ind1(ndim), rind1(ndim), ind2(ndim), rind2(ndim), stat=allocstat)
57  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
58 
59  ALLOCATE(dumb_mat1(ndim,ndim), dumb_mat2(ndim,ndim), stat=allocstat)
60  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
61 
62  ALLOCATE(dumb_mat3(ndim,ndim), dumb_mat4(ndim,ndim), stat=allocstat)
63  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
64 
65  print*, "Initializing the sigma matrices"
66 
68  IF (n2 /= 0) THEN
69  CALL sqrtm_svd(dumb_mat1(1:n2,1:n2),dumb_mat2(1:n2,1:n2),info1,info2,min(max(n2/2,2),64))
71  ELSE
72  sig2=0.d0
73  ENDIF
74 
75  mult=(.not.((tensor_empty(utot)).and.(tensor4_empty(vtot))))
76  q1fill=.true.
78  IF (n1 /= 0) THEN
79 
80  CALL sqrtm_svd(dumb_mat1(1:n1,1:n1),dumb_mat2(1:n1,1:n1),info1,info2,min(max(n1/2,2),64))
82  ELSE
83  q1fill=.false.
84  sig1=0.d0
85  ENDIF
86  sig1r=sig1
87 
88  END SUBROUTINE init_sigma
89 
90  !> Routine to actualize the matrix \f$\sigma_1\f$ based on the state y of the MTV system
91  !> @param y State of the MTV system
92  SUBROUTINE compute_mult_sigma(y)
93  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
94  INTEGER :: info,info2
99  IF (n1 /= 0) THEN
100  CALL sqrtm_svd(dumb_mat1(1:n1,1:n1),dumb_mat2(1:n1,1:n1),info,info2,min(max(n1/2,2),64))
101  ! dumb_mat2=0.D0
102  ! CALL chol(0.5*(dumb_mat1(1:n1,1:n1)+transpose(dumb_mat1(1:n1,1:n1))),dumb_mat2(1:n1,1:n1),info)
103  IF ((.not.any(isnan(dumb_mat2))).and.(info.eq.0).and.(.not.any(dumb_mat2>huge(0.d0)))) THEN
105  ELSE
106  sig1=sig1r
107  ENDIF
108  ELSE
109  sig1=sig1r
110  ENDIF
111  END SUBROUTINE compute_mult_sigma
112 
113 
114 END MODULE sigma
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public chol(A, sqA, info)
Routine to perform a Cholesky decomposition.
Definition: sqrt_mod.f90:386
integer n2
Utility module.
Definition: util.f90:12
The MTV tensors used to integrate the MTV model.
subroutine, public init_sigma(mult, Q1fill)
Subroutine to initialize the sigma matices.
Statistics accumulators.
Definition: stat.f90:14
subroutine, public sqrtm(A, sqA, info, info_triu, bs)
Routine to compute a real square-root of a matrix.
Definition: sqrt_mod.f90:63
integer n1
subroutine, public printmat(A)
Definition: util.f90:200
subroutine, public ireduce(A, Ared, n, ind, rind)
Definition: util.f90:314
Tensor utility module.
Definition: tensor.f90:18
integer, dimension(:), allocatable ind2
real(kind=8), dimension(:,:), allocatable dumb_mat3
Dummy matrix.
The MTV noise sigma matrices used to integrate the MTV model.
subroutine, public compute_mult_sigma(y)
Routine to actualize the matrix based on the state y of the MTV system.
logical function, public tensor_empty(t)
Test if a rank-3 tensor coolist is empty.
Definition: tensor.f90:1288
real(kind=8), dimension(:,:), allocatable, public sig1
state-dependent noise matrix
real(kind=8), dimension(:,:), allocatable, public q2
Constant terms for the state-independent noise covariance matrix.
subroutine, public reduce(A, Ared, n, ind, rind)
Definition: util.f90:286
integer, dimension(:), allocatable rind2
Reduction indices.
real(kind=8), dimension(:,:), allocatable dumb_mat1
Dummy matrix.
real(kind=8), dimension(:,:), allocatable dumb_mat2
Dummy matrix.
subroutine, public sparse_mul3_mat(coolist_ijk, arr_k, res)
Sparse multiplication of a rank-3 tensor coolist with a vector: . Its output is a matrix...
Definition: tensor.f90:948
real(kind=8), dimension(:,:), allocatable, public q1
Constant terms for the state-dependent noise covariance matrix.
subroutine, public init_sqrt
Definition: sqrt_mod.f90:39
type(coolist), dimension(:), allocatable, public utot
Linear terms for the state-dependent noise covariance matrix.
type(coolist4), dimension(:), allocatable, public vtot
Quadratic terms for the state-dependent noise covariance matrix.
The model parameters module.
Definition: params.f90:18
integer, dimension(:), allocatable rind1
real(kind=8), dimension(:,:), allocatable, public sig2
state-independent noise matrix
subroutine, public cprintmat(A)
Definition: util.f90:209
logical function, public tensor4_empty(t)
Test if a rank-4 tensor coolist is empty.
Definition: tensor.f90:1304
integer, dimension(:), allocatable ind1
Utility module with various routine to compute matrix square root.
Definition: sqrt_mod.f90:20
subroutine, public sparse_mul4_mat(coolist_ijkl, arr_k, arr_l, res)
Sparse multiplication of a tensor with two vectors: .
Definition: tensor.f90:998
real(kind=8), dimension(:,:), allocatable dumb_mat4
Dummy matrix.
real(kind=8), dimension(:,:), allocatable, public sig1r
Reduced state-dependent noise matrix.
subroutine, public sqrtm_svd(A, sqA, info, info_triu, bs)
Routine to compute a real square-root of a matrix via a SVD decomposition.
Definition: sqrt_mod.f90:401