horiz_interp/include/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 !> @}