A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
rk4_integrator.f90
Go to the documentation of this file.
1 
2 ! integrator.f90
3 !
4 !> Module with the RK4 integration routines.
5 !
6 !> @copyright
7 !> 2015 Lesley De Cruz & Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !> This module actually contains the RK4 algorithm routines.
14 !> The user can modify it according to its preferred integration scheme.
15 !> For higher-order schemes, additional buffers will probably have to be defined.
16 !
17 !---------------------------------------------------------------------------
18 
19 MODULE integrator
20  USE params, only: ndim
21  USE tensor, only: sparse_mul3
22  USE aotensor_def, only: aotensor
23  IMPLICIT NONE
24 
25  PRIVATE
26 
27  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1 !< Buffer to hold the intermediate position (Heun algorithm)
28  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_ka !< Buffer A to hold tendencies
29  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_kb !< Buffer B to hold tendencies
30 
31  PUBLIC :: init_integrator, step
32 
33 CONTAINS
34 
35  !> Routine to initialise the integration buffers.
36  SUBROUTINE init_integrator
37  INTEGER :: allocstat
38  ALLOCATE(buf_y1(0:ndim),buf_ka(0:ndim),buf_kb(0:ndim) ,stat=allocstat)
39  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
40  END SUBROUTINE init_integrator
41 
42  !> Routine computing the tendencies of the model
43  !> @param t Time at which the tendencies have to be computed. Actually not needed for autonomous systems.
44  !> @param y Point at which the tendencies have to be computed.
45  !> @param res vector to store the result.
46  !> @remark Note that it is NOT safe to pass `y` as a result buffer,
47  !> as this operation does multiple passes.
48  SUBROUTINE tendencies(t,y,res)
49  REAL(KIND=8), INTENT(IN) :: t
50  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
51  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
52  CALL sparse_mul3(aotensor, y, y, res)
53  END SUBROUTINE tendencies
54 
55  !> Routine to perform an integration step (RK4 algorithm). The incremented time is returned.
56  !> @param y Initial point.
57  !> @param t Actual integration time
58  !> @param dt Integration timestep.
59  !> @param res Final point after the step.
60  SUBROUTINE step(y,t,dt,res)
61  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
62  REAL(KIND=8), INTENT(INOUT) :: t
63  REAL(KIND=8), INTENT(IN) :: dt
64  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
65 
66  CALL tendencies(t,y,buf_ka)
67  buf_y1 = y + 0.5*dt*buf_ka
68 
69  CALL tendencies(t+0.5*dt,buf_y1,buf_kb)
70  buf_y1 = y + 0.5*dt*buf_kb
71  buf_ka = buf_ka + 2*buf_kb
72 
73  CALL tendencies(t+0.5*dt,buf_y1,buf_kb)
74  buf_y1 = y + dt*buf_kb
75  buf_ka = buf_ka + 2*buf_kb
76 
77  CALL tendencies(t+dt,buf_y1,buf_kb)
78  buf_ka = buf_ka + buf_kb
79 
80  t=t+dt
81  res=y+buf_ka*dt/6
82  END SUBROUTINE step
83 END MODULE integrator
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public step(y, t, dt, res)
Routine to perform an integration step (Heun algorithm). The incremented time is returned.
Statistics accumulators.
Definition: stat.f90:14
real(kind=8), dimension(:), allocatable buf_kb
Buffer B to hold tendencies.
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
Tensor utility module.
Definition: tensor.f90:18
real(kind=8), dimension(:), allocatable buf_ka
Buffer A to hold tendencies.
Module with the integration routines.
subroutine tendencies(t, y, res)
Routine computing the tendencies of the model.
The model parameters module.
Definition: params.f90:18
real(kind=8), dimension(:), allocatable buf_y1
Buffer to hold the intermediate position (Heun algorithm)
type(coolist), dimension(:), allocatable, public aotensor
- Tensor representation of the tendencies.
subroutine, public init_integrator
Routine to initialise the integration buffers.