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