FMS  2025.04
Flexible Modeling System
cloud_interpolator.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 #define _FLATTEN(A) reshape((A), (/size((A))/) )
19 
20 !> @defgroup cloud_interpolator_mod cloud_interpolator_mod
21 !> @ingroup drifters
22 !! @brief Cloud interpolation routines for use in @ref drifters_mod
23 
24 !> @addtogroup cloud_interpolator_mod
25 !> @{
26 MODULE cloud_interpolator_mod
27 #ifdef use_drifters
28  implicit none
29  private
30 
31  public :: cld_ntrp_linear_cell_interp, cld_ntrp_locate_cell, cld_ntrp_get_cell_values
32  public :: cld_ntrp_expand_index, cld_ntrp_contract_indices
33 
34 ! Include variable "version" to be written to log file.
35 #include<file_version.h>
36 real, parameter :: tol = 10.0*epsilon(1.)
37 
38 CONTAINS
39 
40 !...............................................................................
41 !> Get expanded list of indices from contracted index
42 !> @param Ic contracted index
43 !> @param[out] ie(:) expanded list of indices
44 !> @param[out] ier error flag, non zero if operation unsuccessful
45 pure subroutine cld_ntrp_expand_index(Ic, ie, ier)
46  integer, intent(in) :: Ic
47  integer, intent(out) :: ie(:)
48  integer, intent(out) :: ier
49 
50  integer j, nd
51 
52  ier = 0
53  nd = size(ie) ! dimension
54 
55  if(ic >= 2**nd) then
56  ie = -1
57  ier = 1 ! error
58  return
59  endif
60 
61  do j = 1, nd
62  ie(j) = mod(ic/2**(j-1), 2)
63  end do
64 
65  end subroutine cld_ntrp_expand_index
66 
67 !...............................................................................
68 !> Contract list of indices to an single integer
69 !> @param ie(:) expanded list of indices
70 !> @param[out] Ic contracted index
71 !> @param[out] ier error flag, non zero if operation unsuccessful
72 pure subroutine cld_ntrp_contract_indices(ie, Ic, ier)
73  integer, intent(in) :: ie(:)
74  integer, intent(out) :: Ic
75  integer, intent(out) :: ier
76 
77  integer j, nd
78 
79  ier = 0
80  nd = size(ie) ! dimension
81 
82  ic = ie(nd)
83  do j = nd-1, 1, -1
84  ic = ic * 2
85  ic = ic + ie(j)
86  end do
87 
88  if(ic >= 2**nd) ier = 1
89 
90  end subroutine cld_ntrp_contract_indices
91 
92 
93 !...............................................................................
94 !...............................................................................
95 !> Cloud interpolation for linear cells
96 !> @param fvals values at the cell nodes
97 !> @param ts normalized [0,1]^nd cell coordinates
98 !> @param[out] interpolated value
99 !> @param[out] error flag, non zero if unsucessful
100 pure subroutine cld_ntrp_linear_cell_interp(fvals, ts, f, ier)
101  real, intent(in) :: fvals(0:)
102  real, intent(in) :: ts(:)
103  real, intent(out):: f
104  integer, intent(out) :: ier
105 
106  integer j, nd, Ic, iflag
107  integer ie(size(fvals))
108  real basis
109 
110  ier = 0
111  f = 0.
112  nd = size(ts)
113  if(size(fvals) /= 2**nd) then
114  ier = 1
115  return
116  endif
117 
118  do ic = 0, 2**nd - 1
119  basis = 1.
120  call cld_ntrp_expand_index(ic, ie, iflag)
121  do j = 1, nd
122  basis = basis * ( (1.0-real(ie(j)))*(1.0-ts(j)) + real(ie(j))*ts(j) )
123  end do
124  f = f + fvals(ic)*basis
125  end do
126 
127  end subroutine cld_ntrp_linear_cell_interp
128 
129 !...............................................................................
130 !...............................................................................
131 pure subroutine cld_ntrp_locate_cell(axis, x, index, ier)
132  real, intent(in) :: axis(:) !< axis
133  real, intent(in) :: x !< abscissae
134  integer, intent(out) :: index !< lower-left corner index
135  integer, intent(out) :: ier !< error flag (0=ok)
136 
137  logical down
138  integer n, index1, is
139  real axis_1, axis_n, axis_min, axis_max
140  ier = 0
141  index = -1
142  down = .false.
143  n = size(axis)
144  if(n < 2) then
145  ier = 3
146  return
147  endif
148  axis_1 = axis(1)
149  axis_n = axis(n)
150  axis_min = axis_1
151  axis_max = axis_n
152  if(axis_1 > axis_n) then
153  down = .true.
154  axis_min = axis_n
155  axis_max = axis_1
156  endif
157 
158  if(x < axis_min-tol) then
159  ier = 1
160  return
161  endif
162  if(x > axis_max+tol) then
163  ier = 2
164  return
165  endif
166 
167  index = floor(real(n-1)*(x - axis_1)/(axis_n-axis_1)) + 1
168  index = min(n-1, index)
169  index1 = index+1
170 
171  if(.NOT. down) then
172  if(axis(index) <= x+tol) then
173  if(x <= axis(index1)+tol) then
174  ! axis is uniform, or nearly so. Done!
175  return
176  else
177  ! increase index
178  is = index+1
179  do index = is, n-1
180  index1 = index+1
181  if(axis(index1) >= x-tol) return
182  enddo
183  endif
184  else
185  ! decrease index
186  is = index - 1
187  do index = is, 1, -1
188  if(axis(index) <= x+tol) return
189  enddo
190  endif
191  else
192  ! axis is pointing down
193  if(axis(index) >= x-tol) then
194  if(x >= axis(index1)-tol) then
195  ! axis is uniform, or nearly so. Done!
196  return
197  else
198  ! increase index
199  is = index + 1
200  do index = is, n-1
201  index1 = index+1
202  if(axis(index1) <= x+tol) return
203  enddo
204  endif
205  else
206  ! decrease index
207  is = index - 1
208  do index = is, 1, -1
209  if(axis(index) >= x-tol) return
210  enddo
211  endif
212  endif
213 
214  end subroutine cld_ntrp_locate_cell
215 
216 !...............................................................................
217 !...............................................................................
218 pure subroutine cld_ntrp_get_flat_index(nsizes, indices, flat_index, ier)
219  integer, intent(in) :: nsizes(:) !< size of array along each axis
220  integer, intent(in) :: indices(:) !< cell indices
221  integer, intent(out) :: flat_index !< index into flattened array
222  integer, intent(out) :: ier !< error flag (0=ok)
223 
224  integer nd, id
225 
226  ier = 0
227  flat_index = -1
228  nd = size(nsizes)
229  if(nd /= size(indices)) then
230  ! size mismatch
231  ier = 1
232  return
233  endif
234 
235  flat_index = indices(nd)-1
236  do id = nd-1, 1, -1
237  flat_index = flat_index*nsizes(id) + indices(id)-1
238  enddo
239  flat_index = flat_index + 1
240 
241  end subroutine cld_ntrp_get_flat_index
242 
243 !...............................................................................
244 !...............................................................................
245 pure subroutine cld_ntrp_get_cell_values(nsizes, fnodes, indices, fvals, ier)
246  integer, intent(in) :: nsizes(:) !< size of fnodes along each axis
247  real, intent(in) :: fnodes(:) !< flattened array of node values
248  integer, intent(in) :: indices(:) !< cell indices
249  real, intent(out) :: fvals(0:) !< returned array values in the cell
250  integer, intent(out) :: ier !< error flag (0=ok)
251 
252  integer id, nt, nd, flat_index, Ic, iflag
253  integer, dimension(size(nsizes)) :: cell_indices, node_indices
254  ier = 0
255  fvals = 0.
256 
257  nd = size(nsizes)
258  if(nd /= size(indices)) then
259  ! size mismatch
260  ier = 1
261  return
262  endif
263  if(2**nd > size(fvals)) then
264  ! not enough elements to hold result
265  ier = 2
266  return
267  endif
268  nt = 1
269  do id = 1, nd
270  nt = nt * nsizes(id)
271  enddo
272  if(nt /= size(fnodes)) then
273  ! not enough node values
274  ier = 3
275  return
276  endif
277 
278  do ic = 0, 2**nd-1
279  call cld_ntrp_expand_index(ic, cell_indices, iflag)
280  node_indices = indices + cell_indices
281  call cld_ntrp_get_flat_index(nsizes, node_indices, flat_index, iflag)
282  fvals(ic) = fnodes(flat_index)
283  enddo
284 
285  end subroutine cld_ntrp_get_cell_values
286 
287 #endif
288 end MODULE cloud_interpolator_mod
289 !===============================================================================
290 !> @}
291 ! close documentation grouping