A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
tl_ad_tensor.f90
Go to the documentation of this file.
1 
2 ! tl_ad_tensor.f90
3 !
4 !> Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM.
5 !> Tensors definition module
6 !
7 !> @copyright
8 !> 2016 Lesley De Cruz & Jonathan Demaeyer.
9 !> See LICENSE.txt for license information.
10 !
11 !---------------------------------------------------------------------------!
12 !
13 !> @remark
14 !> The routines of this module should be called only after
15 !> params::init_params() and aotensor_def::init_aotensor() have been called !
16 !
17 !---------------------------------------------------------------------------!
18 
20 
21  !-----------------------------------------------------!
22  ! !
23  ! Preamble and variables declaration !
24  ! !
25  !-----------------------------------------------------!
26 
27  USE params, only:ndim
28  USE aotensor_def, only:aotensor
30  IMPLICIT NONE
31 
32  PRIVATE
33 
34  !> Epsilon to test equality with 0
35  REAL(KIND=8), PARAMETER :: real_eps = 2.2204460492503131e-16
36 
37  !> Vector used to count the tensor elements
38  INTEGER, DIMENSION(:), ALLOCATABLE :: count_elems
39 
40  !> Tensor representation of the Tangent Linear tendencies.
41  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: tltensor
42 
43  !> Tensor representation of the Adjoint tendencies.
44  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: adtensor
45 
47 
48  !-----------------------------------------------------!
49  ! !
50  ! End of preamble !
51  ! !
52  !-----------------------------------------------------!
53 
54 CONTAINS
55 
56  !-----------------------------------------------------!
57  ! !
58  ! Function declarations : !
59  ! !
60  ! These functions should be called only after !
61  ! init_params and init_aotensor !
62  ! !
63  !-----------------------------------------------------!
64 
65  !-----------------------------------------------------!
66  ! !
67  ! Jacobian functions !
68  ! !
69  !-----------------------------------------------------!
70 
71  !> Compute the Jacobian of MAOOAM in point ystar.
72  !> @param ystar array with variables in which the jacobian should be evaluated.
73  !> @return Jacobian in coolist-form (table of tuples {i,j,0,value})
74  FUNCTION jacobian(ystar)
75  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar
76  TYPE(coolist), DIMENSION(ndim) :: jacobian
77  CALL jsparse_mul(aotensor,ystar,jacobian)
78  END FUNCTION jacobian
79 
80  !> Compute the Jacobian of MAOOAM in point ystar.
81  !> @param ystar array with variables in which the jacobian should be evaluated.
82  !> @return Jacobian in matrix form
83  FUNCTION jacobian_mat(ystar)
84  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar
85  REAL(KIND=8), DIMENSION(ndim,ndim) :: jacobian_mat
86  CALL jsparse_mul_mat(aotensor,ystar,jacobian_mat)
87  END FUNCTION jacobian_mat
88 
89  !-----------------------------------------------------!
90  ! !
91  ! Tangent linear model functions !
92  ! !
93  !-----------------------------------------------------!
94 
95  !> Routine to initialize the TL tensor
96  SUBROUTINE init_tltensor
97  INTEGER :: i
98  INTEGER :: AllocStat
99  ALLOCATE(tltensor(ndim),count_elems(ndim), stat=allocstat)
100  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
101  count_elems=0
103 
104  DO i=1,ndim
105  ALLOCATE(tltensor(i)%elems(count_elems(i)), stat=allocstat)
106  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
107  END DO
108 
109  DEALLOCATE(count_elems, stat=allocstat)
110  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
111 
113 
114  CALL simplify(tltensor)
115 
116  END SUBROUTINE init_tltensor
117 
118  !> Routine to compute the TL tensor from the original MAOOAM one
119  !> @param func subroutine used to do the computation
120  SUBROUTINE compute_tltensor(func)
121  EXTERNAL :: func
122  INTERFACE
123  SUBROUTINE func(i,j,k,v)
124  INTEGER, INTENT(IN) :: i,j,k
125  REAL(KIND=8), INTENT(IN) :: v
126  END SUBROUTINE func
127  END INTERFACE
128  INTEGER :: i,j,k,n,ne
129  REAL(KIND=8) :: v
130  DO i=1,ndim
131  ne=aotensor(i)%nelems
132  DO n=1,ne
133  j=aotensor(i)%elems(n)%j
134  k=aotensor(i)%elems(n)%k
135  v=aotensor(i)%elems(n)%v
136  CALL func(i,j,k,v)
137  ENDDO
138  ENDDO
139  END SUBROUTINE compute_tltensor
140 
141  !> Subroutine used to count the number of TL tensor entries
142  !> @param i tensor \f$i\f$ index
143  !> @param j tensor \f$j\f$ index
144  !> @param k tensor \f$k\f$ index
145  !> @param v value that will be added
146  SUBROUTINE tl_add_count(i,j,k,v)
147  INTEGER, INTENT(IN) :: i,j,k
148  REAL(KIND=8), INTENT(IN) :: v
149  IF (abs(v) .ge. real_eps) THEN
150  IF (j /= 0) count_elems(i)=count_elems(i)+1
151  IF (k /= 0) count_elems(i)=count_elems(i)+1
152  ENDIF
153  END SUBROUTINE tl_add_count
154 
155  !> Subroutine used to compute the TL tensor entries
156  !> @param i tensor \f$i\f$ index
157  !> @param j tensor \f$j\f$ index
158  !> @param k tensor \f$k\f$ index
159  !> @param v value to add
160  SUBROUTINE tl_coeff(i,j,k,v)
161  INTEGER, INTENT(IN) :: i,j,k
162  REAL(KIND=8), INTENT(IN) :: v
163  INTEGER :: n
164  IF (.NOT. ALLOCATED(tltensor)) stop "*** tl_coeff routine : tensor not yet allocated ***"
165  IF (.NOT. ALLOCATED(tltensor(i)%elems)) stop "*** tl_coeff routine : tensor not yet allocated ***"
166  IF (abs(v) .ge. real_eps) THEN
167  IF (j /=0) THEN
168  n=(tltensor(i)%nelems)+1
169  tltensor(i)%elems(n)%j=j
170  tltensor(i)%elems(n)%k=k
171  tltensor(i)%elems(n)%v=v
172  tltensor(i)%nelems=n
173  END IF
174  IF (k /=0) THEN
175  n=(tltensor(i)%nelems)+1
176  tltensor(i)%elems(n)%j=k
177  tltensor(i)%elems(n)%k=j
178  tltensor(i)%elems(n)%v=v
179  tltensor(i)%nelems=n
180  END IF
181  END IF
182  END SUBROUTINE tl_coeff
183 
184  !-----------------------------------------------------!
185  ! !
186  ! Adjoint model functions !
187  ! !
188  !-----------------------------------------------------!
189 
190 
191  !> Routine to initialize the AD tensor
192  SUBROUTINE init_adtensor
193  INTEGER :: i
194  INTEGER :: AllocStat
195  ALLOCATE(adtensor(ndim),count_elems(ndim), stat=allocstat)
196  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
197  count_elems=0
199 
200  DO i=1,ndim
201  ALLOCATE(adtensor(i)%elems(count_elems(i)), stat=allocstat)
202  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
203  END DO
204 
205  DEALLOCATE(count_elems, stat=allocstat)
206  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
207 
209 
210  CALL simplify(adtensor)
211 
212  END SUBROUTINE init_adtensor
213 
214  !> Subroutine to compute the AD tensor from the original MAOOAM one
215  !> @param func subroutine used to do the computation
216  SUBROUTINE compute_adtensor(func)
217  EXTERNAL :: func
218  INTERFACE
219  SUBROUTINE func(i,j,k,v)
220  INTEGER, INTENT(IN) :: i,j,k
221  REAL(KIND=8), INTENT(IN) :: v
222  END SUBROUTINE func
223  END INTERFACE
224  INTEGER :: i,j,k,n,ne
225  REAL(KIND=8) :: v
226  DO i=1,ndim
227  ne=aotensor(i)%nelems
228  DO n=1,ne
229  j=aotensor(i)%elems(n)%j
230  k=aotensor(i)%elems(n)%k
231  v=aotensor(i)%elems(n)%v
232  CALL func(i,j,k,v)
233  ENDDO
234  ENDDO
235  END SUBROUTINE compute_adtensor
236 
237  !> Subroutine used to count the number of AD tensor entries
238  !> @param i tensor \f$i\f$ index
239  !> @param j tensor \f$j\f$ index
240  !> @param k tensor \f$k\f$ index
241  !> @param v value that will be added
242  SUBROUTINE ad_add_count(i,j,k,v)
243  INTEGER, INTENT(IN) :: i,j,k
244  REAL(KIND=8), INTENT(IN) :: v
245  IF ((abs(v) .ge. real_eps).AND.(i /= 0)) THEN
246  IF (k /= 0) count_elems(k)=count_elems(k)+1
247  IF (j /= 0) count_elems(j)=count_elems(j)+1
248  ENDIF
249  END SUBROUTINE ad_add_count
250 
251  ! Subroutine used to compute the AD tensor entries
252  !> @param i tensor \f$i\f$ index
253  !> @param j tensor \f$j\f$ index
254  !> @param k tensor \f$k\f$ index
255  !> @param v value to add
256  SUBROUTINE ad_coeff(i,j,k,v)
257  INTEGER, INTENT(IN) :: i,j,k
258  REAL(KIND=8), INTENT(IN) :: v
259  INTEGER :: n
260  IF (.NOT. ALLOCATED(adtensor)) stop "*** ad_coeff routine : tensor not yet allocated ***"
261  IF ((abs(v) .ge. real_eps).AND.(i /=0)) THEN
262  IF (k /=0) THEN
263  IF (.NOT. ALLOCATED(adtensor(k)%elems)) stop "*** ad_coeff routine : tensor not yet allocated ***"
264  n=(adtensor(k)%nelems)+1
265  adtensor(k)%elems(n)%j=i
266  adtensor(k)%elems(n)%k=j
267  adtensor(k)%elems(n)%v=v
268  adtensor(k)%nelems=n
269  END IF
270  IF (j /=0) THEN
271  IF (.NOT. ALLOCATED(adtensor(j)%elems)) stop "*** ad_coeff routine : tensor not yet allocated ***"
272  n=(adtensor(j)%nelems)+1
273  adtensor(j)%elems(n)%j=i
274  adtensor(j)%elems(n)%k=k
275  adtensor(j)%elems(n)%v=v
276  adtensor(j)%nelems=n
277  END IF
278  END IF
279  END SUBROUTINE ad_coeff
280 
281  !-----------------------------------------------------!
282  ! !
283  ! Adjoint model functions : Alternate method !
284  ! using the TL tensor !
285  ! (must be initialized !
286  ! before) !
287  ! !
288  !-----------------------------------------------------!
289 
290  !> Alternate method to initialize the AD tensor from the TL
291  !> tensor
292  !> @remark The #tltensor must be initialised before using this method.
293  SUBROUTINE init_adtensor_ref
294  INTEGER :: i
295  INTEGER :: AllocStat
296  ALLOCATE(adtensor(ndim),count_elems(ndim), stat=allocstat)
297  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
298  count_elems=0
300 
301  DO i=1,ndim
302  ALLOCATE(adtensor(i)%elems(count_elems(i)), stat=allocstat)
303  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
304  END DO
305 
306  DEALLOCATE(count_elems, stat=allocstat)
307  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
308 
310 
311  CALL simplify(adtensor)
312 
313  END SUBROUTINE init_adtensor_ref
314 
315  !> Alternate subroutine to compute the AD tensor from the TL one
316  !> @param func subroutine used to do the computation
317  SUBROUTINE compute_adtensor_ref(func)
318  EXTERNAL :: func
319  INTERFACE
320  SUBROUTINE func(i,j,k,v)
321  INTEGER, INTENT(IN) :: i,j,k
322  REAL(KIND=8), INTENT(IN) :: v
323  END SUBROUTINE func
324  END INTERFACE
325  INTEGER :: i,j,k,n,ne
326  REAL(KIND=8) :: v
327  IF (.NOT. ALLOCATED(tltensor)) stop "*** compute_adtensor_ref routine : tensor TL should be allocated first***"
328  DO i=1,ndim
329  ne=tltensor(i)%nelems
330  DO n=1,ne
331  j=tltensor(i)%elems(n)%j
332  k=tltensor(i)%elems(n)%k
333  v=tltensor(i)%elems(n)%v
334  CALL func(i,j,k,v)
335  ENDDO
336  ENDDO
337  END SUBROUTINE compute_adtensor_ref
338 
339  !> Alternate subroutine used to count the number of AD tensor entries
340  !> from the TL tensor
341  !> @param i tensor \f$i\f$ index
342  !> @param j tensor \f$j\f$ index
343  !> @param k tensor \f$k\f$ index
344  !> @param v value that will be added
345  SUBROUTINE ad_add_count_ref(i,j,k,v)
346  INTEGER, INTENT(IN) :: i,j,k
347  REAL(KIND=8), INTENT(IN) :: v
348  IF ((abs(v) .ge. real_eps).AND.(j /= 0)) count_elems(j)=count_elems(j)+1
349  END SUBROUTINE ad_add_count_ref
350 
351  !> Alternate subroutine used to compute the AD tensor entries from
352  !> the TL tensor
353  !> @param i tensor \f$i\f$ index
354  !> @param j tensor \f$j\f$ index
355  !> @param k tensor \f$k\f$ index
356  !> @param v value to add
357  SUBROUTINE ad_coeff_ref(i,j,k,v)
358  INTEGER, INTENT(IN) :: i,j,k
359  REAL(KIND=8), INTENT(IN) :: v
360  INTEGER :: n
361  IF (.NOT. ALLOCATED(adtensor)) stop "*** ad_coeff_ref routine : tensor not yet allocated ***"
362  IF ((abs(v) .ge. real_eps).AND.(j /=0)) THEN
363  IF (.NOT. ALLOCATED(adtensor(j)%elems)) stop "*** ad_coeff_ref routine : tensor not yet allocated ***"
364  n=(adtensor(j)%nelems)+1
365  adtensor(j)%elems(n)%j=i
366  adtensor(j)%elems(n)%k=k
367  adtensor(j)%elems(n)%v=v
368  adtensor(j)%nelems=n
369  END IF
370  END SUBROUTINE ad_coeff_ref
371 
372  !-----------------------------------------------------!
373  ! !
374  ! Tendencies !
375  ! !
376  !-----------------------------------------------------!
377 
378  !> Tendencies for the AD of MAOOAM in point ystar for perturbation deltay.
379  !> @param t time
380  !> @param ystar vector with the variables (current point in trajectory)
381  !> @param deltay vector with the perturbation of the variables at time t
382  !> @param buf vector (buffer) to store derivatives.
383  SUBROUTINE ad(t,ystar,deltay,buf)
384  REAL(KIND=8), INTENT(IN) :: t
385  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar,deltay
386  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: buf
387  CALL sparse_mul3(adtensor,deltay,ystar,buf)
388  END SUBROUTINE ad
389 
390  !> Tendencies for the TL of MAOOAM in point ystar for perturbation deltay.
391  !> @param t time
392  !> @param ystar vector with the variables (current point in trajectory)
393  !> @param deltay vector with the perturbation of the variables at time t
394  !> @param buf vector (buffer) to store derivatives.
395  SUBROUTINE tl(t,ystar,deltay,buf)
396  REAL(KIND=8), INTENT(IN) :: t
397  REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar,deltay
398  REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: buf
399  CALL sparse_mul3(tltensor,deltay,ystar,buf)
400  END SUBROUTINE tl
401 
402 END MODULE tl_ad_tensor
real(kind=8) function, dimension(ndim, ndim), public jacobian_mat(ystar)
Compute the Jacobian of MAOOAM in point ystar.
real(kind=8), parameter real_eps
Epsilon to test equality with 0.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
Statistics accumulators.
Definition: stat.f90:14
subroutine compute_tltensor(func)
Routine to compute the TL tensor from the original MAOOAM one.
subroutine, public init_tltensor
Routine to initialize the TL tensor.
subroutine, public simplify(tensor)
Routine to simplify a coolist (sparse tensor). For each index , it upper triangularize the matrix ...
Definition: tensor.f90:238
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
subroutine ad_add_count_ref(i, j, k, v)
Alternate subroutine used to count the number of AD tensor entries from the TL tensor.
subroutine ad_add_count(i, j, k, v)
Subroutine used to count the number of AD tensor entries.
subroutine ad_coeff_ref(i, j, k, v)
Alternate subroutine used to compute the AD tensor entries from the TL tensor.
subroutine tl_add_count(i, j, k, v)
Subroutine used to count the number of TL tensor entries.
subroutine, public jsparse_mul(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It's implemented slightly differentl...
Definition: tensor.f90:153
subroutine, public ad(t, ystar, deltay, buf)
Tendencies for the AD of MAOOAM in point ystar for perturbation deltay.
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
subroutine, public jsparse_mul_mat(coolist_ijk, arr_j, jcoo_ij)
Sparse multiplication of two tensors to determine the Jacobian: It's implemented slightly differentl...
Definition: tensor.f90:196
subroutine, public init_adtensor
Routine to initialize the AD tensor.
type(coolist) function, dimension(ndim) jacobian(ystar)
Compute the Jacobian of MAOOAM in point ystar.
subroutine, public init_adtensor_ref
Alternate method to initialize the AD tensor from the TL tensor.
subroutine tl_coeff(i, j, k, v)
Subroutine used to compute the TL tensor entries.
type(coolist), dimension(:), allocatable, public adtensor
Tensor representation of the Adjoint tendencies.
subroutine compute_adtensor(func)
Subroutine to compute the AD tensor from the original MAOOAM one.
subroutine ad_coeff(i, j, k, v)
type(coolist), dimension(:), allocatable, public tltensor
Tensor representation of the Tangent Linear tendencies.
subroutine, public tl(t, ystar, deltay, buf)
Tendencies for the TL of MAOOAM in point ystar for perturbation deltay.
The model parameters module.
Definition: params.f90:18
Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM. Tensors definition module...
type(coolist), dimension(:), allocatable, public aotensor
- Tensor representation of the tendencies.
subroutine compute_adtensor_ref(func)
Alternate subroutine to compute the AD tensor from the TL one.
integer, dimension(:), allocatable count_elems
Vector used to count the tensor elements.