MOM_interpolate.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> This module provides added functionality to the FMS temporal and spatial interpolation routines
7
8use mom_array_transform, only : allocate_rotated_array, rotate_array
9use mom_error_handler, only : mom_error, fatal
10use mom_interp_infra, only : time_interp_extern, init_external_field=>init_extern_field
11use mom_interp_infra, only : time_interp_external_init=>time_interp_extern_init
12use mom_interp_infra, only : horiz_interp_type
13use mom_interp_infra, only : get_external_field_info_infra => get_external_field_info
14use mom_interp_infra, only : run_horiz_interp, build_horiz_interp_weights
15use mom_interp_infra, only : external_field
16use mom_io_infra, only : axistype
17use mom_io_infra, only : get_axis_size, get_axis_data
19use mom_time_manager, only : time_type, set_date, operator(+), operator(<), operator(>)
20
21implicit none ; private
22
23!> Data type used to store information about forcing datasets that are time series
24!! E.g. how do we align the data in the model with the time axis in the file?
25type, public :: forcing_timeseries_dataset
26 character(len=200) :: file_name !< name of file containing river flux forcing
27 logical :: l_time_varying !< .true. => forcing is dependent on model time, .false. => static forcing
28 ! logical :: l_FMS_modulo !< .true. => let FMS handle determining time level to read (e.g. for climatologies)
29 type(time_type) :: data_forcing !< convert data_forcing_year to time type
30 type(time_type) :: data_start !< convert data_start_year to time type
31 type(time_type) :: data_end !< convert data_end_year to time type
32 type(time_type) :: m2d_offset !< add to model time to get data time
34
35public :: time_interp_external, init_external_field, time_interp_external_init
37public :: horiz_interp_type, run_horiz_interp, build_horiz_interp_weights
38public :: external_field
41
42!> Read a field based on model time, and rotate to the model domain.
43interface time_interp_external
44 module procedure time_interp_external_0d
45 module procedure time_interp_external_2d
46 module procedure time_interp_external_3d
47end interface time_interp_external
48
49contains
50
51!> Read a scalar field based on model time.
52subroutine time_interp_external_0d(field, time, data_in, verbose, scale)
53 type(external_field), intent(in) :: field !< Handle for time interpolated field
54 type(time_type), intent(in) :: time !< The target time for the data
55 real, intent(inout) :: data_in !< The interpolated value in arbitrary units [A ~> a]
56 logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging
57 real, optional, intent(in) :: scale !< A scaling factor that new values of data_in are
58 !! multiplied by before it is returned [A a-1 ~> 1]
59 real :: data_in_pre_scale ! The input data before rescaling [a]
60 real :: I_scale ! The inverse of scale [a A-1 ~> 1]
61
62 ! Store the input value in case the scaling factor is perfectly invertable.
63 data_in_pre_scale = data_in
64 i_scale = 1.0
65 if (present(scale)) then ; if ((scale /= 1.0) .and. (scale /= 0.0)) then
66 ! Because time_interp_extern has the ability to only set some values, but no clear
67 ! mechanism to determine which values have been set, the input data has to
68 ! be unscaled so that it will have the right values when it is returned.
69 i_scale = 1.0 / scale
70 data_in = data_in * i_scale
71 endif ; endif
72
73 call time_interp_extern(field, time, data_in, verbose=verbose)
74
75 if (present(scale)) then ; if (scale /= 1.0) then
76 ! Rescale data that has been newly set and restore the scaling of unset data.
77 if (data_in == i_scale * data_in_pre_scale) then
78 data_in = data_in_pre_scale
79 else
80 data_in = scale * data_in
81 endif
82 endif ; endif
83
84end subroutine time_interp_external_0d
85
86!> Read a 2d field from an external based on model time, potentially including horizontal
87!! interpolation and rotation of the data
88subroutine time_interp_external_2d(field, time, data_in, interp, &
89 verbose, horz_interp, mask_out, turns, scale)
90 type(external_field), intent(in) :: field !< Handle for time interpolated field
91 type(time_type), intent(in) :: time !< The target time for the data
92 real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated
93 !! values in arbitrary units [A ~> a]
94 integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method
95 logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging
96 type(horiz_interp_type), &
97 optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation
98 logical, dimension(:,:), &
99 optional, intent(out) :: mask_out !< An array that is true where there is valid data
100 integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data
101 real, optional, intent(in) :: scale !< A scaling factor that new values of data_in are
102 !! multiplied by before it is returned [A a-1 ~> 1]
103
104 real, allocatable :: data_in_pre_scale(:,:) ! The input data before rescaling [a]
105 real, allocatable :: data_pre_rot(:,:) ! The unscaled input data before rotation [a]
106 real :: I_scale ! The inverse of scale [a A-1 ~> 1]
107 integer :: qturns ! The number of quarter turns to rotate the data
108 integer :: i, j
109
110 ! TODO: Mask rotation requires logical array rotation support
111 if (present(mask_out)) &
112 call mom_error(fatal, "Rotation of masked output not yet support")
113
114 if (present(scale)) then ; if ((scale /= 1.0) .and. (scale /= 0.0)) then
115 ! Because time_interp_extern has the ability to only set some values, but no clear mechanism
116 ! to determine which values have been set, the input data has to be unscaled so that it will
117 ! have the right values when it is returned. It may be a problem for some compiler settings
118 ! if there are NaNs in data_in, but they will not spread.
119 if (abs(fraction(scale)) /= 1.0) then
120 ! This scaling factor may not be perfectly invertable, so store the input value
121 allocate(data_in_pre_scale, source=data_in)
122 endif
123 i_scale = 1.0 / scale
124 data_in(:,:) = i_scale * data_in(:,:)
125 endif ; endif
126
127 qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
128
129 if (qturns == 0) then
130 call time_interp_extern(field, time, data_in, interp=interp, &
131 verbose=verbose, horz_interp=horz_interp)
132 else
133 call allocate_rotated_array(data_in, [1,1], -qturns, data_pre_rot)
134 call time_interp_extern(field, time, data_pre_rot, interp=interp, &
135 verbose=verbose, horz_interp=horz_interp)
136 call rotate_array(data_pre_rot, turns, data_in)
137 deallocate(data_pre_rot)
138 endif
139
140 if (present(scale)) then ; if (scale /= 1.0) then
141 ! Rescale data that has been newly set and restore the scaling of unset data.
142 if ((abs(fraction(scale)) /= 1.0) .and. (scale /= 0.0)) then
143 do j=lbound(data_in,2),ubound(data_in,2) ; do i=lbound(data_in,1),ubound(data_in,1)
144 ! This handles the case where scale is not exactly invertable for data
145 ! values that have not been modified by time_interp_extern.
146 if (data_in(i,j) == i_scale * data_in_pre_scale(i,j)) then
147 data_in(i,j) = data_in_pre_scale(i,j)
148 else
149 data_in(i,j) = scale * data_in(i,j)
150 endif
151 enddo ; enddo
152 else
153 data_in(:,:) = scale * data_in(:,:)
154 endif
155 endif ; endif
156
157end subroutine time_interp_external_2d
158
159
160!> Read a 3d field based on model time, and rotate to the model grid
161subroutine time_interp_external_3d(field, time, data_in, interp, &
162 verbose, horz_interp, mask_out, turns, scale)
163 type(external_field), intent(in) :: field !< Handle for time interpolated field
164 type(time_type), intent(in) :: time !< The target time for the data
165 real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated
166 !! values in arbitrary units [A ~> a]
167 integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method
168 logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging
169 type(horiz_interp_type), &
170 optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation
171 logical, dimension(:,:,:), &
172 optional, intent(out) :: mask_out !< An array that is true where there is valid data
173 integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data
174 real, optional, intent(in) :: scale !< A scaling factor that new values of data_in are
175 !! multiplied by before it is returned [A a-1 ~> 1]
176
177 real, allocatable :: data_in_pre_scale(:,:,:) ! The input data before rescaling [a]
178 real, allocatable :: data_pre_rot(:,:,:) ! The unscaled input data before rotation [a]
179 real :: I_scale ! The inverse of scale [a A-1 ~> 1]
180 integer :: qturns ! The number of quarter turns to rotate the data
181 integer :: i, j, k
182
183 ! TODO: Mask rotation requires logical array rotation support
184 if (present(mask_out)) &
185 call mom_error(fatal, "Rotation of masked output not yet support")
186
187 if (present(scale)) then ; if ((scale /= 1.0) .and. (scale /= 0.0)) then
188 ! Because time_interp_extern has the ability to only set some values, but no clear mechanism
189 ! to determine which values have been set, the input data has to be unscaled so that it will
190 ! have the right values when it is returned. It may be a problem for some compiler settings
191 ! if there are NaNs in data_in, but they will not spread.
192 if (abs(fraction(scale)) /= 1.0) then
193 ! This scaling factor may not be perfectly invertable, so store the input value
194 allocate(data_in_pre_scale, source=data_in)
195 endif
196 i_scale = 1.0 / scale
197 data_in(:,:,:) = i_scale * data_in(:,:,:)
198 endif ; endif
199
200 qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
201
202 if (qturns == 0) then
203 call time_interp_extern(field, time, data_in, interp=interp, &
204 verbose=verbose, horz_interp=horz_interp)
205 else
206 call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre_rot)
207 call time_interp_extern(field, time, data_pre_rot, interp=interp, &
208 verbose=verbose, horz_interp=horz_interp)
209 call rotate_array(data_pre_rot, turns, data_in)
210 deallocate(data_pre_rot)
211 endif
212
213 if (present(scale)) then ; if (scale /= 1.0) then
214 ! Rescale data that has been newly set and restore the scaling of unset data.
215 if ((abs(fraction(scale)) /= 1.0) .and. (scale /= 0.0)) then
216 do k=lbound(data_in,3),ubound(data_in,3)
217 do j=lbound(data_in,2),ubound(data_in,2)
218 do i=lbound(data_in,1),ubound(data_in,1)
219 ! This handles the case where scale is not exactly invertable for data
220 ! values that have not been modified by time_interp_extern.
221 if (data_in(i,j,k) == i_scale * data_in_pre_scale(i,j,k)) then
222 data_in(i,j,k) = data_in_pre_scale(i,j,k)
223 else
224 data_in(i,j,k) = scale * data_in(i,j,k)
225 endif
226 enddo
227 enddo
228 enddo
229 else
230 data_in(:,:,:) = scale * data_in(:,:,:)
231 endif
232 endif ; endif
233
234end subroutine time_interp_external_3d
235
236!> Set time_type variables in forcing_timeseries_dataset type based on integer input
237!! TODO: make this part of forcing_timeseries_dataset class if OO is okay in MOM6?
238subroutine forcing_timeseries_set_time_type_vars(data_start_year, data_end_year, data_ref_year, &
239 model_ref_year, data_forcing_year, forcing_dataset)
240
241 integer, intent(in) :: data_start_year !< first year of data to read
242 !! (this is ignored for static forcing)
243 integer, intent(in) :: data_end_year !< last year of data to read
244 !! (this is ignored for static forcing)
245 integer, intent(in) :: data_ref_year !< for time-varying forcing, align
246 !! data_ref_year in file with
247 !! model_ref_year in model
248 integer, intent(in) :: model_ref_year !< for time-varying forcing, align
249 !! data_ref_year in file with
250 !! model_ref_year in model
251 integer, intent(in) :: data_forcing_year !< for static forcing, read file at this
252 !! date (this is ignored for time-varying
253 !! forcing)
254 type(forcing_timeseries_dataset), intent(inout) :: forcing_dataset !< information about forcing file
255
256 if (forcing_dataset%l_time_varying) then
257 forcing_dataset%data_start = set_date(data_start_year, 1, 1)
258 forcing_dataset%data_end = set_date(data_end_year, 1, 1)
259 forcing_dataset%m2d_offset = set_date(data_ref_year - model_ref_year, 1, 1)
260 else
261 forcing_dataset%data_forcing = set_date(data_forcing_year, 1, 1)
262 endif
263
265
266!> If necessary, apply an offset to convert from model time to forcing time and then
267!! ensure result is within acceptable bounds
268function map_model_time_to_forcing_time(Time, forcing_dataset)
269
270 type(time_type), intent(in) :: time !< Model time
271 type(forcing_timeseries_dataset), intent(in) :: forcing_dataset !< information about forcing file
272 type(time_type) :: map_model_time_to_forcing_time !< time to read forcing file
273
274 if (forcing_dataset%l_time_varying) then
275 map_model_time_to_forcing_time = time + forcing_dataset%m2d_offset
276 ! If Time + offset is not between data_start and data_end, use whichever of those values is closer
277 if (map_model_time_to_forcing_time < forcing_dataset%data_start) &
278 map_model_time_to_forcing_time = forcing_dataset%data_start
279 if (map_model_time_to_forcing_time > forcing_dataset%data_end) &
280 map_model_time_to_forcing_time = forcing_dataset%data_end
281 else
282 map_model_time_to_forcing_time = forcing_dataset%data_forcing
283 endif
284
286
287
288subroutine get_external_field_info(field, size, axes, missing)
289 type(external_field), intent(in) :: field
290 !< Handle for time interpolated external field returned from a previous
291 !! call to init_external_field()
292 integer, optional, intent(inout) :: size(4)
293 !< Dimension sizes for the input data
294 type(axis_info), optional, intent(inout) :: axes(4)
295 !< Axis types for the input data
296 real, optional, intent(inout) :: missing
297 !< Missing value for the input data
298
299 type(axistype) :: axes_infra(4)
300 ! Axis as represented in the infra
301 character(len=256) :: axis_name
302 ! Axis name
303 real, allocatable :: ax_data(:)
304 ! Axis points
305
306 integer :: n
307 ! Axis index
308 integer :: ax_size
309 ! Axis size
310
311 if (present(axes)) then
312 call get_external_field_info_infra(field, size=size, axes=axes_infra, &
313 missing=missing)
314 ! TODO: Most of these methods were written to expect four dimensions.
315 do n=1,4
316 ! Convert axistype to axis_info
317 ax_size = get_axis_size(axes_infra(n))
318 allocate(ax_data(ax_size))
319 call get_axis_data(axes_infra(n), axis_name, ax_data)
320 call set_axis_info(axes(n), trim(axis_name), ax_data=ax_data)
321 deallocate(ax_data)
322 enddo
323 else
324 call get_external_field_info_infra(field, size=size, missing=missing)
325 endif
326end subroutine get_external_field_info
327
328
329end module mom_interpolate