A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
rk2_MTV_integrator.f90
Go to the documentation of this file.
1 
2 ! rk2_MTV_integrator.f90
3 !
4 !> Module with the MTV rk2 integration routines.
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !> This module actually contains the Heun algorithm routines.
14 !
15 !---------------------------------------------------------------------------
16 
18  USE params, only: ndim
19  USE stoch_params, only: q_ar,q_au,q_or,q_ou,mnuti
21  USE mtv_int_tensor, only: ltot,htot,mtot,btot
22  USE aotensor_def, only:aotensor
23  USE sigma
24  USE stoch_mod
26  IMPLICIT NONE
27 
28  PRIVATE
29 
30  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1,buf_f0,buf_f1 !< Integration buffers
31  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dw, dwmult !< Standard gaussian noise buffers
32  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dwar,dwau,dwor,dwou !< Standard gaussian noise buffers
33  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: anoise,noise !< Additive noise term
34  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: noisemult !< Multiplicative noise term
35  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: g !< G term of the MTV tendencies
36  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_g !< Buffer for the G term computation
37 
38  LOGICAL :: mult !< Logical indicating if the sigma1 matrix must be computed for every state change
39  LOGICAL :: q1fill !< Logical indicating if the matrix Q1 is non-zero
40  LOGICAL :: compute_mult !< Logical indicating if the Gaussian noise for the multiplicative noise must be computed
41 
42  REAL(KIND=8), PARAMETER :: sq2 = sqrt(2.d0) !< Hard coded square root of 2
43 
44  PUBLIC :: init_integrator, step, full_step
45 
46 CONTAINS
47 
48  !> Subroutine to initialize the MTV rk2 integrator
49  SUBROUTINE init_integrator
50  INTEGER :: AllocStat
51 
52  CALL init_ss_integrator ! Initialize the uncoupled resolved dynamics
53 
54  ALLOCATE(buf_y1(0:ndim),buf_f0(0:ndim),buf_f1(0:ndim),stat=allocstat)
55  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
56 
57  buf_y1=0.d0
58  buf_f1=0.d0
59  buf_f0=0.d0
60 
61  print*, 'Initializing the integrator ...'
62  CALL init_sigma(mult,q1fill)
63  CALL init_noise
64  CALL init_g
65  END SUBROUTINE init_integrator
66 
67  !> Routine to initialize the noise vectors and buffers
68  SUBROUTINE init_noise
69  INTEGER :: AllocStat
70  ALLOCATE(dw(0:ndim), dwmult(0:ndim), stat=allocstat)
71  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
72 
73  ALLOCATE(dwar(0:ndim),dwau(0:ndim),dwor(0:ndim),dwou(0:ndim),stat=allocstat)
74  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
75 
76  ALLOCATE(anoise(0:ndim), noise(0:ndim), noisemult(0:ndim), stat=allocstat)
77  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
78 
79  dw=0.d0
80  dwmult=0.d0
81 
82  dwar=0.d0
83  dwor=0.d0
84  dwau=0.d0
85  dwou=0.d0
86 
87  anoise=0.d0
88  noise=0.d0
89  noisemult=0.d0
90 
91  compute_mult=((q1fill).OR.(mult))
92 
93  END SUBROUTINE init_noise
94 
95  !> Routine to initialize the G term
96  SUBROUTINE init_g
97  INTEGER :: AllocStat
98  ALLOCATE(g(0:ndim), buf_g(0:ndim), stat=allocstat)
99  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
100  END SUBROUTINE init_g
101 
102  !> Routine to actualize the G term based on the state y of the MTV system
103  !> @param y State of the MTV system
104  SUBROUTINE compg(y)
105  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
106 
107  g=htot
108  CALL sparse_mul2_k(ltot,y,buf_g)
109  g=g+buf_g
110  CALL sparse_mul3(btot,y,y,buf_g)
111  g=g+buf_g
112  CALL sparse_mul4(mtot,y,y,y,buf_g)
113  g=g+buf_g
114  END SUBROUTINE compg
115 
116  !> Routine to perform an integration step (Heun algorithm) of the MTV system. The incremented time is returned.
117  !> @param y Initial point.
118  !> @param t Actual integration time
119  !> @param dt Integration timestep.
120  !> @param dtn Stochastic integration timestep (normally square-root of dt).
121  !> @param res Final point after the step.
122  !> @param tend Partial or full tendencies used to perform the step (used for debugging).
123  SUBROUTINE step(y,t,dt,dtn,res,tend)
124  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
125  REAL(KIND=8), INTENT(INOUT) :: t
126  REAL(KIND=8), INTENT(IN) :: dt,dtn
127  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res,tend
128 
129  CALL compg(y)
130 
131  CALL stoch_atm_res_vec(dwar)
132  CALL stoch_oc_res_vec(dwor)
134  CALL stoch_vec(dw)
135  IF (compute_mult) CALL stoch_vec(dwmult)
136  noise(1:ndim)=matmul(sig2,dw(1:ndim))
137  IF ((mult).and.(mod(t,mnuti)<dt)) CALL compute_mult_sigma(y)
138  IF (compute_mult) noisemult(1:ndim)=matmul(sig1,dwmult(1:ndim))
139 
140  CALL tendencies(t,y,buf_f0)
141  buf_y1 = y+dt*(buf_f0+g)+(anoise+sq2*(noise+noisemult))*dtn
142 
143  buf_f1=g
144  CALL compg(buf_y1)
145  g=0.5*(g+buf_f1)
146 
147  IF ((mult).and.(mod(t,mnuti)<dt)) CALL compute_mult_sigma(buf_y1)
148  IF (compute_mult) THEN
149  buf_f1(1:ndim)=matmul(sig1,dwmult(1:ndim))
150  noisemult(1:ndim)=0.5*(noisemult(1:ndim)+buf_f1(1:ndim))
151  ENDIF
152 
153 
154  CALL tendencies(t,buf_y1,buf_f1)
155  buf_f0=0.5*(buf_f0+buf_f1)
156  res=y+dt*(buf_f0+g)+(anoise+sq2*(noise+noisemult))*dtn
157  ! tend=G+sq2*(noise+noisemult)/dtn
158  tend=sq2*noisemult/dtn
159  t=t+dt
160 
161  END SUBROUTINE step
162 
163  !> Routine to perform an integration step (Heun algorithm) of the full stochastic system. The incremented time is returned.
164  !> @param y Initial point.
165  !> @param t Actual integration time
166  !> @param dt Integration timestep.
167  !> @param dtn Stochastoc integration timestep (normally square-root of dt).
168  !> @param res Final point after the step.
169  SUBROUTINE full_step(y,t,dt,dtn,res)
170  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
171  REAL(KIND=8), INTENT(INOUT) :: t
172  REAL(KIND=8), INTENT(IN) :: dt,dtn
173  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
174  CALL stoch_atm_res_vec(dwar)
176  CALL stoch_oc_res_vec(dwor)
179  CALL sparse_mul3(aotensor,y,y,buf_f0)
180  buf_y1 = y+dt*buf_f0+anoise
182  res=y+0.5*(buf_f0+buf_f1)*dt+anoise
183  t=t+dt
184  END SUBROUTINE full_step
185 
186 END MODULE rk2_mtv_integrator
The stochastic models parameters module.
subroutine, public stoch_oc_unres_vec(dW)
routine to fill the unresolved oceanic component of a vector with standard gaussian noise process val...
Definition: stoch_mod.f90:120
real(kind=8) q_au
Atmospheric unresolved component noise amplitude.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public stoch_atm_unres_vec(dW)
routine to fill the unresolved atmospheric component of a vector with standard gaussian noise process...
Definition: stoch_mod.f90:88
The MTV tensors used to integrate the MTV model.
subroutine, public init_sigma(mult, Q1fill)
Subroutine to initialize the sigma matices.
real(kind=8), dimension(:), allocatable anoise
Statistics accumulators.
Definition: stat.f90:14
subroutine init_noise
Routine to initialize the noise vectors and buffers.
logical mult
Logical indicating if the sigma1 matrix must be computed for every state change.
subroutine, public tendencies(t, y, res)
Routine computing the tendencies of the uncoupled resolved model.
real(kind=8) q_or
Oceanic resolved component noise amplitude.
subroutine init_g
Routine to initialize the G term.
subroutine, public init_integrator
Subroutine to initialize the MTV rk2 integrator.
subroutine, public step(y, t, dt, dtn, res, tend)
Routine to perform an integration step (Heun algorithm) of the MTV system. The incremented time is re...
real(kind=8), dimension(:), allocatable dwmult
Standard gaussian noise buffers.
subroutine, public stoch_atm_res_vec(dW)
routine to fill the resolved atmospheric component of a vector with standard gaussian noise process v...
Definition: stoch_mod.f90:77
type(coolist4), dimension(:), allocatable, public mtot
Tensor for the cubic terms.
subroutine, public full_step(y, t, dt, dtn, res)
Routine to perform an integration step (Heun algorithm) of the full stochastic system. The incremented time is returned.
The equation tensor for the coupled ocean-atmosphere model with temperature which allows for an exten...
subroutine, public sparse_mul3(coolist_ijk, arr_j, arr_k, res)
Sparse multiplication of a tensor with two vectors: .
Definition: tensor.f90:129
subroutine, public init_ss_integrator
Subroutine to initialize the uncoupled resolved rk2 integrator.
Tensor utility module.
Definition: tensor.f90:18
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.
Module with the stochastic uncoupled resolved nonlinear and tangent linear rk2 dynamics integration r...
type(coolist), dimension(:), allocatable, public btot
Total quadratic tensor.
real(kind=8), dimension(:), allocatable buf_y1
real(kind=8), dimension(:), allocatable, public htot
Total constant vector.
Module with the MTV rk2 integration routines.
subroutine, public stoch_oc_res_vec(dW)
routine to fill the resolved oceanic component of a vector with standard gaussian noise process value...
Definition: stoch_mod.f90:109
type(coolist), dimension(:), allocatable, public ltot
Total linear tensor.
real(kind=8), dimension(:,:), allocatable, public sig1
state-dependent noise matrix
real(kind=8), dimension(:), allocatable dwau
real(kind=8), dimension(:), allocatable dwou
Standard gaussian noise buffers.
real(kind=8) q_ou
Oceanic unresolved component noise amplitude.
logical compute_mult
Logical indicating if the Gaussian noise for the multiplicative noise must be computed.
real(kind=8), parameter sq2
Hard coded square root of 2.
subroutine, public stoch_vec(dW)
Routine to fill a vector with standard Gaussian noise process values.
Definition: stoch_mod.f90:57
logical q1fill
Logical indicating if the matrix Q1 is non-zero.
subroutine, public sparse_mul2_k(coolist_ijk, arr_k, res)
Sparse multiplication of a rank-3 sparse tensor coolist with a vector: .
Definition: tensor.f90:1045
The model parameters module.
Definition: params.f90:18
Utility module containing the stochastic related routines.
Definition: stoch_mod.f90:15
real(kind=8), dimension(:), allocatable buf_g
Buffer for the G term computation.
real(kind=8), dimension(:,:), allocatable, public sig2
state-independent noise matrix
subroutine compg(y)
Routine to actualize the G term based on the state y of the MTV system.
real(kind=8), dimension(:), allocatable buf_f0
real(kind=8), dimension(:), allocatable buf_f1
Integration buffers.
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 aotensor
- Tensor representation of the tendencies.
real(kind=8), dimension(:), allocatable dwar
real(kind=8), dimension(:), allocatable dwor
real(kind=8), dimension(:), allocatable noise
Additive noise term.
real(kind=8), dimension(:), allocatable dw
real(kind=8), dimension(:), allocatable g
G term of the MTV tendencies.
real(kind=8) mnuti
Multiplicative noise update time interval.
real(kind=8) q_ar
Atmospheric resolved component noise amplitude.
real(kind=8), dimension(:), allocatable noisemult
Multiplicative noise term.