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