interface for MOM_interp_infra (#1310)

This PR provides new explicit interfaces to the horizontal interpolation routines with thin-layer wrappers of the underlying FMS infrastructure routines. It also renames the key routines, from horiz_interp_new to build_horiz_interp_weights and from horiz_interp to run_horiz_interp, to more accurately reflect what these routines do.  Successful testing has verified that all answers are bitwise identical, and there are no changes to any output files.

src/framework/MOM_horizontal_regridding.F90

@@ -12,7 +12,7 @@ module MOM_horizontal_regridding
 use MOM_file_parser,   only : get_param, log_param, log_version, param_file_type
 use MOM_grid,          only : ocean_grid_type
 use MOM_interpolate,   only : time_interp_external, get_external_field_info, horiz_interp_init
-use MOM_interpolate,   only : horiz_interp_new, horiz_interp, horiz_interp_type
+use MOM_interpolate,   only : build_horiz_interp_weights, run_horiz_interp, horiz_interp_type
 use MOM_io_infra,      only : axistype, get_axis_data
 use MOM_time_manager,  only : time_type
 
@@ -526,8 +526,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
 !   call fms routine horiz_interp to interpolate input level data to model horizontal grid
     if (.not. is_ongrid) then
       if (k == 1) then
-        call horiz_interp_new(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
-              interp_method='bilinear', src_modulo=.true.)
+        call build_horiz_interp_weights(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
+                                        interp_method='bilinear', src_modulo=.true.)
       endif
 
       if (debug) then
@@ -539,7 +539,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
     if (is_ongrid) then
       tr_out(is:ie,js:je)=tr_in(is:ie,js:je)
     else
-      call horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.)
+      call run_horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), &
+                            missing_value=missing_value, new_missing_handle=.true.)
     endif
 
     mask_out=1.0
@@ -810,8 +811,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
       ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
       if (k == 1) then
-        call horiz_interp_new(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
-                              interp_method='bilinear', src_modulo=.true.)
+        call build_horiz_interp_weights(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
+                                        interp_method='bilinear', src_modulo=.true.)
       endif
 
       if (debug) then
@@ -820,8 +821,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
       tr_out(:,:) = 0.0
 
-      call horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, &
-                        new_missing_handle=.true.)
+      call run_horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, &
+                            new_missing_handle=.true.)
 
       mask_out(:,:) = 1.0
       do j=js,je ; do i=is,ie

src/framework/MOM_interp_infra.F90

@@ -14,9 +14,10 @@ module MOM_interp_infra
 
 implicit none ; private
 
+public :: horiz_interp_type, horiz_interp_init
 public :: time_interp_extern, init_extern_field, time_interp_external_init
 public :: get_external_field_info
-public :: horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new
+public :: run_horiz_interp, build_horiz_interp_weights
 
 !> Read a field based on model time, and rotate to the model domain.
 interface time_interp_extern
@@ -25,8 +26,131 @@ module MOM_interp_infra
   module procedure time_interp_extern_3d
 end interface time_interp_extern
 
+!> perform horizontal interpolation of field
+interface run_horiz_interp
+  module procedure horiz_interp_from_weights_field2d
+  module procedure horiz_interp_from_weights_field3d
+end interface
+
+!> build weights for horizontal interpolation of field
+interface build_horiz_interp_weights
+  module procedure build_horiz_interp_weights_2d_to_2d
+end interface build_horiz_interp_weights
+
 contains
 
+!> perform horizontal interpolation of a 2d field using pre-computed weights
+!! source and destination coordinates are 2d
+subroutine horiz_interp_from_weights_field2d(Interp, data_in, data_out, verbose, &
+                                             mask_in, mask_out, missing_value, missing_permit, &
+                                             err_msg, new_missing_handle)
+
+  type(horiz_interp_type), intent(in) :: Interp           !< type containing interpolation
+                                                          !! options/weights
+  real, intent(in), dimension(:,:) :: data_in             !< input data
+  real, intent(out), dimension(:,:) :: data_out           !< output data
+  integer, intent(in), optional :: verbose                !< verbosity level
+  real, intent(in), dimension(:,:), optional :: mask_in   !< mask for input data
+  real, intent(out), dimension(:,:), optional :: mask_out !< mask for output data
+  real, intent(in), optional :: missing_value             !< missing value
+  integer, intent(in), optional :: missing_permit         !< number of allowed points with missing value
+                                                          !! for interpolation (0-3)
+  character(len=*), intent(out), optional :: err_msg      !< error message
+  logical, intent(in), optional :: new_missing_handle     !< unknown
+
+  call horiz_interp(Interp, data_in, data_out, verbose, &
+                    mask_in, mask_out, missing_value, missing_permit, &
+                    err_msg, new_missing_handle )
+
+end subroutine horiz_interp_from_weights_field2d
+
+
+!> perform horizontal interpolation of a 3d field using pre-computed weights
+!! source and destination coordinates are 2d
+subroutine horiz_interp_from_weights_field3d(Interp, data_in, data_out, verbose, mask_in, mask_out, &
+                                             missing_value, missing_permit, err_msg)
+
+  type(horiz_interp_type), intent(in) :: Interp             !< type containing interpolation
+                                                            !! options/weights
+  real, intent(in), dimension(:,:,:) :: data_in             !< input data
+  real, intent(out), dimension(:,:,:) :: data_out           !< output data
+  integer, intent(in), optional :: verbose                  !< verbosity level
+  real, intent(in), dimension(:,:,:), optional :: mask_in   !< mask for input data
+  real, intent(out), dimension(:,:,:), optional :: mask_out !< mask for output data
+  real, intent(in), optional :: missing_value               !< missing value
+  integer, intent(in), optional :: missing_permit           !< number of allowed points with missing value
+                                                            !! for interpolation (0-3)
+  character(len=*), intent(out), optional :: err_msg        !< error message
+
+  call horiz_interp(Interp, data_in, data_out, verbose, mask_in, mask_out, &
+                    missing_value, missing_permit, err_msg)
+
+end subroutine horiz_interp_from_weights_field3d
+
+
+!> build horizontal interpolation weights from source grid defined by 2d lon/lat to destination grid
+!! defined by 2d lon/lat
+subroutine build_horiz_interp_weights_2d_to_2d(Interp, lon_in, lat_in, lon_out, lat_out, &
+                                               verbose, interp_method, num_nbrs, max_dist, &
+                                               src_modulo, mask_in, mask_out, &
+                                               is_latlon_in, is_latlon_out)
+
+  type(horiz_interp_type), intent(inout) :: Interp         !< type containing interpolation options/weights
+  real, intent(in), dimension(:,:) :: lon_in               !< input longitude 2d
+  real, intent(in), dimension(:,:) :: lat_in               !< input latitude 2d
+  real, intent(in), dimension(:,:) :: lon_out              !< output longitude 2d
+  real, intent(in), dimension(:,:) :: lat_out              !< output latitude 2d
+  integer, intent(in), optional :: verbose                 !< verbosity level
+  character(len=*), intent(in), optional :: interp_method  !< interpolation method
+  integer, intent(in), optional :: num_nbrs                !< number of nearest neighbors
+  real, intent(in), optional :: max_dist                   !< maximum region of influence
+  logical, intent(in), optional :: src_modulo              !< periodicity of E-W boundary
+  real, intent(in), dimension(:,:), optional :: mask_in    !< mask for input data
+  real, intent(inout),dimension(:,:), optional :: mask_out !< mask for output data
+  logical, intent(in), optional :: is_latlon_in            !< input grid is regular lat/lon grid
+  logical, intent(in), optional :: is_latlon_out           !< output grid is regular lat/lon grid
+
+  call horiz_interp_new(Interp, lon_in, lat_in, lon_out, lat_out, &
+                        verbose, interp_method, num_nbrs, max_dist, &
+                        src_modulo, mask_in, mask_out, &
+                        is_latlon_in, is_latlon_out)
+
+end subroutine build_horiz_interp_weights_2d_to_2d
+
+
+!> get size of an external field from field index
+function get_extern_field_size(index)
+
+  integer :: index                     !< field index
+  integer :: get_extern_field_size(4)  !< field size
+
+  get_extern_field_size = get_external_field_size(index)
+
+end function get_extern_field_size
+
+
+!> get axes of an external field from field index
+function get_extern_field_axes(index)
+
+  integer :: index                                      !< field index
+  type(axistype), dimension(4) :: get_extern_field_axes !< field axes
+
+  get_extern_field_axes = get_external_field_axes(index)
+
+end function get_extern_field_axes
+
+
+!> get missing value of an external field from field index
+function get_extern_field_missing(index)
+
+  integer :: index                 !< field index
+  real :: get_extern_field_missing !< field missing value
+
+  get_extern_field_missing = get_external_field_missing(index)
+
+end function get_extern_field_missing
+
+
 !> Get information about the external fields.
 subroutine get_external_field_info(field_id, size, axes, missing)
   integer,                                intent(in)    :: field_id !< The integer index of the external
@@ -37,15 +161,15 @@ subroutine get_external_field_info(field_id, size, axes, missing)
   real,                         optional, intent(inout) :: missing !< Missing value for the input data
 
   if (present(size)) then
-    size(1:4) = get_external_field_size(field_id)
+    size(1:4) = get_extern_field_size(field_id)
   endif
 
   if (present(axes)) then
-    axes(1:4) = get_external_field_axes(field_id)
+    axes(1:4) = get_extern_field_axes(field_id)
   endif
 
   if (present(missing)) then
-    missing = get_external_field_missing(field_id)
+    missing = get_extern_field_missing(field_id)
   endif
 
 end subroutine get_external_field_info
@@ -62,6 +186,7 @@ subroutine time_interp_extern_0d(field_id, time, data_in, verbose)
   call time_interp_external(field_id, time, data_in, verbose=verbose)
 end subroutine time_interp_extern_0d
 
+
 !> Read a 2d field from an external based on model time, potentially including horizontal
 !! interpolation and rotation of the data
 subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out)
@@ -98,6 +223,8 @@ subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_
                             horz_interp=horz_interp, mask_out=mask_out)
 end subroutine time_interp_extern_3d
 
+
+!> initialize an external field
 integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, &
                                    threading, ierr, ignore_axis_atts )
 

src/framework/MOM_interpolate.F90

@@ -7,14 +7,15 @@ module MOM_interpolate
 use MOM_error_handler,   only : MOM_error, FATAL
 use MOM_interp_infra,    only : time_interp_extern, init_external_field=>init_extern_field
 use MOM_interp_infra,    only : time_interp_external_init, get_external_field_info
-use MOM_interp_infra,    only : horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new
+use MOM_interp_infra,    only : horiz_interp_type, horiz_interp_init
+use MOM_interp_infra,    only : run_horiz_interp, build_horiz_interp_weights
 use MOM_io_infra,        only : axistype
 use MOM_time_manager,    only : time_type
 
 implicit none ; private
 
 public :: time_interp_external, init_external_field, time_interp_external_init, get_external_field_info
-public :: horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new
+public :: horiz_interp_type, horiz_interp_init, run_horiz_interp, build_horiz_interp_weights
 
 !> Read a field based on model time, and rotate to the model domain.
 interface time_interp_external