A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
rk2_ss_integrator.f90
Go to the documentation of this file.
1 
2 ! rk2_ss_integrator.f90
3 !
4 !> Module with the stochastic uncoupled resolved nonlinear and tangent linear
5 !> rk2 dynamics integration routines.
6 !
7 !> @copyright
8 !> 2018 Jonathan Demaeyer.
9 !> See LICENSE.txt for license information.
10 !
11 !---------------------------------------------------------------------------!
12 !
13 !> @remark
14 !> This module actually contains the Heun algorithm routines.
15 !
16 !---------------------------------------------------------------------------
17 
19  USE params, only: ndim,natm
20  USE stoch_params, only: q_ar,q_au,q_or,q_ou
21  USE tensor, only: sparse_mul3
23  USE stoch_mod
24  IMPLICIT NONE
25 
26  PRIVATE
27 
28  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dwar,dwor !< Standard gaussian noise buffers
29 
30  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: anoise !< Additive noise term
31 
32  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1,buf_f0,buf_f1 !< Integration buffers
33 
35 
36 CONTAINS
37 
38  !> Subroutine to initialize the uncoupled resolved rk2 integrator
39  SUBROUTINE init_ss_integrator
40  INTEGER :: AllocStat
41 
42  ALLOCATE(buf_y1(0:ndim),buf_f0(0:ndim),buf_f1(0:ndim),stat=allocstat)
43  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
44 
45  ALLOCATE(anoise(0:ndim),stat=allocstat)
46  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
47 
48  ALLOCATE(dwar(0:ndim),dwor(0:ndim),stat=allocstat)
49  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
50 
51  dwar=0.d0
52  dwor=0.d0
53 
54  END SUBROUTINE init_ss_integrator
55 
56  !> Routine computing the tendencies of the uncoupled resolved model
57  !> @param t Time at which the tendencies have to be computed. Actually not needed for autonomous systems.
58  !> @param y Point at which the tendencies have to be computed.
59  !> @param res vector to store the result.
60  !> @remark Note that it is NOT safe to pass `y` as a result buffer,
61  !> as this operation does multiple passes.
62  SUBROUTINE tendencies(t,y,res)
63  REAL(KIND=8), INTENT(IN) :: t
64  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
65  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
66  CALL sparse_mul3(ss_tensor, y, y, res)
67  END SUBROUTINE tendencies
68 
69  !> Tendencies for the tangent linear model of the uncoupled resolved dynamics
70  !> in point ystar for perturbation deltay.
71  !> @param t time
72  !> @param y point of the tangent space at which the tendencies have to be computed.
73  !> @param ys point in trajectory to which the tangent space belongs.
74  !> @param res vector to store the result.
75  SUBROUTINE tl_tendencies(t,y,ys,res)
76  REAL(KIND=8), INTENT(IN) :: t
77  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y,ys
78  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
79  CALL sparse_mul3(ss_tl_tensor, y, ys, res)
80  END SUBROUTINE tl_tendencies
81 
82  !> Routine to perform a stochastic integration step of the unresolved
83  !> uncoupled dynamics (Heun algorithm).
84  !> The incremented time is returned.
85  !> @param y Initial point.
86  !> @param ys Dummy argument for compatibility.
87  !> @param t Actual integration time
88  !> @param dt Integration timestep.
89  !> @param dtn Stochastic integration timestep (normally square-root of dt).
90  !> @param res Final point after the step.
91  SUBROUTINE ss_step(y,ys,t,dt,dtn,res)
92  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y,ys
93  REAL(KIND=8), INTENT(INOUT) :: t
94  REAL(KIND=8), INTENT(IN) :: dt,dtn
95  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
96 
99  anoise=(q_ar*dwar+q_or*dwor)*dtn
100  CALL tendencies(t,y,buf_f0)
101  buf_y1 = y+dt*buf_f0+anoise
102  CALL tendencies(t,buf_y1,buf_f1)
103  res=y+0.5*(buf_f0+buf_f1)*dt+anoise
104  t=t+dt
105  END SUBROUTINE ss_step
106 
107  !> Routine to perform a stochastic integration step of the unresolved
108  !> uncoupled tangent linear dynamics (Heun algorithm).
109  !> The incremented time is returned.
110  !> @param y Initial point.
111  !> @param ys point in trajectory to which the tangent space belongs.
112  !> @param t Actual integration time
113  !> @param dt Integration timestep.
114  !> @param dtn Stochastic integration timestep (normally square-root of dt).
115  !> @param res Final point after the step.
116  SUBROUTINE ss_tl_step(y,ys,t,dt,dtn,res)
117  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y,ys
118  REAL(KIND=8), INTENT(INOUT) :: t
119  REAL(KIND=8), INTENT(IN) :: dt,dtn
120  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
121 
122  CALL stoch_atm_res_vec(dwar)
123  CALL stoch_oc_res_vec(dwor)
124  anoise=(q_ar*dwar+q_or*dwor)*dtn
125  CALL tl_tendencies(t,y,ys,buf_f0)
126  buf_y1 = y+dt*buf_f0+anoise
127  CALL tl_tendencies(t,buf_y1,ys,buf_f1)
128  res=y+0.5*(buf_f0+buf_f1)*dt+anoise
129  t=t+dt
130  END SUBROUTINE ss_tl_step
131 
132 END MODULE rk2_ss_integrator
The stochastic models parameters module.
real(kind=8) q_au
Atmospheric unresolved component noise amplitude.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
Statistics accumulators.
Definition: stat.f90:14
The resolved-unresolved components decomposition of the tensor.
Definition: dec_tensor.f90:16
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, 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
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 ...
real(kind=8), dimension(:), allocatable dwor
Standard gaussian noise buffers.
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
subroutine, public init_ss_integrator
Subroutine to initialize the uncoupled resolved rk2 integrator.
Tensor utility module.
Definition: tensor.f90:18
real(kind=8), dimension(:), allocatable buf_f1
Integration buffers.
real(kind=8), dimension(:), allocatable buf_f0
Module with the stochastic uncoupled resolved nonlinear and tangent linear rk2 dynamics integration r...
type(coolist), dimension(:), allocatable, public ss_tl_tensor
Tensor of the tangent linear model tendencies of the resolved component alone.
Definition: dec_tensor.f90:50
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
real(kind=8), dimension(:), allocatable dwar
subroutine, public tl_tendencies(t, y, ys, res)
Tendencies for the tangent linear model of the uncoupled resolved dynamics in point ystar for perturb...
real(kind=8) q_ou
Oceanic unresolved component noise amplitude.
integer natm
Number of atmospheric basis functions.
Definition: params.f90:83
real(kind=8), dimension(:), allocatable buf_y1
type(coolist), dimension(:), allocatable, public ss_tensor
Tensor holding the part of the resolved tensor involving only resolved variables. ...
Definition: dec_tensor.f90:33
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 anoise
Additive noise term.
real(kind=8) q_ar
Atmospheric resolved component noise amplitude.