FMS  2025.04
Flexible Modeling System
horiz_interp.inc
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 !> @addtogroup horiz_interp_mod
19 !> @{
20  !> @brief Creates a 1D @ref horiz_interp_type with the given parameters
21  subroutine horiz_interp_new_1d_ (Interp, lon_in, lat_in, lon_out, lat_out, verbose, &
22  interp_method, num_nbrs, max_dist, src_modulo, &
23  grid_at_center, mask_in, mask_out)
24 
25  !-----------------------------------------------------------------------
26  type(horiz_interp_type), intent(inout) :: Interp
27  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
28  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
29  integer, intent(in), optional :: verbose
30  character(len=*), intent(in), optional :: interp_method
31  integer, intent(in), optional :: num_nbrs
32  real(FMS_HI_KIND_), intent(in), optional :: max_dist
33  logical, intent(in), optional :: src_modulo
34  logical, intent(in), optional :: grid_at_center
35  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in !< dummy variable
36  real(FMS_HI_KIND_), intent(inout),dimension(:,:), optional :: mask_out !< dummy variable
37  !-----------------------------------------------------------------------
38  real(FMS_HI_KIND_), dimension(:,:), allocatable :: lon_src, lat_src, lon_dst, lat_dst
39  real(FMS_HI_KIND_), dimension(:), allocatable :: lon_src_1d, lat_src_1d, lon_dst_1d, lat_dst_1d
40  integer :: i, j, nlon_in, nlat_in, nlon_out, nlat_out
41  logical :: center
42  character(len=40) :: method
43  integer, parameter :: kindl = fms_hi_kind_ !> real kind size currently compiling
44  !-----------------------------------------------------------------------
45  call horiz_interp_init
46 
47  method = 'conservative'
48  if(present(interp_method)) method = interp_method
49 
50  select case (trim(method))
51  case ("conservative")
52  interp%interp_method = conserve
53  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose)
54  case ("bilinear")
55  interp%interp_method = bilinear
56  center = .false.
57  if(present(grid_at_center) ) center = grid_at_center
58  if(center) then
59  nlon_out = size(lon_out(:)); nlat_out = size(lat_out(:))
60  allocate(lon_dst(nlon_out,nlat_out), lat_dst(nlon_out,nlat_out))
61  do i = 1, nlon_out
62  lon_dst(i,:) = lon_out(i)
63  enddo
64  do j = 1, nlat_out
65  lat_dst(:,j) = lat_out(j)
66  enddo
67 
68  call horiz_interp_bilinear_new ( interp, lon_in, lat_in, lon_dst, lat_dst, &
69  verbose, src_modulo)
70  deallocate(lon_dst, lat_dst)
71  else
72  nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
73  nlon_out = size(lon_out(:))-1; nlat_out = size(lat_out(:))-1
74  allocate(lon_src_1d(nlon_in), lat_src_1d(nlat_in))
75  allocate(lon_dst(nlon_out,nlat_out), lat_dst(nlon_out,nlat_out))
76  do i = 1, nlon_in
77  lon_src_1d(i) = (lon_in(i) + lon_in(i+1)) * 0.5_kindl
78  enddo
79  do j = 1, nlat_in
80  lat_src_1d(j) = (lat_in(j) + lat_in(j+1)) * 0.5_kindl
81  enddo
82  do i = 1, nlon_out
83  lon_dst(i,:) = (lon_out(i) + lon_out(i+1)) * 0.5_kindl
84  enddo
85  do j = 1, nlat_out
86  lat_dst(:,j) = (lat_out(j) + lat_out(j+1)) * 0.5_kindl
87  enddo
88  call horiz_interp_bilinear_new ( interp, lon_src_1d, lat_src_1d, lon_dst, lat_dst, &
89  verbose, src_modulo)
90  deallocate(lon_src_1d, lat_src_1d, lon_dst, lat_dst)
91  endif
92  case ("bicubic")
93  interp%interp_method = bicubic
94  center = .false.
95  if(present(grid_at_center) ) center = grid_at_center
96  !No need to expand to 2d, horiz_interp_bicubic_new does 1d-1d
97  if(center) then
98  call horiz_interp_bicubic_new ( interp, lon_in, lat_in, lon_out, lat_out, &
99  verbose, src_modulo)
100  else
101  nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
102  nlon_out = size(lon_out(:))-1; nlat_out = size(lat_out(:))-1
103  allocate(lon_src_1d(nlon_in), lat_src_1d(nlat_in))
104  allocate(lon_dst_1d(nlon_out), lat_dst_1d(nlat_out))
105  do i = 1, nlon_in
106  lon_src_1d(i) = (lon_in(i) + lon_in(i+1)) * 0.5_kindl
107  enddo
108  do j = 1, nlat_in
109  lat_src_1d(j) = (lat_in(j) + lat_in(j+1)) * 0.5_kindl
110  enddo
111  do i = 1, nlon_out
112  lon_dst_1d(i) = (lon_out(i) + lon_out(i+1)) * 0.5_kindl
113  enddo
114  do j = 1, nlat_out
115  lat_dst_1d(j) = (lat_out(j) + lat_out(j+1)) * 0.5_kindl
116  enddo
117  call horiz_interp_bicubic_new ( interp, lon_src_1d, lat_src_1d, lon_dst_1d, lat_dst_1d, &
118  verbose, src_modulo)
119  deallocate(lon_src_1d, lat_src_1d, lon_dst_1d, lat_dst_1d)
120  endif
121  case ("spherical")
122  interp%interp_method = spherical
123  nlon_in = size(lon_in(:)); nlat_in = size(lat_in(:))
124  nlon_out = size(lon_out(:)); nlat_out = size(lat_out(:))
125  allocate(lon_src(nlon_in,nlat_in), lat_src(nlon_in,nlat_in))
126  allocate(lon_dst(nlon_out,nlat_out), lat_dst(nlon_out,nlat_out))
127  do i = 1, nlon_in
128  lon_src(i,:) = lon_in(i)
129  enddo
130  do j = 1, nlat_in
131  lat_src(:,j) = lat_in(j)
132  enddo
133  do i = 1, nlon_out
134  lon_dst(i,:) = lon_out(i)
135  enddo
136  do j = 1, nlat_out
137  lat_dst(:,j) = lat_out(j)
138  enddo
139  call horiz_interp_spherical_new ( interp, lon_src, lat_src, lon_dst, lat_dst, &
140  num_nbrs, max_dist, src_modulo)
141  deallocate(lon_src, lat_src, lon_dst, lat_dst)
142  case default
143  call mpp_error(fatal,'horiz_interp_mod: interp_method should be conservative, bilinear, bicubic, spherical')
144  end select
145 
146  !-----------------------------------------------------------------------
147  interp% HI_KIND_TYPE_ % is_allocated = .true.
148  interp%I_am_initialized = .true.
149 
150  end subroutine horiz_interp_new_1d_
151 
152 !#######################################################################
153 
154  subroutine horiz_interp_new_1d_src_ (Interp, lon_in, lat_in, lon_out, lat_out, &
155  verbose, interp_method, num_nbrs, max_dist, &
156  src_modulo, grid_at_center, mask_in, mask_out, is_latlon_out )
157 
158  type(horiz_interp_type), intent(inout) :: Interp
159  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
160  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
161  integer, intent(in), optional :: verbose
162  character(len=*), intent(in), optional :: interp_method
163  integer, intent(in), optional :: num_nbrs !< minimum number of neighbors
164  real(FMS_HI_KIND_), intent(in), optional :: max_dist
165  logical, intent(in), optional :: src_modulo
166  logical, intent(in), optional :: grid_at_center
167  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
168  real(FMS_HI_KIND_), intent(out),dimension(:,:), optional :: mask_out
169  logical, intent(in), optional :: is_latlon_out
170 
171  real(FMS_HI_KIND_), dimension(:,:), allocatable :: lon_src, lat_src
172  real(FMS_HI_KIND_), dimension(:), allocatable :: lon_src_1d, lat_src_1d
173  integer :: i, j, nlon_in, nlat_in
174  character(len=40) :: method
175  logical :: center
176  logical :: dst_is_latlon
177  integer, parameter :: kindl = fms_hi_kind_ !< real kind size currently compiling
178  !-----------------------------------------------------------------------
179  call horiz_interp_init
180 
181  method = 'conservative'
182  if(present(interp_method)) method = interp_method
183 
184  select case (trim(method))
185  case ("conservative")
186  interp%interp_method = conserve
187  !--- check to see if the source grid is regular lat-lon grid or not.
188  if(PRESENT(is_latlon_out)) then
189  dst_is_latlon = is_latlon_out
190  else
191  dst_is_latlon = is_lat_lon(lon_out, lat_out)
192  end if
193  if(dst_is_latlon ) then
194  if(present(mask_in)) then
195  if ( any(mask_in < -.0001_kindl) .or. any(mask_in > 1.0001_kindl)) &
196  call mpp_error(fatal, &
197  'horiz_interp_conserve_new_1d_src(horiz_interp_conserve_mod): input mask not between 0,1')
198  allocate(interp%HI_KIND_TYPE_%mask_in(size(mask_in,1), size(mask_in,2)) )
199  interp%HI_KIND_TYPE_%mask_in = mask_in
200  end if
201  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out(:,1), lat_out(1,:), &
202  verbose=verbose )
203  else
204  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out, lat_out, &
205  verbose=verbose, mask_in=mask_in, mask_out=mask_out )
206  end if
207  case ("bilinear")
208  interp%interp_method = bilinear
209  center = .false.
210  if(present(grid_at_center) ) center = grid_at_center
211  if(center) then
212  call horiz_interp_bilinear_new ( interp, lon_in, lat_in, lon_out, lat_out, &
213  verbose, src_modulo )
214  else
215  nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
216  allocate(lon_src_1d(nlon_in), lat_src_1d(nlat_in))
217  do i = 1, nlon_in
218  lon_src_1d(i) = (lon_in(i) + lon_in(i+1)) * 0.5_kindl
219  enddo
220  do j = 1, nlat_in
221  lat_src_1d(j) = (lat_in(j) + lat_in(j+1)) * 0.5_kindl
222  enddo
223  call horiz_interp_bilinear_new ( interp, lon_src_1d, lat_src_1d, lon_out, lat_out, &
224  verbose, src_modulo )
225  deallocate(lon_src_1d,lat_src_1d)
226  endif
227  case ("bicubic")
228  interp%interp_method = bicubic
229  center = .false.
230  if(present(grid_at_center) ) center = grid_at_center
231  if(center) then
232  call horiz_interp_bicubic_new ( interp, lon_in, lat_in, lon_out, lat_out, &
233  verbose, src_modulo )
234  else
235  nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
236  allocate(lon_src_1d(nlon_in), lat_src_1d(nlat_in))
237  do i = 1, nlon_in
238  lon_src_1d(i) = (lon_in(i) + lon_in(i+1)) * 0.5_kindl
239  enddo
240  do j = 1, nlat_in
241  lat_src_1d(j) = (lat_in(j) + lat_in(j+1)) * 0.5_kindl
242  enddo
243  call horiz_interp_bicubic_new ( interp, lon_src_1d, lat_src_1d, lon_out, lat_out, &
244  verbose, src_modulo )
245  deallocate(lon_src_1d,lat_src_1d)
246  endif
247  case ("spherical")
248  interp%interp_method = spherical
249  nlon_in = size(lon_in(:)); nlat_in = size(lat_in(:))
250  allocate(lon_src(nlon_in,nlat_in), lat_src(nlon_in,nlat_in))
251  do i = 1, nlon_in
252  lon_src(i,:) = lon_in(i)
253  enddo
254  do j = 1, nlat_in
255  lat_src(:,j) = lat_in(j)
256  enddo
257  call horiz_interp_spherical_new ( interp, lon_src, lat_src, lon_out, lat_out, &
258  num_nbrs, max_dist, src_modulo)
259  deallocate(lon_src, lat_src)
260  case default
261  call mpp_error(fatal,'interp_method should be conservative, bilinear, bicubic, spherical')
262  end select
263 
264  !-----------------------------------------------------------------------
265  interp% HI_KIND_TYPE_ % is_allocated = .true.
266  interp%I_am_initialized = .true.
267 
268  end subroutine horiz_interp_new_1d_src_
269 
270 !#######################################################################
271 
272  subroutine horiz_interp_new_2d_ (Interp, lon_in, lat_in, lon_out, lat_out, &
273  verbose, interp_method, num_nbrs, max_dist, &
274  src_modulo, mask_in, mask_out, is_latlon_in, is_latlon_out )
275  type(horiz_interp_type), intent(inout) :: Interp
276  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
277  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
278  integer, intent(in), optional :: verbose
279  character(len=*), intent(in), optional :: interp_method
280  integer, intent(in), optional :: num_nbrs
281  real(FMS_HI_KIND_), intent(in), optional :: max_dist
282  logical, intent(in), optional :: src_modulo
283  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
284  real(FMS_HI_KIND_), intent(out),dimension(:,:), optional :: mask_out
285  logical, intent(in), optional :: is_latlon_in, is_latlon_out
286  logical :: src_is_latlon, dst_is_latlon
287  character(len=40) :: method
288  integer, parameter :: kindl = fms_hi_kind_ !< real kind size currently compiling
289 !-----------------------------------------------------------------------
290  call horiz_interp_init
291 
292  method = 'bilinear'
293  if(present(interp_method)) method = interp_method
294 
295  select case (trim(method))
296  case ("conservative")
297  interp%interp_method = conserve
298  if(PRESENT(is_latlon_in)) then
299  src_is_latlon = is_latlon_in
300  else
301  src_is_latlon = is_lat_lon(lon_in, lat_in)
302  end if
303  if(PRESENT(is_latlon_out)) then
304  dst_is_latlon = is_latlon_out
305  else
306  dst_is_latlon = is_lat_lon(lon_out, lat_out)
307  end if
308  if(src_is_latlon .AND. dst_is_latlon) then
309  if(present(mask_in)) then
310  if ( any(mask_in < -0.0001_kindl) .or. any(mask_in > 1.0001_kindl)) then
311  call mpp_error(fatal, 'horiz_interp_conserve_new_2d(horiz_interp_conserve_mod):' // &
312  ' input mask not between 0,1')
313  endif
314  allocate(interp%HI_KIND_TYPE_%mask_in(size(mask_in,1), size(mask_in,2)) )
315  interp%HI_KIND_TYPE_%mask_in = mask_in
316  end if
317  call horiz_interp_conserve_new ( interp, lon_in(:,1), lat_in(1,:), lon_out(:,1), lat_out(1,:), &
318  verbose=verbose )
319  else if(src_is_latlon) then
320  call horiz_interp_conserve_new ( interp, lon_in(:,1), lat_in(1,:), lon_out, lat_out, &
321  verbose=verbose, mask_in=mask_in, mask_out=mask_out )
322  else if(dst_is_latlon) then
323  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out(:,1), lat_out(1,:), &
324  verbose=verbose, mask_in=mask_in, mask_out=mask_out )
325  else
326  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out, lat_out, &
327  verbose=verbose, mask_in=mask_in, mask_out=mask_out )
328  end if
329 
330  case ("spherical")
331  interp%interp_method = spherical
332  call horiz_interp_spherical_new ( interp, lon_in, lat_in, lon_out, lat_out, &
333  num_nbrs, max_dist, src_modulo )
334  case ("bilinear")
335  interp%interp_method = bilinear
336  call horiz_interp_bilinear_new ( interp, lon_in, lat_in, lon_out, lat_out, &
337  verbose, src_modulo )
338  case default
339  call mpp_error(fatal,'when source grid are 2d, interp_method should be spherical or bilinear')
340  end select
341 
342  interp% HI_KIND_TYPE_ % is_allocated = .true.
343  interp%I_am_initialized = .true.
344 
345  end subroutine horiz_interp_new_2d_
346 
347 !#######################################################################
348  subroutine horiz_interp_new_1d_dst_ (Interp, lon_in, lat_in, lon_out, lat_out, &
349  verbose, interp_method, num_nbrs, max_dist, src_modulo, mask_in, mask_out, is_latlon_in )
350  type(horiz_interp_type), intent(inout) :: Interp
351  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
352  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
353  integer, intent(in), optional :: verbose
354  character(len=*), intent(in), optional :: interp_method
355  integer, intent(in), optional :: num_nbrs
356  real(FMS_HI_KIND_), intent(in), optional :: max_dist
357  logical, intent(in), optional :: src_modulo
358  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
359  real(FMS_HI_KIND_), intent(out),dimension(:,:), optional :: mask_out
360  logical, intent(in), optional :: is_latlon_in
361 
362  character(len=40) :: method
363  integer, parameter :: kindl = fms_hi_kind_ !< real kind size currently compiling
364  !-------------some local variables-----------------------------------------------
365  integer :: i, j, nlon_out, nlat_out
366  real(FMS_HI_KIND_), dimension(:,:), allocatable :: lon_dst, lat_dst
367  logical :: src_is_latlon
368  !-----------------------------------------------------------------------
369  call horiz_interp_init
370 
371  method = 'bilinear'
372  if(present(interp_method)) method = interp_method
373 
374  nlon_out = size(lon_out(:)); nlat_out = size(lat_out(:))
375  allocate(lon_dst(nlon_out,nlat_out), lat_dst(nlon_out,nlat_out))
376  do i = 1, nlon_out
377  lon_dst(i,:) = lon_out(i)
378  enddo
379  do j = 1, nlat_out
380  lat_dst(:,j) = lat_out(j)
381  enddo
382 
383  select case (trim(method))
384  case ("conservative")
385  interp%interp_method = conserve
386  if(PRESENT(is_latlon_in)) then
387  src_is_latlon = is_latlon_in
388  else
389  src_is_latlon = is_lat_lon(lon_in, lat_in)
390  end if
391 
392  if(src_is_latlon) then
393  if(present(mask_in)) then
394  if ( any(mask_in < -0.0001_kindl) .or. any(mask_in > 1.0001_kindl)) &
395  call mpp_error(fatal, &
396  'horiz_interp_conserve_new_1d_dst(horiz_interp_conserve_mod): input mask not between 0,1')
397  allocate(interp%HI_KIND_TYPE_%mask_in(size(mask_in,1), size(mask_in,2)) )
398  interp%HI_KIND_TYPE_%mask_in = mask_in
399  end if
400  call horiz_interp_conserve_new ( interp, lon_in(:,1), lat_in(1,:), lon_out, lat_out, &
401  verbose=verbose)
402  else
403  call horiz_interp_conserve_new ( interp, lon_in, lat_in, lon_out, lat_out, &
404  verbose=verbose, mask_in=mask_in, mask_out=mask_out )
405  end if
406  case ("bilinear")
407  interp%interp_method = bilinear
408  call horiz_interp_bilinear_new ( interp, lon_in, lat_in, lon_dst, lat_dst, &
409  verbose, src_modulo )
410  case ("spherical")
411  interp%interp_method = spherical
412  call horiz_interp_spherical_new ( interp, lon_in, lat_in, lon_dst, lat_dst, &
413  num_nbrs, max_dist, src_modulo)
414  case default
415  call mpp_error(fatal,'when source grid are 2d, interp_method should be spherical or bilinear')
416  end select
417 
418  deallocate(lon_dst,lat_dst)
419 
420  interp% HI_KIND_TYPE_ % is_allocated = .true.
421  interp%I_am_initialized = .true.
422 
423  end subroutine horiz_interp_new_1d_dst_
424 
425 !#######################################################################
426 
427  subroutine horiz_interp_base_2d_ ( Interp, data_in, data_out, verbose, &
428  mask_in, mask_out, missing_value, missing_permit, &
429  err_msg, new_missing_handle )
430 !-----------------------------------------------------------------------
431  type (horiz_interp_type), intent(in) :: Interp
432  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
433  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
434  integer, intent(in), optional :: verbose
435  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
436  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
437  real(FMS_HI_KIND_), intent(in), optional :: missing_value
438  integer, intent(in), optional :: missing_permit
439  character(len=*), intent(out), optional :: err_msg
440  logical, intent(in), optional :: new_missing_handle
441 !-----------------------------------------------------------------------
442  if(present(err_msg)) err_msg = ''
443  if(.not.interp%I_am_initialized) then
444  if(fms_error_handler('horiz_interp','The horiz_interp_type variable is not initialized',err_msg)) return
445  endif
446 
447  select case(interp%interp_method)
448  case(conserve)
449  call horiz_interp_conserve(interp,data_in, data_out, verbose, mask_in, mask_out)
450  case(bilinear)
451  call horiz_interp_bilinear(interp,data_in, data_out, verbose, mask_in, mask_out, &
452  missing_value, missing_permit, new_missing_handle )
453  case(bicubic)
454  call horiz_interp_bicubic(interp,data_in, data_out, verbose, mask_in, mask_out, &
455  missing_value, missing_permit )
456  case(spherical)
457  call horiz_interp_spherical(interp,data_in, data_out, verbose, mask_in, mask_out, &
458  missing_value )
459  case default
460  call mpp_error(fatal,'interp_method should be conservative, bilinear, bicubic, spherical')
461  end select
462 
463  return
464 
465  end subroutine horiz_interp_base_2d_
466 
467 !#######################################################################
468 
469  !> Overload of interface HORIZ_INTERP_BASE_2D_
470  !! uses 3d arrays for data and mask
471  !! this allows for multiple interpolations with one call
472  subroutine horiz_interp_base_3d_ ( Interp, data_in, data_out, verbose, mask_in, mask_out, &
473  missing_value, missing_permit, err_msg )
474  !-----------------------------------------------------------------------
475  ! overload of interface HORIZ_INTERP_BASE_2D_
476  ! uses 3d arrays for data and mask
477  ! this allows for multiple interpolations with one call
478  !-----------------------------------------------------------------------
479  type (horiz_interp_type), intent(in) :: Interp
480  real(FMS_HI_KIND_), intent(in), dimension(:,:,:) :: data_in
481  real(FMS_HI_KIND_), intent(out), dimension(:,:,:) :: data_out
482  integer, intent(in), optional :: verbose
483  real(FMS_HI_KIND_), intent(in), dimension(:,:,:), optional :: mask_in
484  real(FMS_HI_KIND_), intent(out), dimension(:,:,:), optional :: mask_out
485  real(FMS_HI_KIND_), intent(in), optional :: missing_value
486  integer, intent(in), optional :: missing_permit
487  character(len=*), intent(out), optional :: err_msg
488  !-----------------------------------------------------------------------
489  integer :: n
490 
491  if(present(err_msg)) err_msg = ''
492  if(.not.interp%I_am_initialized) then
493  if(fms_error_handler('horiz_interp','The horiz_interp_type variable is not initialized',err_msg)) return
494  endif
495 
496  do n = 1, size(data_in,3)
497  if (present(mask_in))then
498  if(present(mask_out)) then
499  call horiz_interp( interp, data_in(:,:,n), data_out(:,:,n), &
500  verbose, mask_in(:,:,n), mask_out(:,:,n), &
501  missing_value, missing_permit )
502  else
503  call horiz_interp( interp, data_in(:,:,n), data_out(:,:,n), &
504  verbose, mask_in(:,:,n), missing_value = missing_value, &
505  missing_permit = missing_permit )
506  endif
507  else
508  if(present(mask_out)) then
509  call horiz_interp( interp, data_in(:,:,n), data_out(:,:,n), &
510  verbose, mask_out=mask_out(:,:,n), missing_value = missing_value, &
511  missing_permit = missing_permit )
512  else
513  call horiz_interp( interp, data_in(:,:,n), data_out(:,:,n), &
514  verbose, missing_value = missing_value, &
515  missing_permit = missing_permit )
516  endif
517  endif
518  enddo
519 
520  return
521 !-----------------------------------------------------------------------
522  end subroutine horiz_interp_base_3d_
523 
524 !#######################################################################
525 
526 !> Interpolates from a rectangular grid to rectangular grid.
527 !! interp_method can be the value conservative, bilinear or spherical.
528 !! horiz_interp_new don't need to be called before calling this routine.
529  subroutine horiz_interp_solo_1d_ ( data_in, lon_in, lat_in, lon_out, lat_out, &
530  data_out, verbose, mask_in, mask_out, &
531  interp_method, missing_value, missing_permit, &
532  num_nbrs, max_dist,src_modulo, grid_at_center )
533 !-----------------------------------------------------------------------
534  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
535  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
536  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
537  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
538  integer, intent(in), optional :: verbose
539  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
540  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
541  character(len=*), intent(in), optional :: interp_method
542  real(FMS_HI_KIND_), intent(in), optional :: missing_value
543  integer, intent(in), optional :: missing_permit
544  integer, intent(in), optional :: num_nbrs
545  real(FMS_HI_KIND_), intent(in), optional :: max_dist
546  logical, intent(in), optional :: src_modulo
547  logical, intent(in), optional :: grid_at_center
548 !-----------------------------------------------------------------------
549  type (horiz_interp_type) :: Interp
550 !-----------------------------------------------------------------------
551  call horiz_interp_init
552 
553  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
554  interp_method, num_nbrs, max_dist, src_modulo, grid_at_center )
555 
556  call horiz_interp ( interp, data_in, data_out, verbose, &
557  mask_in, mask_out, missing_value, missing_permit )
558 
559  call horiz_interp_del ( interp )
560 !-----------------------------------------------------------------------
561 
562  end subroutine horiz_interp_solo_1d_
563 
564 !#######################################################################
565 
566 !> Interpolates from a uniformly spaced grid to any output grid.
567 !! interp_method can be the value "onservative","bilinear" or "spherical".
568 !! horiz_interp_new don't need to be called before calling this routine.
569  subroutine horiz_interp_solo_1d_src_ ( data_in, lon_in, lat_in, lon_out, lat_out, &
570  data_out, verbose, mask_in, mask_out, &
571  interp_method, missing_value, missing_permit, &
572  num_nbrs, max_dist, src_modulo, grid_at_center )
573 !-----------------------------------------------------------------------
574  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
575  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
576  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
577  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
578  integer, intent(in), optional :: verbose
579  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
580  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
581  character(len=*), intent(in), optional :: interp_method
582  real(FMS_HI_KIND_), intent(in), optional :: missing_value
583  integer, intent(in), optional :: missing_permit
584  integer, intent(in), optional :: num_nbrs
585  real(FMS_HI_KIND_), intent(in), optional :: max_dist
586  logical, intent(in), optional :: src_modulo
587  logical, intent(in), optional :: grid_at_center
588 
589 !-----------------------------------------------------------------------
590  type (horiz_interp_type) :: Interp
591  logical :: dst_is_latlon
592  character(len=128) :: method
593 !-----------------------------------------------------------------------
594  call horiz_interp_init
595  method = 'conservative'
596  if(present(interp_method)) method = interp_method
597  dst_is_latlon = .true.
598  if(trim(method) == 'conservative') dst_is_latlon = is_lat_lon(lon_out, lat_out)
599 
600  if(dst_is_latlon) then
601  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
602  interp_method, num_nbrs, max_dist, src_modulo, &
603  grid_at_center, is_latlon_out = dst_is_latlon )
604  call horiz_interp ( interp, data_in, data_out, verbose, &
605  mask_in, mask_out, missing_value, missing_permit )
606  else
607  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
608  interp_method, num_nbrs, max_dist, src_modulo, &
609  grid_at_center, mask_in, mask_out, is_latlon_out = dst_is_latlon)
610 
611  call horiz_interp ( interp, data_in, data_out, verbose, &
612  missing_value=missing_value, missing_permit=missing_permit )
613  end if
614 
615  call horiz_interp_del ( interp )
616 
617 !-----------------------------------------------------------------------
618 
619  end subroutine horiz_interp_solo_1d_src_
620 
621 
622 !#######################################################################
623 
624 !> Interpolates from any grid to any grid. interp_method should be "spherical"
625 !! horiz_interp_new don't need to be called before calling this routine.
626  subroutine horiz_interp_solo_2d_ ( data_in, lon_in, lat_in, lon_out, lat_out, data_out, &
627  verbose, mask_in, mask_out, interp_method, missing_value,&
628  missing_permit, num_nbrs, max_dist, src_modulo )
629 !-----------------------------------------------------------------------
630  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
631  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
632  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
633  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
634  integer, intent(in), optional :: verbose
635  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
636  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
637  character(len=*), intent(in), optional :: interp_method
638  real(FMS_HI_KIND_), intent(in), optional :: missing_value
639  integer, intent(in), optional :: missing_permit
640  integer, intent(in), optional :: num_nbrs
641  real(FMS_HI_KIND_), intent(in), optional :: max_dist
642  logical, intent(in), optional :: src_modulo
643 !-----------------------------------------------------------------------
644  type (horiz_interp_type) :: Interp
645  logical :: dst_is_latlon, src_is_latlon
646  character(len=128) :: method
647 !-----------------------------------------------------------------------
648  call horiz_interp_init
649 
650  method = 'conservative'
651  if(present(interp_method)) method = interp_method
652  dst_is_latlon = .true.
653  src_is_latlon = .true.
654  if(trim(method) == 'conservative') then
655  dst_is_latlon = is_lat_lon(lon_out, lat_out)
656  src_is_latlon = is_lat_lon(lon_in, lat_in)
657  end if
658 
659  if(dst_is_latlon .and. src_is_latlon) then
660  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
661  interp_method, num_nbrs, max_dist, src_modulo, &
662  is_latlon_in=dst_is_latlon, is_latlon_out = dst_is_latlon )
663  call horiz_interp ( interp, data_in, data_out, verbose, &
664  mask_in, mask_out, missing_value, missing_permit )
665  else
666  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
667  interp_method, num_nbrs, max_dist, src_modulo, &
668  mask_in, mask_out, &
669  is_latlon_in=dst_is_latlon, is_latlon_out = dst_is_latlon)
670  call horiz_interp ( interp, data_in, data_out, verbose, &
671  missing_value=missing_value, missing_permit=missing_permit )
672  end if
673 
674  call horiz_interp_del ( interp )
675 
676 !-----------------------------------------------------------------------
677 
678  end subroutine horiz_interp_solo_2d_
679 
680 !#######################################################################
681 
682 !> interpolates from any grid to rectangular longitude/latitude grid.
683 !! interp_method should be "spherical".
684 !! horiz_interp_new don't need to be called before calling this routine.
685  subroutine horiz_interp_solo_1d_dst_ ( data_in, lon_in, lat_in, lon_out, lat_out, data_out, &
686  verbose, mask_in, mask_out,interp_method,missing_value, &
687  missing_permit, num_nbrs, max_dist, src_modulo)
688 !-----------------------------------------------------------------------
689  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
690  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
691  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
692  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
693  integer, intent(in), optional :: verbose
694  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
695  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
696  character(len=*), intent(in), optional :: interp_method
697  real(FMS_HI_KIND_), intent(in), optional :: missing_value
698  integer, intent(in), optional :: missing_permit
699  integer, intent(in), optional :: num_nbrs
700  real(FMS_HI_KIND_), intent(in), optional :: max_dist
701  logical, intent(in), optional :: src_modulo
702 !-----------------------------------------------------------------------
703  type (horiz_interp_type) :: Interp
704  logical :: src_is_latlon
705  character(len=128) :: method
706 !-----------------------------------------------------------------------
707  call horiz_interp_init
708 
709  method = 'conservative'
710  if(present(interp_method)) method = interp_method
711  src_is_latlon = .true.
712  if(trim(method) == 'conservative') src_is_latlon = is_lat_lon(lon_in, lat_in)
713 
714  if(src_is_latlon) then
715  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
716  interp_method, num_nbrs, max_dist, src_modulo, &
717  is_latlon_in = src_is_latlon )
718  call horiz_interp ( interp, data_in, data_out, verbose, &
719  mask_in, mask_out, missing_value, missing_permit )
720  else
721  call horiz_interp_new ( interp, lon_in, lat_in, lon_out, lat_out, verbose, &
722  interp_method, num_nbrs, max_dist, src_modulo, &
723  mask_in, mask_out, is_latlon_in = src_is_latlon)
724 
725  call horiz_interp ( interp, data_in, data_out, verbose, &
726  missing_value=missing_value, missing_permit=missing_permit )
727  end if
728 
729  call horiz_interp_del ( interp )
730 
731 !-----------------------------------------------------------------------
732 
733  end subroutine horiz_interp_solo_1d_dst_
734 
735 !#######################################################################
736 
737 !> Overloaded version of interface horiz_interp_solo_2
738  subroutine horiz_interp_solo_old_ (data_in, wb, sb, dx, dy, &
739  lon_out, lat_out, data_out, &
740  verbose, mask_in, mask_out)
741 
742 !-----------------------------------------------------------------------
743  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in !< Global input data stored from west to east
744  !! (1st dimension), south to north (2nd dimension)
745  real(FMS_HI_KIND_), intent(in) :: wb !< Longitude (radians) that correspond to western-most
746  !! boundary of grid box j=1 in array data_in
747  real(FMS_HI_KIND_), intent(in) :: sb !< Latitude (radians) that correspond to western-most
748  !! boundary of grid box j=1 in array data_in
749  real(FMS_HI_KIND_), intent(in) :: dx !< Grid spacing (in radians) for the longitude axis
750  !! (first dimension) for the input data
751  real(FMS_HI_KIND_), intent(in) :: dy !< Grid spacing (in radians) for the latitude axis
752  !! (first dimension) for the input data
753  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out !< The longitude edges (in radians) for output
754  !! data grid boxes. The values are for adjacent grid boxes
755  !! and must increase in value. If there are MLON grid boxes
756  !! there must be MLON+1 edge values
757  real(FMS_HI_KIND_), intent(in), dimension(:) :: lat_out !< The latitude edges (in radians) for output
758  !! data grid boxes. The values are for adjacent grid boxes
759  !! and may increase or decrease in value. If there are NLAT
760  !! grid boxes there must be NLAT+1 edge values
761  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out !< Output data on the output grid defined by grid box
762  integer, intent(in), optional :: verbose
763  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
764  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
765 !-----------------------------------------------------------------------
766  real(FMS_HI_KIND_), dimension(size(data_in,1)+1) :: blon_in
767  real(FMS_HI_KIND_), dimension(size(data_in,2)+1) :: blat_in
768  integer :: i, j, nlon_in, nlat_in
769  real(FMS_HI_KIND_) :: tpi
770  integer, parameter :: kindl = fms_hi_kind_ !< real size at compile time
771 !-----------------------------------------------------------------------
772  call horiz_interp_init
773 
774  tpi = 2.0_kindl * real(pi, fms_hi_kind_)
775  nlon_in = size(data_in,1)
776  nlat_in = size(data_in,2)
777 
778  do i = 1, nlon_in+1
779  blon_in(i) = wb + real(i-1, fms_hi_kind_)*dx
780  enddo
781  if (abs(blon_in(nlon_in+1)-blon_in(1)-tpi) < epsilon(blon_in)) &
782  blon_in(nlon_in+1)=blon_in(1)+tpi
783 
784  do j = 2, nlat_in
785  blat_in(j) = sb + real(j-1, fms_hi_kind_)*dy
786  enddo
787  blat_in(1) = -0.5_kindl * real(pi, fms_hi_kind_)
788  blat_in(nlat_in+1) = 0.5_kindl * real(pi, fms_hi_kind_)
789 
790 
791  call horiz_interp_solo_1d (data_in, blon_in, blat_in, &
792  lon_out, lat_out, data_out, &
793  verbose, mask_in, mask_out )
794 
795 !-----------------------------------------------------------------------
796 
797  end subroutine horiz_interp_solo_old_
798 
799 !#######################################################################
800 
801 
802  !####################################################################
803  function is_lat_lon_(lon, lat)
804  real(FMS_HI_KIND_), dimension(:,:), intent(in) :: lon, lat
805  logical :: IS_LAT_LON_
806  integer :: i, j, nlon, nlat, num
807 
808  is_lat_lon_ = .true.
809  nlon = size(lon,1)
810  nlat = size(lon,2)
811  loop_lat: do j = 1, nlat
812  do i = 2, nlon
813  if(lat(i,j) .NE. lat(1,j)) then
814  is_lat_lon_ = .false.
815  exit loop_lat
816  end if
817  end do
818  end do loop_lat
819 
820  if(is_lat_lon_) then
821  loop_lon: do i = 1, nlon
822  do j = 2, nlat
823  if(lon(i,j) .NE. lon(i,1)) then
824  is_lat_lon_ = .false.
825  exit loop_lon
826  end if
827  end do
828  end do loop_lon
829  end if
830 
831  num = 0
832  if(is_lat_lon_) num = 1
833  call mpp_min(num)
834  if(num == 1) then
835  is_lat_lon_ = .true.
836  else
837  is_lat_lon_ = .false.
838  end if
839 
840  return
841  end function is_lat_lon_
842 
843  !> Subroutine for reading a weight file and use it to fill in the horiz interp type
844 !! for the bilinear interpolation method.
845  subroutine horiz_interp_read_weights_(Interp, weight_filename, lon_out, lat_out, lon_in, lat_in, &
846  weight_file_source, interp_method, isw, iew, jsw, jew, nglon, nglat)
847  type(horiz_interp_type), intent(inout) :: Interp !< Horiz interp time to fill
848  character(len=*), intent(in) :: weight_filename !< Name of the weight file
849  real(FMS_HI_KIND_), intent(in) :: lat_out(:,:) !< Output (model) latitude
850  real(FMS_HI_KIND_), intent(in) :: lon_out(:,:) !< Output (model) longitude
851  real(FMS_HI_KIND_), intent(in) :: lat_in(:) !< Input (data) latitude
852  real(FMS_HI_KIND_), intent(in) :: lon_in(:) !< Input (data) longitude
853  character(len=*), intent(in) :: weight_file_source !< Source of the weight file
854  character(len=*), intent(in) :: interp_method !< The interp method to use
855  integer, intent(in) :: isw, iew, jsw, jew !< Starting and ending indices of the compute domain
856  integer, intent(in) :: nglon !< Number of longitudes in the global domain
857  integer, intent(in) :: nglat !< Number of latitudes in the globl domain
858 
859  integer :: i, j !< For do loops
860  integer :: nlon_in !< Number of longitude in the data
861  integer :: nlat_in !< Number of latitude in the data grid
862  real(FMS_HI_KIND_), allocatable :: lon_src_1d(:) !< Center points of the longitude data grid
863  real(FMS_HI_KIND_), allocatable :: lat_src_1d(:) !< Center points of the lattiude data grid
864  integer, parameter :: kindl = fms_hi_kind_ !< real kind size currently compiling
865 
866  select case (trim(interp_method))
867  case ("bilinear")
868  !! This is to reproduce the behavior in horiz_interp_new
869  !! The subroutine assumes that the data grid (lon_in, lat_in) are
870  !! the edges and not the centers.
871  !! Data_override passes in the edges, which are calculated using the axis_edges subroutine
872  nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
873  allocate(lon_src_1d(nlon_in), lat_src_1d(nlat_in))
874  do i = 1, nlon_in
875  lon_src_1d(i) = (lon_in(i) + lon_in(i+1)) * 0.5_kindl
876  enddo
877  do j = 1, nlat_in
878  lat_src_1d(j) = (lat_in(j) + lat_in(j+1)) * 0.5_kindl
879  enddo
880 
881  call horiz_interp_read_weights_bilinear(interp, weight_filename, lon_out, lat_out, &
882  lon_src_1d, lat_src_1d, weight_file_source, interp_method, &
883  isw, iew, jsw, jew, nglon, nglat)
884  deallocate(lon_src_1d,lat_src_1d)
885  case default
886  call mpp_error(fatal, "Reading weight from file is not supported for the "//&
887  trim(interp_method)//" method. It is currently only supported for bilinear")
888  end select
889  end subroutine horiz_interp_read_weights_
890 !> @}
horiz_interp_solo_1d_dst_
subroutine horiz_interp_solo_1d_dst_(data_in, lon_in, lat_in, lon_out, lat_out, data_out, verbose, mask_in, mask_out, interp_method, missing_value, missing_permit, num_nbrs, max_dist, src_modulo)
interpolates from any grid to rectangular longitude/latitude grid. interp_method should be "spherical...
Definition: horiz_interp.inc:688
horiz_interp_new_1d_src_
subroutine horiz_interp_new_1d_src_(Interp, lon_in, lat_in, lon_out, lat_out, verbose, interp_method, num_nbrs, max_dist, src_modulo, grid_at_center, mask_in, mask_out, is_latlon_out)
Definition: horiz_interp.inc:157
horiz_interp_solo_1d_
subroutine horiz_interp_solo_1d_(data_in, lon_in, lat_in, lon_out, lat_out, data_out, verbose, mask_in, mask_out, interp_method, missing_value, missing_permit, num_nbrs, max_dist, src_modulo, grid_at_center)
Interpolates from a rectangular grid to rectangular grid. interp_method can be the value conservative...
Definition: horiz_interp.inc:533
horiz_interp_base_3d_
subroutine horiz_interp_base_3d_(Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, missing_permit, err_msg)
Overload of interface HORIZ_INTERP_BASE_2D_ uses 3d arrays for data and mask this allows for multiple...
Definition: horiz_interp.inc:474
horiz_interp_read_weights_
subroutine horiz_interp_read_weights_(Interp, weight_filename, lon_out, lat_out, lon_in, lat_in, weight_file_source, interp_method, isw, iew, jsw, jew, nglon, nglat)
Subroutine for reading a weight file and use it to fill in the horiz interp type for the bilinear int...
Definition: horiz_interp.inc:847
horiz_interp_new_1d_
subroutine horiz_interp_new_1d_(Interp, lon_in, lat_in, lon_out, lat_out, verbose, interp_method, num_nbrs, max_dist, src_modulo, grid_at_center, mask_in, mask_out)
Definition: horiz_interp.inc:24
horiz_interp_solo_2d_
subroutine horiz_interp_solo_2d_(data_in, lon_in, lat_in, lon_out, lat_out, data_out, verbose, mask_in, mask_out, interp_method, missing_value, missing_permit, num_nbrs, max_dist, src_modulo)
Interpolates from any grid to any grid. interp_method should be "spherical" horiz_interp_new don't ne...
Definition: horiz_interp.inc:629
horiz_interp_solo_1d_src_
subroutine horiz_interp_solo_1d_src_(data_in, lon_in, lat_in, lon_out, lat_out, data_out, verbose, mask_in, mask_out, interp_method, missing_value, missing_permit, num_nbrs, max_dist, src_modulo, grid_at_center)
Interpolates from a uniformly spaced grid to any output grid. interp_method can be the value "onserva...
Definition: horiz_interp.inc:573
horiz_interp_solo_old_
subroutine horiz_interp_solo_old_(data_in, wb, sb, dx, dy, lon_out, lat_out, data_out, verbose, mask_in, mask_out)
Overloaded version of interface horiz_interp_solo_2.
Definition: horiz_interp.inc:741