A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
aotensor_def.f90
Go to the documentation of this file.
1 
2 ! aotensor_def.f90
3 !
4 !> The equation tensor for the coupled ocean-atmosphere model
5 !> with temperature which allows for an extensible set of modes
6 !> in the ocean and in the atmosphere.
7 !
8 !> @copyright
9 !> 2015 Lesley De Cruz & Jonathan Demaeyer.
10 !> See LICENSE.txt for license information.
11 !
12 !---------------------------------------------------------------------------!
13 !
14 !> @remark
15 !> Generated Fortran90/95 code
16 !> from aotensor.lua
17 !
18 !---------------------------------------------------------------------------!
19 
20 
22 
23  !-----------------------------------------------------!
24  ! !
25  ! Preamble and variables declaration !
26  ! !
27  !-----------------------------------------------------!
28 
29  USE params
30  USE inprod_analytic
31  USE tensor, only:coolist,simplify
32  IMPLICIT NONE
33 
34  PRIVATE
35 
36  !> Vector used to count the tensor elements
37  INTEGER, DIMENSION(:), ALLOCATABLE :: count_elems
38 
39  !> Epsilon to test equality with 0
40  REAL(KIND=8), PARAMETER :: real_eps = 2.2204460492503131e-16
41 
42  PUBLIC :: init_aotensor
43 
44  !> \f$\mathcal{T}_{i,j,k}\f$ - Tensor representation of the tendencies.
45  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: aotensor
46 
47 
48  !-----------------------------------------------------!
49  ! !
50  ! End of preamble !
51  ! !
52  !-----------------------------------------------------!
53 
54 CONTAINS
55 
56  !-----------------------------------------------------!
57  ! !
58  ! Function declarations !
59  ! !
60  !-----------------------------------------------------!
61 
62  !> Translate the \f$\psi_{a,i}\f$ coefficients into effective coordinates
63  FUNCTION psi(i)
64  INTEGER :: i,psi
65  psi = i
66  END FUNCTION psi
67 
68  !> Translate the \f$\theta_{a,i}\f$ coefficients into effective coordinates
69  FUNCTION theta(i)
70  INTEGER :: i,theta
71  theta = i + natm
72  END FUNCTION theta
73 
74  !> Translate the \f$\psi_{o,i}\f$ coefficients into effective coordinates
75  FUNCTION a(i)
76  INTEGER :: i,A
77  a = i + 2 * natm
78  END FUNCTION a
79 
80  !> Translate the \f$\delta T_{o,i}\f$ coefficients into effective coordinates
81  FUNCTION t(i)
82  INTEGER :: i,T
83  t = i + 2 * natm + noc
84  END FUNCTION t
85 
86  !> Kronecker delta function
87  FUNCTION kdelta(i,j)
88  INTEGER :: i,j,kdelta
89  kdelta=0
90  IF (i == j) kdelta = 1
91  END FUNCTION kdelta
92 
93  !> Subroutine to add element in the #aotensor \f$\mathcal{T}_{i,j,k}\f$ structure.
94  !> @param i tensor \f$i\f$ index
95  !> @param j tensor \f$j\f$ index
96  !> @param k tensor \f$k\f$ index
97  !> @param v value to add
98  SUBROUTINE coeff(i,j,k,v)
99  INTEGER, INTENT(IN) :: i,j,k
100  REAL(KIND=8), INTENT(IN) :: v
101  INTEGER :: n
102  IF (.NOT. ALLOCATED(aotensor)) stop "*** coeff routine : tensor not yet allocated ***"
103  IF (.NOT. ALLOCATED(aotensor(i)%elems)) stop "*** coeff routine : tensor not yet allocated ***"
104  IF (abs(v) .ge. real_eps) THEN
105  n=(aotensor(i)%nelems)+1
106  IF (j .LE. k) THEN
107  aotensor(i)%elems(n)%j=j
108  aotensor(i)%elems(n)%k=k
109  ELSE
110  aotensor(i)%elems(n)%j=k
111  aotensor(i)%elems(n)%k=j
112  END IF
113  aotensor(i)%elems(n)%v=v
114  aotensor(i)%nelems=n
115  END IF
116  END SUBROUTINE coeff
117 
118  !> Subroutine to count the elements of the #aotensor \f$\mathcal{T}_{i,j,k}\f$. Add +1 to count_elems(i) for each value that is added to the tensor i-th component.
119  !> @param i tensor \f$i\f$ index
120  !> @param j tensor \f$j\f$ index
121  !> @param k tensor \f$k\f$ index
122  !> @param v value that will be added
123  SUBROUTINE add_count(i,j,k,v)
124  INTEGER, INTENT(IN) :: i,j,k
125  REAL(KIND=8), INTENT(IN) :: v
126  IF (abs(v) .ge. real_eps) count_elems(i)=count_elems(i)+1
127  END SUBROUTINE add_count
128 
129  !> Subroutine to compute the tensor #aotensor
130  !> @param func External function to be used
131  SUBROUTINE compute_aotensor(func)
132  EXTERNAL :: func
133  INTERFACE
134  SUBROUTINE func(i,j,k,v)
135  INTEGER, INTENT(IN) :: i,j,k
136  REAL(KIND=8), INTENT(IN) :: v
137  END SUBROUTINE func
138  END INTERFACE
139  INTEGER :: i,j,k
140  CALL func(theta(1),0,0,(cpa / (1 - atmos%a(1,1) * sig0)))
141  DO i = 1, natm
142  DO j = 1, natm
143  CALL func(psi(i),psi(j),0,-(((atmos%c(i,j) * betp) / atmos%a(i,i))) -&
144  &(kd * kdelta(i,j)) / 2 + atmos%a(i,j)*nuap)
145  CALL func(theta(i),psi(j),0,(atmos%a(i,j) * kd * sig0) / (-2 + 2 * atmos%a(i,i) * sig0))
146  CALL func(psi(i),theta(j),0,(kd * kdelta(i,j)) / 2)
147  CALL func(theta(i),theta(j),0,(-((sig0 * (2. * atmos%c(i,j) * betp +&
148  & atmos%a(i,j) * (kd + 4. * kdp)))) + 2. * (lsbpa + sc * lpa) &
149  &* kdelta(i,j)) / (-2. + 2. * atmos%a(i,i) * sig0))
150  DO k = 1, natm
151  CALL func(psi(i),psi(j),psi(k),-((atmos%b(i,j,k) / atmos%a(i,i))))
152  CALL func(psi(i),theta(j),theta(k),-((atmos%b(i,j,k) / atmos%a(i,i))))
153  CALL func(theta(i),psi(j),theta(k),(atmos%g(i,j,k) -&
154  & atmos%b(i,j,k) * sig0) / (-1 + atmos%a(i,i) *&
155  & sig0))
156  CALL func(theta(i),theta(j),psi(k),(atmos%b(i,j,k) * sig0) / (1 - atmos%a(i,i) * sig0))
157  END DO
158  END DO
159  DO j = 1, noc
160  CALL func(psi(i),a(j),0,kd * atmos%d(i,j) / (2 * atmos%a(i,i)))
161  CALL func(theta(i),a(j),0,kd * (atmos%d(i,j) * sig0) / (2 - 2 * atmos%a(i,i) * sig0))
162  CALL func(theta(i),t(j),0,atmos%s(i,j) * (2 * lsbpo + lpa) / (2 - 2 * atmos%a(i,i) * sig0))
163  END DO
164  END DO
165  DO i = 1, noc
166  DO j = 1, natm
167  CALL func(a(i),psi(j),0,ocean%K(i,j) * dp / (ocean%M(i,i) + g))
168  CALL func(a(i),theta(j),0,-(ocean%K(i,j)) * dp / (ocean%M(i,i) + g))
169  END DO
170  DO j = 1, noc
171  CALL func(a(i),a(j),0,-((ocean%N(i,j) * betp + ocean%M(i,i) * (rp + dp) * kdelta(i,j)&
172  & - ocean%M(i,j)**2*nuop)) / (ocean%M(i,i) + g))
173  DO k = 1, noc
174  CALL func(a(i),a(j),a(k),-(ocean%C(i,j,k)) / (ocean%M(i,i) + g))
175  END DO
176  END DO
177  END DO
178  DO i = 1, noc
179  CALL func(t(i),0,0,cpo * ocean%W(i,1))
180  DO j = 1, natm
181  CALL func(t(i),theta(j),0,ocean%W(i,j) * (2 * sc * lpo + sbpa))
182  END DO
183  DO j = 1, noc
184  CALL func(t(i),t(j),0,-((lpo + sbpo)) * kdelta(i,j))
185  DO k = 1, noc
186  CALL func(t(i),a(j),t(k),-(ocean%O(i,j,k)))
187  END DO
188  END DO
189  END DO
190  END SUBROUTINE compute_aotensor
191 
192  !-----------------------------------------------------!
193  ! !
194  ! Initialisation routine !
195  ! !
196  !-----------------------------------------------------!
197 
198  !> Subroutine to initialise the #aotensor tensor
199  !> @remark This procedure will also call params::init_params() and inprod_analytic::init_inprod() .
200  !> It will finally call inprod_analytic::deallocate_inprod() to remove the inner products, which are not needed
201  !> anymore at this point.
202  SUBROUTINE init_aotensor
203  INTEGER :: i
204  INTEGER :: AllocStat
205 
206  CALL init_params ! Iniatialise the parameter
207 
208  CALL init_inprod ! Initialise the inner product tensors
209 
210  ALLOCATE(aotensor(ndim),count_elems(ndim), stat=allocstat)
211  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
212  count_elems=0
213 
215 
216  DO i=1,ndim
217  ALLOCATE(aotensor(i)%elems(count_elems(i)), stat=allocstat)
218  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
219  END DO
220 
221  DEALLOCATE(count_elems, stat=allocstat)
222  IF (allocstat /= 0) stop "*** Deallocation problem ! ***"
223 
224  CALL compute_aotensor(coeff)
225 
226  CALL simplify(aotensor)
227 
228  END SUBROUTINE init_aotensor
229 END MODULE aotensor_def
230 
231 
232 
real(kind=8) nuap
Non-dimensional dissipation in the atmosphere.
Definition: params.f90:72
real(kind=8) cpa
- Non-dimensional constant short-wave radiation of the atmosphere.
Definition: params.f90:58
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
integer noc
Number of oceanic basis functions.
Definition: params.f90:84
subroutine, public init_inprod
Initialisation of the inner product.
real(kind=8) sig0
- Non-dimensional static stability of the atmosphere.
Definition: params.f90:27
real(kind=8) betp
- Non-dimensional beta parameter
Definition: params.f90:67
Inner products between the truncated set of basis functions for the ocean and atmosphere streamfuncti...
Statistics accumulators.
Definition: stat.f90:14
subroutine coeff(i, j, k, v)
Subroutine to add element in the aotensor structure.
integer function a(i)
Translate the coefficients into effective coordinates.
subroutine, public simplify(tensor)
Routine to simplify a coolist (sparse tensor). For each index , it upper triangularize the matrix ...
Definition: tensor.f90:238
real(kind=8) nuop
Non-dimensional dissipation in the ocean.
Definition: params.f90:73
real(kind=8) kdp
- Non-dimensional internal atmospheric friction coefficient.
Definition: params.f90:54
real(kind=8) kd
- Non-dimensional bottom atmospheric friction coefficient.
Definition: params.f90:53
The equation tensor for the coupled ocean-atmosphere model with temperature which allows for an exten...
integer function kdelta(i, j)
Kronecker delta function.
subroutine compute_aotensor(func)
Subroutine to compute the tensor aotensor.
Tensor utility module.
Definition: tensor.f90:18
real(kind=8) sc
Ratio of surface to atmosphere temperature.
Definition: params.f90:65
real(kind=8) g
Definition: params.f90:50
subroutine add_count(i, j, k, v)
Subroutine to count the elements of the aotensor . Add +1 to count_elems(i) for each value that is ad...
real(kind=8) sbpo
- Long wave radiation lost by ocean to atmosphere & space.
Definition: params.f90:60
real(kind=8) dp
- Non-dimensional mechanical coupling parameter between the ocean and the atmosphere.
Definition: params.f90:52
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
integer function psi(i)
Translate the coefficients into effective coordinates.
real(kind=8), parameter real_eps
Epsilon to test equality with 0.
integer natm
Number of atmospheric basis functions.
Definition: params.f90:83
type(atm_tensors), public atmos
Atmospheric tensors.
real(kind=8) lpo
- Non-dimensional sensible + turbulent heat exchange from ocean to atmosphere.
Definition: params.f90:57
integer function t(i)
Translate the coefficients into effective coordinates.
real(kind=8) sbpa
- Long wave radiation from atmosphere absorbed by ocean.
Definition: params.f90:61
integer, dimension(:), allocatable count_elems
Vector used to count the tensor elements.
The model parameters module.
Definition: params.f90:18
subroutine init_params
Parameters initialisation routine.
Definition: params.f90:138
real(kind=8) lsbpa
- Long wave radiation lost by atmosphere to space & ocean.
Definition: params.f90:63
real(kind=8) lpa
- Non-dimensional sensible + turbulent heat exchange from atmosphere to ocean.
Definition: params.f90:59
type(coolist), dimension(:), allocatable, public aotensor
- Tensor representation of the tendencies.
real(kind=8) rp
- Frictional coefficient at the bottom of the ocean.
Definition: params.f90:51
integer function theta(i)
Translate the coefficients into effective coordinates.
type(ocean_tensors), public ocean
Oceanic tensors.
real(kind=8) cpo
- Non-dimensional constant short-wave radiation of the ocean.
Definition: params.f90:56
real(kind=8) lsbpo
- Long wave radiation from ocean absorbed by atmosphere.
Definition: params.f90:62
subroutine, public init_aotensor
Subroutine to initialise the aotensor tensor.