A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
rk2_stoch_integrator.f90
Go to the documentation of this file.
1 
2 ! rk2_stoch_integrator.f90
3 !
4 !> Module with the stochastic 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 !> There are four modes for this integrator:
15 !> - full: use the full dynamics
16 !> - ures: use the intrinsic unresolved dynamics
17 !> - qfst: use the quadratic terms of the unresolved tendencies
18 !> - reso: use the resolved dynamics alone
19 !
20 !---------------------------------------------------------------------------
21 
23  USE params, only: ndim,natm
24  USE stoch_params, only: q_ar,q_au,q_or,q_ou,mode
26  USE aotensor_def, only: aotensor
28  USE stoch_mod
29  IMPLICIT NONE
30 
31  PRIVATE
32 
33  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dwar,dwau,dwor,dwou !< Standard gaussian noise buffers
34 
35  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1,buf_f0,buf_f1 !< Integration buffers
36 
37  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: anoise !< Additive noise term
38 
39  TYPE(coolist), DIMENSION(:), ALLOCATABLE :: int_tensor !< Dummy tensor that will hold the tendencies tensor
40 
41  PUBLIC :: init_integrator, step
42 
43 CONTAINS
44 
45  !> Subroutine to initialize the integrator
46  !> @param force Parameter to force the mode of the integrator
47  SUBROUTINE init_integrator(force)
48  INTEGER :: AllocStat
49  CHARACTER*4, INTENT(IN), OPTIONAL :: force
50  CHARACTER*4 :: test
51 
52  ALLOCATE(buf_y1(0:ndim),buf_f0(0:ndim),buf_f1(0:ndim),stat=allocstat)
53  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
54 
55  ALLOCATE(anoise(0:ndim),stat=allocstat)
56  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
57 
58  ALLOCATE(dwar(0:ndim),dwau(0:ndim),dwor(0:ndim),dwou(0:ndim),stat=allocstat)
59  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
60 
61  ALLOCATE(int_tensor(ndim),stat=allocstat)
62  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
63 
64  dwar=0.d0
65  dwor=0.d0
66  dwau=0.d0
67  dwou=0.d0
68 
69  IF (PRESENT(force)) THEN
70  test=force
71  ELSE
72  test=mode
73  ENDIF
74 
75  SELECT CASE (test)
76  CASE('full')
78  CASE('ures')
80  CASE('qfst')
82  CASE('reso')
84  CASE DEFAULT
85  stop '*** MODE variable not properly defined ***'
86  END SELECT
87 
88  END SUBROUTINE init_integrator
89 
90  !> Routine computing the tendencies of the selected model
91  !> @param t Time at which the tendencies have to be computed. Actually not needed for autonomous systems.
92  !> @param y Point at which the tendencies have to be computed.
93  !> @param res vector to store the result.
94  !> @remark Note that it is NOT safe to pass `y` as a result buffer,
95  !> as this operation does multiple passes.
96  SUBROUTINE tendencies(t,y,res)
97  REAL(KIND=8), INTENT(IN) :: t
98  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
99  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
100  CALL sparse_mul3(int_tensor, y, y, res)
101  END SUBROUTINE tendencies
102 
103  !> Routine to perform a stochastic step of the selected dynamics (Heun algorithm).
104  !> The incremented time is returned.
105  !> @param y Initial point.
106  !> @param t Actual integration time
107  !> @param dt Integration timestep.
108  !> @param dtn Stochastic integration timestep (normally square-root of dt).
109  !> @param res Final point after the step.
110  !> @param tend Partial or full tendencies used to perform the step (used for debugging).
111  SUBROUTINE step(y,t,dt,dtn,res,tend)
112  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
113  REAL(KIND=8), INTENT(INOUT) :: t
114  REAL(KIND=8), INTENT(IN) :: dt,dtn
115  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res,tend
116 
117  CALL stoch_atm_res_vec(dwar)
119  CALL stoch_oc_res_vec(dwor)
122  CALL tendencies(t,y,buf_f0)
123  CALL sparse_mul3(int_tensor,y,y,tend)
124  buf_y1 = y+dt*buf_f0+anoise
126  tend=0.5*(tend+buf_f1)
127  CALL tendencies(t,buf_y1,buf_f1)
128  res=y+0.5*(buf_f0+buf_f1)*dt+anoise
129  t=t+dt
130  END SUBROUTINE step
131 
132 END MODULE rk2_stoch_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
subroutine, public copy_coo(src, dst)
Routine to copy a coolist.
Definition: tensor.f90:72
Statistics accumulators.
Definition: stat.f90:14
The resolved-unresolved components decomposition of the tensor.
Definition: dec_tensor.f90:16
type(coolist), dimension(:), allocatable, public byyy
Tensor holding the quadratic part of the unresolved tendencies involving unresolved variables...
Definition: dec_tensor.f90:48
real(kind=8) q_or
Oceanic resolved component noise amplitude.
subroutine, public init_integrator(force)
Subroutine to initialize the integrator.
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
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
real(kind=8), dimension(:), allocatable buf_y1
Tensor utility module.
Definition: tensor.f90:18
real(kind=8), dimension(:), allocatable dwou
Standard gaussian noise buffers.
real(kind=8), dimension(:), allocatable dwar
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
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
real(kind=8), dimension(:), allocatable anoise
Additive noise term.
type(coolist), dimension(:), allocatable, public ff_tensor
Tensor holding the part of the unresolved tensor involving only unresolved variables.
Definition: dec_tensor.f90:31
real(kind=8) q_ou
Oceanic unresolved component noise amplitude.
type(coolist), dimension(:), allocatable int_tensor
Dummy tensor that will hold the tendencies tensor.
real(kind=8), dimension(:), allocatable dwor
subroutine, public step(y, t, dt, dtn, res, tend)
Routine to perform a stochastic step of the selected dynamics (Heun algorithm). The incremented time ...
integer natm
Number of atmospheric basis functions.
Definition: params.f90:83
type(coolist), dimension(:), allocatable, public ss_tensor
Tensor holding the part of the resolved tensor involving only resolved variables. ...
Definition: dec_tensor.f90:33
real(kind=8), dimension(:), allocatable buf_f1
Integration buffers.
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_f0
type(coolist), dimension(:), allocatable, public aotensor
- Tensor representation of the tendencies.
character(len=4) mode
Stochastic mode parameter.
subroutine tendencies(t, y, res)
Routine computing the tendencies of the selected model.
real(kind=8) q_ar
Atmospheric resolved component noise amplitude.
Module with the stochastic rk2 integration routines.
real(kind=8), dimension(:), allocatable dwau