A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
MAR.f90
Go to the documentation of this file.
1 
2 ! MAR.f90
3 !
4 !> Multidimensional Autoregressive module to generate the correlation for the
5 !> WL parameterization.
6 !
7 !> @copyright
8 !> 2018 Jonathan Demaeyer.
9 !> See LICENSE.txt for license information.
10 !
11 !---------------------------------------------------------------------------!
12 !
13 !> @remark
14 !> Based on the equation
15 !> \f$y_n = \sum_{i=1}^m y_{n-i} \cdot W_i + Q \cdot \xi_n
16 !> for an order \f$m\f$ MAR with \xi_n a Gaussian random variable.
17 !---------------------------------------------------------------------------!
18 
19 MODULE mar
20  USE params, only: ndim
21  USE sf_def, only: n_unres,ind,rind
22  USE sqrt_mod, only: init_sqrt,sqrtm
23  USE util, only: ireduce
24  USE stoch_mod, only: gasdev
25  IMPLICIT NONE
26 
27  PRIVATE
28 
29  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: q !< Square root of the noise covariance matrix
30  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: qred !< Reduce version of Q
31  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: rred !< Covariance matrix of the noise
32  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: w !< W_i matrix
33  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: wred !< Reduce W_i matrix
34  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y,dw
35 
36  INTEGER :: ms !< order of the MAR
37 
39  PUBLIC :: wred,qred,rred
40 
41 CONTAINS
42 
43  !> Subroutine to initialise the MAR
44  SUBROUTINE init_mar
45  INTEGER :: AllocStat,nf,i,info,info2
46  INTEGER, DIMENSION(3) :: s
47 
48  print*, 'Initializing the MAR integrator...'
49 
50  print*, 'Loading the MAR config from files...'
51 
52  OPEN(20,file='MAR_R_params.def',status='old')
53  READ(20,*) nf,ms
54  IF (nf /= n_unres) stop "*** Dimension in files MAR_R_params.def and sf.nml do not correspond ! ***"
55  ALLOCATE(qred(n_unres,n_unres),rred(n_unres,n_unres),wred(ms,n_unres,n_unres), stat=allocstat)
56  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
57  ALLOCATE(q(ndim,ndim),w(ms,ndim,ndim), stat=allocstat)
58  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
59  ALLOCATE(buf_y(0:ndim), dw(ndim), stat=allocstat)
60  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
61  READ(20,*) rred
62  CLOSE(20)
63 
64  OPEN(20,file='MAR_W_params.def',status='old')
65  READ(20,*) nf,ms
66  s=shape(wred)
67  IF (nf /= n_unres) stop "*** Dimension in files MAR_W_params.def and sf.nml do not correspond ! ***"
68  IF (s(1) /= ms) stop "*** MAR order in files MAR_R_params.def and MAR_W_params.def do not correspond ! ***"
69  DO i=1,ms
70  READ(20,*) wred(i,:,:)
71  ENDDO
72  CLOSE(20)
73 
74  CALL init_sqrt
75  CALL sqrtm(rred,qred,info,info2)
77 
78  DO i=1,ms
79  CALL ireduce(w(i,:,:),wred(i,:,:),n_unres,ind,rind)
80  ENDDO
81 
82  ! Kept for internal testing - Uncomment if not needed
83  ! DEALLOCATE(Wred,Rred,Qred, STAT=AllocStat)
84  ! IF (AllocStat /= 0) STOP "*** Deallocation problem ! ***"
85 
86  print*, 'MAR of order',ms,'found!'
87 
88  END SUBROUTINE init_mar
89 
90  !> Routine to generate one step of the MAR
91  !> @param x State vector of the MAR (store the \f$y_i\f$)
92  SUBROUTINE mar_step(x)
93  REAL(KIND=8), DIMENSION(0:ndim,ms), INTENT(INOUT) :: x
94  INTEGER :: j
95 
96  CALL stoch_vec(dw)
97  buf_y=0.d0
98  buf_y(1:ndim)=matmul(q,dw)
99  DO j=1,ms
100  buf_y(1:ndim)=buf_y(1:ndim)+matmul(x(1:ndim,j),w(j,:,:))
101  ENDDO
102  x=eoshift(x,shift=-1,boundary=buf_y,dim=2)
103  END SUBROUTINE mar_step
104 
105 
106  !> Routine to generate one step of the reduce MAR
107  !> @param xred State vector of the MAR (store the \f$y_i\f$)
108  !> @remark For debugging purpose only
109  SUBROUTINE mar_step_red(xred)
110  REAL(KIND=8), DIMENSION(0:ndim,ms), INTENT(INOUT) :: xred
111  INTEGER :: j
112 
113  CALL stoch_vec(dw)
114  buf_y=0.d0
115  buf_y(1:n_unres)=matmul(qred,dw(1:n_unres))
116  DO j=1,ms
117  buf_y(1:n_unres)=buf_y(1:n_unres)+matmul(xred(1:n_unres,j),wred(j,:,:))
118  ENDDO
119  xred=eoshift(xred,shift=-1,boundary=buf_y,dim=2)
120  END SUBROUTINE mar_step_red
121 
122 
123 
124  SUBROUTINE stoch_vec(dW)
125  REAL(KIND=8), DIMENSION(ndim), INTENT(INOUT) :: dW
126  INTEGER :: i
127  DO i=1,ndim
128  dw(i)=gasdev()
129  ENDDO
130  END SUBROUTINE stoch_vec
131 
132 
133 
134 END MODULE mar
subroutine, public mar_step(x)
Routine to generate one step of the MAR.
Definition: MAR.f90:93
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
real(kind=8), dimension(:,:), allocatable, public q
Square root of the noise covariance matrix.
Definition: MAR.f90:29
Utility module.
Definition: util.f90:12
real(kind=8), dimension(:,:,:), allocatable, public w
W_i matrix.
Definition: MAR.f90:32
real(kind=8) function, public gasdev()
Definition: stoch_mod.f90:32
Statistics accumulators.
Definition: stat.f90:14
real(kind=8), dimension(:), allocatable dw
Definition: MAR.f90:34
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, dimension(:), allocatable, public rind
Unresolved reduction indices.
Definition: sf_def.f90:24
subroutine, public ireduce(A, Ared, n, ind, rind)
Definition: util.f90:314
subroutine, public mar_step_red(xred)
Routine to generate one step of the reduce MAR.
Definition: MAR.f90:110
real(kind=8), dimension(:,:,:), allocatable, public wred
Reduce W_i matrix.
Definition: MAR.f90:33
subroutine stoch_vec(dW)
Definition: MAR.f90:125
subroutine, public init_sqrt
Definition: sqrt_mod.f90:39
integer, public n_unres
Number of unresolved variables.
Definition: sf_def.f90:26
real(kind=8), dimension(:,:), allocatable, public qred
Reduce version of Q.
Definition: MAR.f90:30
integer, dimension(:), allocatable, public ind
Definition: sf_def.f90:24
The model parameters module.
Definition: params.f90:18
Utility module containing the stochastic related routines.
Definition: stoch_mod.f90:15
Multidimensional Autoregressive module to generate the correlation for the WL parameterization.
Definition: MAR.f90:19
real(kind=8), dimension(:,:), allocatable, public rred
Covariance matrix of the noise.
Definition: MAR.f90:31
Module to select the resolved-unresolved components.
Definition: sf_def.f90:12
Utility module with various routine to compute matrix square root.
Definition: sqrt_mod.f90:20
real(kind=8), dimension(:), allocatable buf_y
Definition: MAR.f90:34
integer, public ms
order of the MAR
Definition: MAR.f90:36
subroutine, public init_mar
Subroutine to initialise the MAR.
Definition: MAR.f90:45