A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
memory.f90
Go to the documentation of this file.
1 
2 ! memory.f90
3 !
4 !> Module that compute the memory term \f$M_3\f$ of the WL parameterization.
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !
14 !---------------------------------------------------------------------------!
15 
16 MODULE memory
17 
18  USE tensor
19  USE wl_tensor
20  USE params, only: ndim
21  USE stoch_params, only: muti,mems,x_int_mode
23 
24  IMPLICIT NONE
25 
26  PRIVATE
27 
28  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: x !< Array storing the previous state of the system
29  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: xs !< Array storing the resolved time evolution of the previous state of the system
30  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zs !< Dummy array to replace Xs in case where the evolution is not stored
31  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_m !< Dummy vector
32  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_m3 !< Dummy vector to store the \f$M_3\f$ integrand
33 
34  INTEGER :: t_index !< Integer storing the time index (current position in the arrays)
35 
36  PROCEDURE(ss_step), POINTER :: step !< Procedural pointer pointing on the resolved dynamics step routine
37 
39 
40 
41 CONTAINS
42 
43  !> Subroutine to initialise the memory
44  SUBROUTINE init_memory
45  INTEGER :: AllocStat
46 
48 
49  ALLOCATE(x(0:ndim,mems), stat=allocstat)
50  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
51 
52  x=0.d0
53 
54  IF (b23def) THEN
55  ALLOCATE(xs(0:ndim,mems), zs(0:ndim,mems), stat=allocstat)
56  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
57 
58  xs=0.d0
59  ENDIF
60 
61  ALLOCATE(buf_m3(0:ndim), buf_m(0:ndim), stat=allocstat)
62  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
63 
64  SELECT CASE (x_int_mode)
65  CASE('reso')
66  step => ss_step
67  CASE('tang')
68  step => ss_tl_step
69  CASE DEFAULT
70  stop '*** X_INT_MODE variable not properly defined in stoch_params.nml ***'
71  END SELECT
72 
73  END SUBROUTINE init_memory
74 
75  !> Compute the integrand of \f$M_3\f$ at each time in the past and integrate
76  !> to get the memory term
77  !> @param y current state
78  !> @param dt timestep
79  !> @param dtn stochastic timestep
80  !> @param savey set if the state is stored in X at the end
81  !> @param save_ev set if the result of the resolved time evolution is stored in Xs at the end
82  !> @param evolve set if the resolved time evolution is performed
83  !> @param inter set over which time interval the resolved time evolution must be computed
84  !> @param h_int result of the integration - give the memory term
85  SUBROUTINE compute_m3(y,dt,dtn,savey,save_ev,evolve,inter,h_int)
86  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
87  REAL(KIND=8), INTENT(IN) :: dt,dtn
88  LOGICAL, INTENT(IN) :: savey,save_ev,evolve
89  REAL(KIND=8), INTENT(IN) :: inter
90  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: h_int
91  REAL(KIND=8) :: t
92  INTEGER :: i,j
93 
94  x(:,t_index)=y
95  IF (b23def) THEN
96  xs(:,t_index)=y
97  zs(:,t_index)=y
98  DO i=1,mems-1
99  j=modulo(t_index+i-1,mems)+1
100  zs(:,j)=xs(:,j)
101  IF (evolve) THEN
102  IF (dt.lt.inter) THEN
103  t=0.d0
104  DO WHILE (t+dt<inter)
105  CALL step(zs(:,j),y,t,dt,dtn,zs(:,j))
106  ENDDO
107  CALL step(zs(:,j),y,t,inter-t,sqrt(inter-t),zs(:,j))
108  ELSE
109  CALL step(zs(:,j),y,t,inter,sqrt(inter),zs(:,j))
110  ENDIF
111  ENDIF
112  IF (save_ev) xs(:,j)=zs(:,j)
113  ENDDO
114  ENDIF
115 
116 
117  ! Computing the integral
118  h_int=0.d0
119 
120  DO i=1,mems
121  j=modulo(t_index+i-2,mems)+1
122  buf_m3=0.d0
123  IF (ldef) THEN
124  CALL sparse_mul3(ltot(:,i),y,x(:,j),buf_m)
126  ENDIF
127  IF (b14def) THEN
128  CALL sparse_mul3(b14(:,i),x(:,j),x(:,j),buf_m)
130  ENDIF
131  IF (b23def) THEN
132  CALL sparse_mul3(b23(:,i),x(:,j),zs(:,j),buf_m)
134  ENDIF
135  IF (mdef) THEN
136  CALL sparse_mul4(mtot(:,i),x(:,j),x(:,j),zs(:,j),buf_m)
138  ENDIF
139  IF ((i.eq.1).or.(i.eq.mems)) THEN
140  h_int=h_int+0.5*buf_m3
141  ELSE
142  h_int=h_int+buf_m3
143  ENDIF
144  ENDDO
145 
146  h_int=muti*h_int
147  IF (savey) THEN
148  t_index=t_index-1
149  IF (t_index.eq.0) t_index=mems
150  ENDIF
151  END SUBROUTINE compute_m3
152 
153  !> Routine to test the #compute_M3 routine
154  !> @param y current state
155  !> @param dt timestep
156  !> @param dtn stochastic timestep
157  !> @param h_int result of the integration - give the memory term
158  SUBROUTINE test_m3(y,dt,dtn,h_int)
159  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
160  REAL(KIND=8), INTENT(IN) :: dt,dtn
161  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: h_int
162  INTEGER :: i,j
163 
164  CALL compute_m3(y,dt,dtn,.true.,.true.,.true.,muti,h_int)
165  print*, t_index
166  print*, 'X'
167  DO i=1,mems
168  j=modulo(t_index+i-1,mems)+1
169  print*, i,j,x(1,j)
170  ENDDO
171 
172  IF (b23def) THEN
173  print*, 'Xs'
174  DO i=1,mems
175  j=modulo(t_index+i-1,mems)+1
176  print*, i,j,xs(1,j)
177  ENDDO
178  ENDIF
179  print*, 'h_int',h_int
180  END SUBROUTINE test_m3
181 
182 END MODULE memory
183 
The stochastic models parameters module.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
real(kind=8), dimension(:,:), allocatable x
Array storing the previous state of the system.
Definition: memory.f90:28
real(kind=8) muti
Memory update time interval.
logical, public b14def
Definition: WL_tensor.f90:75
real(kind=8), dimension(:), allocatable buf_m
Dummy vector.
Definition: memory.f90:31
logical, public b23def
Definition: WL_tensor.f90:75
Statistics accumulators.
Definition: stat.f90:14
logical, public ldef
Definition: WL_tensor.f90:75
procedure(ss_step), pointer step
Procedural pointer pointing on the resolved dynamics step routine.
Definition: memory.f90:36
subroutine, public ss_tl_step(y, ys, t, dt, dtn, res)
Routine to perform a stochastic integration step of the unresolved uncoupled tangent linear dynamics ...
subroutine, public ss_step(y, ys, t, dt, dtn, res)
Routine to perform a stochastic integration step of the unresolved uncoupled dynamics (Heun algorithm...
subroutine, public sparse_mul3(coolist_ijk, arr_j, arr_k, res)
Sparse multiplication of a tensor with two vectors: .
Definition: tensor.f90:129
character(len=4) x_int_mode
Integration mode for the resolved component.
Tensor utility module.
Definition: tensor.f90:18
type(coolist), dimension(:,:), allocatable, public b14
Joint 1st and 4th tensors.
Definition: WL_tensor.f90:64
Module with the stochastic uncoupled resolved nonlinear and tangent linear rk2 dynamics integration r...
real(kind=8), dimension(:,:), allocatable xs
Array storing the resolved time evolution of the previous state of the system.
Definition: memory.f90:29
real(kind=8), dimension(:), allocatable buf_m3
Dummy vector to store the integrand.
Definition: memory.f90:32
subroutine, public test_m3(y, dt, dtn, h_int)
Routine to test the #compute_M3 routine.
Definition: memory.f90:159
The WL tensors used to integrate the model.
Definition: WL_tensor.f90:18
type(coolist4), dimension(:,:), allocatable, public mtot
Tensor for the cubic terms.
Definition: WL_tensor.f90:67
real(kind=8), dimension(:,:), allocatable zs
Dummy array to replace Xs in case where the evolution is not stored.
Definition: memory.f90:30
integer t_index
Integer storing the time index (current position in the arrays)
Definition: memory.f90:34
subroutine, public init_memory
Subroutine to initialise the memory.
Definition: memory.f90:45
The model parameters module.
Definition: params.f90:18
subroutine, public compute_m3(y, dt, dtn, savey, save_ev, evolve, inter, h_int)
Compute the integrand of at each time in the past and integrate to get the memory term...
Definition: memory.f90:86
type(coolist), dimension(:,:), allocatable, public b23
Joint 2nd and 3rd tensors.
Definition: WL_tensor.f90:65
logical, public mdef
Boolean to (de)activate the computation of the terms.
Definition: WL_tensor.f90:75
Module that compute the memory term of the WL parameterization.
Definition: memory.f90:16
subroutine, public sparse_mul4(coolist_ijkl, arr_j, arr_k, arr_l, res)
Sparse multiplication of a rank-4 tensor coolist with three vectors: .
Definition: tensor.f90:974
type(coolist), dimension(:,:), allocatable, public ltot
Total linear tensor.
Definition: WL_tensor.f90:57
integer mems
Number of steps in the memory kernel integral.