A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
rk2_WL_integrator.f90
Go to the documentation of this file.
1 
2 ! rk2_WL_integrator.f90
3 !
4 !> Module with the WL 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,natm
21  USE wl_tensor
22  USE aotensor_def, only:aotensor
23  USE mar, only: init_mar,mar_step,ms
24  USE stoch_mod
25  USE memory
27  IMPLICIT NONE
28 
29  PRIVATE
30 
31  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1,buf_f0,buf_f1 !< Integration buffers
32  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_m2,buf_m1,buf_m3,buf_m,buf_m3s !< Dummy buffers holding the terms /f$M_i\f$ of the parameterization
33  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: anoise !< Additive noise term
34  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dwar,dwau,dwor,dwou !< Standard gaussian noise buffers
35 
36  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: x1 !< Buffer holding the subsequent states of the first MAR
37  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: x2 !< Buffer holding the subsequent states of the second MAR
38 
39  PUBLIC :: init_integrator, step, full_step
40 
41 CONTAINS
42  !> Subroutine that initialize the MARs, the memory unit and the integration buffers
43  SUBROUTINE init_integrator
44  INTEGER :: AllocStat,i
45 
47 
48  print*, 'Initializing the integrator ...'
49 
50  IF (mode.ne.'ures') THEN
51  print*, '*** Mode set to ',mode,' in stoch_params.nml ***'
52  print*, '*** WL configuration only support unresolved mode ***'
53  stop "*** Please change to 'ures' and perform the configuration again ! ***"
54  ENDIF
55 
56  ALLOCATE(buf_y1(0:ndim),buf_f0(0:ndim),buf_f1(0:ndim),stat=allocstat)
57  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
58 
59  ALLOCATE(buf_m1(0:ndim), buf_m2(0:ndim), buf_m3(0:ndim), buf_m(0:ndim), buf_m3s(0:ndim), stat=allocstat)
60  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
61 
62  ALLOCATE(dwar(0:ndim),dwau(0:ndim),dwor(0:ndim),dwou(0:ndim),stat=allocstat)
63  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
64 
65  ALLOCATE(anoise(0:ndim), stat=allocstat)
66  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
67 
68  buf_y1=0.d0
69  buf_f1=0.d0
70  buf_f0=0.d0
71 
72  dwar=0.d0
73  dwor=0.d0
74  dwau=0.d0
75  dwou=0.d0
76 
77  buf_m1=0.d0
78  buf_m2=0.d0
79  buf_m3=0.d0
80  buf_m3s=0.d0
81  buf_m=0.d0
82 
83  print*, 'Initializing the MARs ...'
84 
85  CALL init_mar
86 
87  ALLOCATE(x1(0:ndim,ms), x2(0:ndim,ms), stat=allocstat)
88 
89  x1=0.d0
90  DO i=1,50000
91  CALL mar_step(x1)
92  ENDDO
93 
94  x2=0.d0
95  DO i=1,50000
96  CALL mar_step(x2)
97  ENDDO
98 
99  CALL init_memory
100 
101  END SUBROUTINE init_integrator
102 
103  !> Routine to compute the \f$M_1\f$ term
104  !> @param y Present state of the WL system
105  SUBROUTINE compute_m1(y)
106  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
107  buf_m1=0.d0
108  IF (m12def) CALL sparse_mul2_k(m12, y, buf_m1)
110  END SUBROUTINE compute_m1
111 
112  !> Routine to compute the \f$M_2\f$ term
113  !> @param y Present state of the WL system
114  SUBROUTINE compute_m2(y)
115  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
116  buf_m=0.d0
117  buf_m2=0.d0
118  IF (m21def) CALL sparse_mul3(m21, y, x1(0:ndim,1), buf_m)
119  IF (m22def) CALL sparse_mul3(m22, x2(0:ndim,1), x2(0:ndim,1), buf_m2)
121  END SUBROUTINE compute_m2
122 
123  !> Routine to perform an integration step (Heun algorithm) of the WL system.
124  !> The incremented time is returned.
125  !> @param y Initial point.
126  !> @param t Actual integration time
127  !> @param dt Integration timestep.
128  !> @param dtn Stochastic integration timestep (normally square-root of dt).
129  !> @param res Final point after the step.
130  !> @param tend Partial or full tendencies used to perform the step (used for debugging).
131  SUBROUTINE step(y,t,dt,dtn,res,tend)
132  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
133  REAL(KIND=8), INTENT(INOUT) :: t
134  REAL(KIND=8), INTENT(IN) :: dt,dtn
135  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res,tend
136  INTEGER :: i
137 
138  IF (mod(t,muti)<dt) THEN
139  CALL compute_m3(y,dts,dtsn,.true.,.true.,.true.,muti/2,buf_m3s)
141  DO i=1,1
142  CALL compute_m3(y,dts,dtsn,.false.,.true.,.true.,muti/2,buf_m3s)
144  ENDDO
145  !DO i=1,2
146  ! CALL compute_M3(y,dts,dtsn,.false.,.true.,.true.,muti/2,buf_M3s)
147  ! buf_M3=buf_M3+buf_M3s
148  !ENDDO
149  buf_m3=buf_m3/2
150  ENDIF
151 
152 
153  CALL stoch_atm_res_vec(dwar)
154  CALL stoch_oc_res_vec(dwor)
155  anoise=(q_ar*dwar+q_or*dwor)*dtn
156 
157  CALL tendencies(t,y,buf_f0)
158  CALL mar_step(x1)
159  CALL mar_step(x2)
160  CALL compute_m1(y)
161  CALL compute_m2(y)
163  buf_y1 = y+dt*buf_f0+anoise
164 
165  CALL tendencies(t+dt,buf_y1,buf_f1)
166  CALL compute_m1(buf_y1)
167  CALL compute_m2(buf_y1)
168  !IF (mod(t,muti)<dt) CALL compute_M3(buf_y1,dts,dtsn,.false.,.true.,buf_M3)
169 
171  res=y+dt*buf_f0+anoise
172 
173  tend=buf_m3
174  t=t+dt
175 
176  END SUBROUTINE step
177 
178  !> Routine to perform an integration step (Heun algorithm) of the full stochastic system. The incremented time is returned.
179  !> @param y Initial point.
180  !> @param t Actual integration time
181  !> @param dt Integration timestep.
182  !> @param dtn Stochastoc integration timestep (normally square-root of dt).
183  !> @param res Final point after the step.
184  SUBROUTINE full_step(y,t,dt,dtn,res)
185  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y
186  REAL(KIND=8), INTENT(INOUT) :: t
187  REAL(KIND=8), INTENT(IN) :: dt,dtn
188  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res
189  CALL stoch_atm_res_vec(dwar)
191  CALL stoch_oc_res_vec(dwor)
194  CALL sparse_mul3(aotensor,y,y,buf_f0)
195  buf_y1 = y+dt*buf_f0+anoise
197  res=y+0.5*(buf_f0+buf_f1)*dt+anoise
198  t=t+dt
199  END SUBROUTINE full_step
200 
201 END MODULE rk2_wl_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.
subroutine, public mar_step(x)
Routine to generate one step of the MAR.
Definition: MAR.f90:93
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
real(kind=8), dimension(:,:), allocatable x1
Buffer holding the subsequent states of the first MAR.
real(kind=8) muti
Memory update time interval.
real(kind=8), dimension(:), allocatable buf_m
Dummy vector.
Definition: memory.f90:31
Statistics accumulators.
Definition: stat.f90:14
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.
type(coolist), dimension(:), allocatable, public m12
Second component of the M1 term.
Definition: WL_tensor.f90:43
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
procedure(ss_step), pointer step
Procedural pointer pointing on the resolved dynamics step routine.
Definition: memory.f90:36
real(kind=8), dimension(:), allocatable dwor
Standard gaussian noise buffers.
logical, public m21def
Definition: WL_tensor.f90:75
The equation tensor for the coupled ocean-atmosphere model with temperature which allows for an exten...
real(kind=8), dimension(:,:), allocatable x2
Buffer holding the subsequent states of the second MAR.
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...
real(kind=8), dimension(:), allocatable buf_m3
Dummy vector to store the integrand.
Definition: memory.f90:32
logical, public m12def
Definition: WL_tensor.f90:75
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 m22
Second tensor of the M2 term.
Definition: WL_tensor.f90:49
real(kind=8), dimension(:), allocatable dwar
real(kind=8), dimension(:), allocatable buf_m3s
Dummy buffers holding the terms /f$M_i.
subroutine, public init_integrator
Subroutine that initialize the MARs, the memory unit and the integration buffers. ...
real(kind=8) q_ou
Oceanic unresolved component noise amplitude.
The WL tensors used to integrate the model.
Definition: WL_tensor.f90:18
subroutine compute_m2(y)
Routine to compute the term.
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.
integer natm
Number of atmospheric basis functions.
Definition: params.f90:83
real(kind=8), dimension(:), allocatable buf_m2
real(kind=8), dimension(:), allocatable buf_y1
real(kind=8) dtsn
Square root of the intrisic resolved dynamics time step.
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
subroutine, public init_memory
Subroutine to initialise the memory.
Definition: memory.f90:45
Module with the WL rk2 integration routines.
The model parameters module.
Definition: params.f90:18
Utility module containing the stochastic related routines.
Definition: stoch_mod.f90:15
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
real(kind=8), dimension(:), allocatable anoise
Additive noise term.
logical, public m22def
Definition: WL_tensor.f90:75
Multidimensional Autoregressive module to generate the correlation for the WL parameterization.
Definition: MAR.f90:19
Module that compute the memory term of the WL parameterization.
Definition: memory.f90:16
real(kind=8), dimension(:), allocatable buf_m1
type(coolist), dimension(:), allocatable, public aotensor
- Tensor representation of the tendencies.
character(len=4) mode
Stochastic mode parameter.
real(kind=8), dimension(:), allocatable dwou
Standard gaussian noise buffers.
real(kind=8) q_ar
Atmospheric resolved component noise amplitude.
type(coolist), dimension(:), allocatable, public m21
First tensor of the M2 term.
Definition: WL_tensor.f90:48
integer, public ms
order of the MAR
Definition: MAR.f90:36
subroutine, public init_mar
Subroutine to initialise the MAR.
Definition: MAR.f90:45
real(kind=8), dimension(:), allocatable dwau
real(kind=8), dimension(:), allocatable, public m1tot
Total vector.
Definition: WL_tensor.f90:45
real(kind=8) dts
Intrisic resolved dynamics time step.
subroutine compute_m1(y)
Routine to compute the term.