axis_utils/include/axis_utils2.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 !> @defgroup axis_utils2_mod axis_utils2_mod
  20 !> @ingroup axis_utils
  21 !> @brief A set of utilities for manipulating axes and extracting axis attributes.
  22 !! FMS2_IO equivalent version of @ref axis_utils_mod.
  23 !> @author M.J. Harrison
  24
  25 !> @addtogroup axis_utils2_mod
  26 !> @{
  27
  28   !> get axis edge data from a given file
  29   subroutine AXIS_EDGES_(fileobj, name, edge_data, reproduce_null_char_bug_flag)
  30
  31   class(FmsNetcdfFile_t), intent(in)            :: fileobj  !< File object to read from
  32   character(len=*), intent(in)                  :: name  !< Name of a given axis
  33   real(FMS_AU_KIND_), dimension(:), intent(out) :: edge_data  !< Returned edge data from given axis name
  34   logical, intent(in), optional                 :: reproduce_null_char_bug_flag  !< Flag indicating to reproduce
  35                                      !! the mpp_io bug where the null characters were not removed
  36                                      !! after reading a string attribute
  37
  38   integer                                              :: ndims
  39   character(len=128)                                   :: buffer
  40   integer,                 dimension(:),   allocatable :: dim_sizes
  41   real(kind=FMS_AU_KIND_), dimension(:),   allocatable :: r_var
  42   real(kind=FMS_AU_KIND_), dimension(:,:), allocatable :: r2d
  43   integer                                              :: i
  44   integer                                              :: n
  45   logical                                              :: reproduce_null_char_bug  !< Local flag
  46                                               !! indicating to reproduce the mpp_io bug where
  47                                               !! the null characters were not removed after reading a string attribute
  48   integer, parameter                                   :: lkind = FMS_AU_KIND_
  49   integer :: edge_index(2) !< Index to use when reading the edges from the file
  50                            !! (/1, 2/) if the axis data is monotonically increasing
  51                            !! (/2, 1/) if the axis data is monotonically decreasing
  52
  53   ndims = get_variable_num_dimensions(fileobj, name)
  54   allocate(dim_sizes(ndims))
  55
  56   call get_variable_size(fileobj, name, dim_sizes)
  57
  58   n = dim_sizes(1)
  59   if (size(edge_data) .ne. n+1) then
  60     call mpp_error(FATAL, "axis_edge: incorrect size of edge_data array.")
  61   endif
  62   deallocate(dim_sizes)
  63
  64   reproduce_null_char_bug = .false.
  65   if (present(reproduce_null_char_bug_flag)) reproduce_null_char_bug = reproduce_null_char_bug_flag
  66
  67   buffer = ""
  68   if (variable_att_exists(fileobj, name, "edges")) then
  69    !! If the reproduce_null_char_bug flag is turned on fms2io will not remove the null character
  70     call get_variable_attribute(fileobj, name, "edges", buffer(1:128), &
  71         reproduce_null_char_bug_flag=reproduce_null_char_bug)
  72
  73    !! Check for a null character here, if it exists *_bnds will be calculated instead of read in
  74     if (reproduce_null_char_bug) then
  75         i = 0
  76         i = index(buffer, char(0))
  77         if (i > 0) buffer = ""
  78     endif
  79   elseif (variable_att_exists(fileobj, name, "bounds")) then
  80    !! If the reproduce_null_char_bug flag is turned on fms2io will not remove the null character
  81     call get_variable_attribute(fileobj, name, "bounds", buffer(1:128), &
  82         reproduce_null_char_bug_flag=reproduce_null_char_bug)
  83
  84     !! Check for a null character here, if it exists *_bnds will be calculated instead of read in
  85     if (reproduce_null_char_bug) then
  86         i = 0
  87         i = index(buffer, char(0))
  88         if (i > 0) buffer = ""
  89     endif
  90   endif
  91   if (trim(buffer) .ne. "") then
  92     ndims = get_variable_num_dimensions(fileobj, buffer)
  93     allocate(dim_sizes(ndims))
  94
  95     call get_variable_size(fileobj, buffer, dim_sizes)
  96
  97     if (size(dim_sizes) .eq. 1) then
  98       if (dim_sizes(1) .ne. n+1) then
  99         call mpp_error(FATAL, "axis_edges: incorrect size of edge data.")
 100       endif
 101
 102       call read_data(fileobj, buffer, edge_data)
 103
 104     elseif (size(dim_sizes) .eq. 2) then
 105       if (dim_sizes(1) .ne. 2) then
 106         call mpp_error(FATAL, "axis_edges: first dimension of edge must be of size 2")
 107       endif
 108       if (dim_sizes(2) .ne. n) then
 109         call mpp_error(FATAL, "axis_edges: incorrect size of edge data.")
 110       endif
 111
 112       allocate(r2d(dim_sizes(1), dim_sizes(2)))
 113       call read_data(fileobj, buffer, r2d)
 114       edge_index = (/1, 2/)
 115       if (r2d(1,1) .gt. r2d(1,2)) edge_index = (/2, 1 /)
 116       edge_data(1:dim_sizes(2)) = r2d(edge_index(1),:)
 117       edge_data(dim_sizes(2)+1) = r2d(edge_index(2),dim_sizes(2))
 118       deallocate(r2d)
 119     endif
 120     deallocate(dim_sizes)
 121   else
 122       allocate(r_var(n))
 123
 124       call read_data(fileobj, name, r_var)
 125
 126       do i = 2, n
 127          edge_data(i) = r_var(i-1) + 0.5_lkind*(r_var(i) - r_var(i-1))
 128       enddo
 129       edge_data(1)    = r_var(1) - 0.5_lkind*(r_var(2) - r_var(1))
 130       if (abs(edge_data(1)) .lt. 1.e-10_lkind) then
 131          edge_data(1) = 0.0_lkind
 132       endif
 133       edge_data(n+1)  = r_var(n) + 0.5_lkind*(r_var(n) - r_var(n-1))
 134       deallocate(r_var)
 135   endif
 136   end subroutine AXIS_EDGES_
 137
 138   !> @brief Returns lon_strt <= longitude <= lon_strt+360
 139   !! @return real lon_in_range */
 140
 141   function LON_IN_RANGE_(lon, l_strt)
 142     real(kind=FMS_AU_KIND_), intent(in) :: lon, l_strt
 143     real(kind=FMS_AU_KIND_)             :: LON_IN_RANGE_
 144     real(kind=FMS_AU_KIND_)             :: l_end
 145     integer, parameter                  :: lkind = FMS_AU_KIND_
 146
 147     LON_IN_RANGE_ = lon
 148     l_end = l_strt + 360.0_lkind
 149
 150     if (abs(LON_IN_RANGE_ - l_strt) < 1.e-4_lkind) then
 151       LON_IN_RANGE_ = l_strt
 152       return
 153     endif
 154
 155     if (abs(LON_IN_RANGE_ - l_end) < 1.e-4_lkind) then
 156       LON_IN_RANGE_ = l_strt
 157       return
 158     endif
 159
 160     do
 161       if (LON_IN_RANGE_ < l_strt) then
 162         LON_IN_RANGE_ = real(LON_IN_RANGE_, FMS_AU_KIND_) + real(f360, FMS_AU_KIND_)
 163       else if (LON_IN_RANGE_  >  l_end) then
 164         LON_IN_RANGE_ = real(LON_IN_RANGE_, FMS_AU_KIND_) - real(f360, FMS_AU_KIND_)
 165       else
 166         exit
 167       end if
 168     end do
 169
 170   end function LON_IN_RANGE_
 171
 172   !> @brief Returns monotonic array of longitudes s.t., lon_strt <= lon(:) < lon_strt+360.
 173   !!
 174   !! This may require that entries be moved from the beginning of the array to
 175   !! the end. If no entries are moved (i.e., if lon(:) is already monotonic in
 176   !! the range from lon_start to lon_start + 360), then istrt is set to 0. If
 177   !! any entries are moved, then istrt is set to the original index of the entry
 178   !! which becomes lon(1).
 179   !!
 180   !! e.g.,
 181   !!
 182   !!        lon = 0 1 2 3 4 5  ...  358 359; lon_strt = 3
 183   !!    ==> lon = 3 4 5 6 7 8  ...  359 360 361 362; istrt = 4
 184   !!
 185   subroutine TRANLON_(lon, lon_start, istrt)
 186     real(kind=FMS_AU_KIND_), intent(inout), dimension(:) :: lon
 187     real(kind=FMS_AU_KIND_), intent(in)                  :: lon_start
 188     integer, intent(out)                                 :: istrt
 189     integer                 :: len, i
 190     real(kind=FMS_AU_KIND_) :: lon_strt, tmp(size(lon(:))-1)
 191
 192     len = size(lon(:))
 193
 194     do i = 1, len
 195        lon(i) = lon_in_range(lon(i),lon_start)
 196     enddo
 197
 198     istrt = 0
 199     do i = 1,len-1
 200        if (lon(i+1) < lon(i)) then
 201           istrt = i+1
 202           exit
 203        endif
 204     enddo
 205
 206     if (istrt>1) then ! grid is not monotonic
 207        if (abs(lon(len)-lon(1)) < real(epsln, FMS_AU_KIND_)) then
 208           tmp = cshift(lon(1:len-1),istrt-1)
 209           lon(1:len-1) = tmp
 210           lon(len)     = lon(1)
 211        else
 212           lon = cshift(lon,istrt-1)
 213        endif
 214
 215        lon_strt = lon(1)
 216        do i=2,len
 217           lon(i)   = lon_in_range(lon(i),lon_strt)
 218           lon_strt = lon(i)
 219        enddo
 220     endif
 221
 222     return
 223   end subroutine TRANLON_
 224
 225
 226   function FRAC_INDEX_(rval, array)
 227
 228     integer                               :: ia, i, ii, iunit
 229     real(kind=FMS_AU_KIND_)               :: rval !< arbitrary data...same units as elements in "array"
 230     real(kind=FMS_AU_KIND_)               :: FRAC_INDEX_
 231     real(kind=FMS_AU_KIND_), dimension(:) :: array !< array of data points  (must be monotonically increasing)
 232     logical                               :: keep_going
 233     integer, parameter                    :: lkind = FMS_AU_KIND_
 234     ia = size(array(:))
 235
 236     do i = 2, ia
 237        if (array(i) < array(i-1)) then
 238           iunit = stdout()
 239           write (iunit,*) '=> Error: "frac_index" array must be monotonically' &
 240           & // 'increasing when searching for nearest value to ', rval
 241           write (iunit,*) '          array(i) < array(i-1) for i=',i
 242           write (iunit,*) '          array(i) for i=1..ia follows:'
 243           do ii = 1, ia
 244              write (iunit,*) 'i=',ii, ' array(i)=',array(ii)
 245           enddo
 246           call mpp_error(FATAL,' "frac_index" array must be monotonically increasing.')
 247        endif
 248     enddo
 249
 250     if (rval < array(1) .or. rval > array(ia)) then
 251         FRAC_INDEX_ = -1.0_lkind
 252     else
 253        i = 1
 254        keep_going = .true.
 255        do while (i <= ia .and. keep_going)
 256           i = i+1
 257           if (rval <= array(i)) then
 258              FRAC_INDEX_ = real((i-1), lkind) + (rval-array(i-1)) / (array(i) - array(i-1))
 259              keep_going  = .false.
 260           endif
 261        enddo
 262     endif
 263   end function FRAC_INDEX_
 264
 265   !> @brief Return index of nearest point along axis
 266   !!
 267   !>     nearest_index = index of nearest data point within "array" corresponding to
 268   !!            "value".
 269   !!
 270   !!     inputs:
 271   !!
 272   !!     rval   = arbitrary data...same units as elements in "array"
 273   !!     array  = array of data points  (must be monotonically)
 274   !!     ia     = dimension of "array"
 275   !!
 276   !!     output:
 277   !!
 278   !!     nearest_index =  index of nearest data point to "value"
 279   !!             if "value" is outside the domain of "array" then nearest_index = 1
 280   !!             or "ia" depending on whether array(1) or array(ia) is
 281   !!             closest to "value"
 282   !!
 283   !!             note: if "array" is dimensioned array(0:ia) in the calling
 284   !!                   program, then the returned index should be reduced
 285   !!                   by one to account for the zero base.
 286   !!
 287   !!     example:
 288   !!
 289   !!     let model depths be defined by the following:
 290   !!     parameter (km=5)
 291   !!     dimension z(km)
 292   !!     data z /5.0, 10.0, 50.0, 100.0, 250.0/
 293   !!
 294   !!     k1 = nearest_index (12.5, z, km)
 295   !!     k2 = nearest_index (0.0, z, km)
 296   !!
 297   !!     k1 would be set to 2, and k2 would be set to 1 so that
 298   !!     z(k1) would be the nearest data point to 12.5 and z(k2) would
 299   !!     be the nearest data point to 0.0
 300   !! @return integer nearest_index
 301   function NEAREST_INDEX_(rval, array)
 302     real(kind=FMS_AU_KIND_), intent(in)               :: rval !< arbitrary data...same units as elements in "array"
 303     real(kind=FMS_AU_KIND_), intent(in), dimension(:) :: array !< array of data points  (must be monotonic)
 304
 305     integer                               :: NEAREST_INDEX_
 306     integer                               :: ia !< dimension of "array"
 307     integer                               :: i !< For looping through "array"
 308
 309     logical :: increasing !< .True. if the array is increasing
 310
 311     ia = SIZE(array(:))
 312
 313     ! check if array is increasing
 314     increasing = .true.
 315     DO i = 2, ia-1
 316        IF( array(i) .lt. array(i-1)) then
 317          increasing = .false.
 318          exit
 319        endif
 320     END DO
 321
 322     if (.not. increasing) then
 323       ! if not increasing, check that it is decreasing
 324       DO i = 2, ia-1
 325          IF( array(i) .gt. array(i-1)) &
 326             call mpp_error(FATAL, 'axis_utils2::nearest_index array is NOT monotonously ordered')
 327       END DO
 328     endif
 329
 330     array_is_increasing: if (increasing) then
 331       !< Check if the rval is outside the range of the array
 332       if (rval .le. array(1)) then
 333         NEAREST_INDEX_ = 1
 334         return
 335       elseif (rval .ge. array(ia)) then
 336         NEAREST_INDEX_ = ia
 337         return
 338       endif
 339
 340       DO i = 2, ia
 341         if (rval .le. array(i)) then
 342           NEAREST_INDEX_ = i
 343           if (array(i) -rval .gt. rval - array(i-1)) NEAREST_INDEX_ = i - 1
 344           return
 345         endif
 346       END DO
 347     else !array_is_decreasing
 348       !< Check if the rval is outside the range of the array
 349       if (rval .le. array(ia)) then
 350         NEAREST_INDEX_ = ia
 351         return
 352       elseif (rval .gt. array(1)) then
 353         NEAREST_INDEX_ = 1
 354         return
 355       endif
 356
 357       DO i = 2, ia
 358         if (rval .ge. array(i)) then
 359           NEAREST_INDEX_ = i
 360           if (rval - array(i) .gt. array(i-1) -rval ) NEAREST_INDEX_ = i - 1
 361             return
 362         endif
 363       END DO
 364     endif array_is_increasing
 365   end function NEAREST_INDEX_
 366
 367   !#############################################################################
 368
 369   subroutine INTERP_1D_LINEAR_(grid1,grid2,data1,data2)
 370
 371     real(kind=FMS_AU_KIND_), dimension(:),    intent(in) :: grid1, data1, grid2
 372     real(kind=FMS_AU_KIND_), dimension(:), intent(inout) :: data2
 373
 374     integer                 :: n1, n2, i, n
 375     real(kind=FMS_AU_KIND_) :: w
 376     integer, parameter      :: lkind = FMS_AU_KIND_
 377
 378     n1 = size(grid1(:))
 379     n2 = size(grid2(:))
 380
 381
 382     do i = 2, n1
 383        if (grid1(i) <= grid1(i-1)) call mpp_error(FATAL, 'grid1 not monotonic')
 384     enddo
 385
 386     do i = 2, n2
 387        if (grid2(i) <= grid2(i-1)) call mpp_error(FATAL, 'grid2 not monotonic')
 388     enddo
 389
 390     if (grid1(1) > grid2(1) ) call mpp_error(FATAL, 'grid2 lies outside grid1')
 391     if (grid1(n1) < grid2(n2) ) call mpp_error(FATAL, 'grid2 lies outside grid1')
 392
 393     do i = 1, n2
 394        n = nearest_index(grid2(i),grid1)
 395
 396        if (grid1(n) < grid2(i)) then
 397           w        = (grid2(i)-grid1(n))/(grid1(n+1)-grid1(n))
 398           data2(i) = (1.0_lkind-w)*data1(n) + w*data1(n+1)
 399        else
 400           if(n==1) then
 401              data2(i) = data1(n)
 402           else
 403              w        = (grid2(i)-grid1(n-1))/(grid1(n)-grid1(n-1))
 404              data2(i) = (1.0_lkind-w)*data1(n-1) + w*data1(n)
 405           endif
 406        endif
 407     enddo
 408
 409
 410     return
 411
 412   end subroutine INTERP_1D_LINEAR_
 413
 414   !###################################################################
 415   subroutine INTERP_1D_CUBIC_SPLINE_(grid1, grid2, data1, data2, yp1, ypn)
 416
 417     real(kind=FMS_AU_KIND_), dimension(:),    intent(in)    :: grid1, grid2, data1
 418     real(kind=FMS_AU_KIND_), dimension(:),    intent(inout) :: data2
 419     real(kind=FMS_AU_KIND_),                  intent(in)    :: yp1, ypn
 420
 421     real(kind=FMS_AU_KIND_), dimension(size(grid1))         :: y2, u
 422     real(kind=FMS_AU_KIND_)                                 :: sig, p, qn, un, h, a ,b
 423     integer                                                 :: n, m, i, k, klo, khi
 424     integer, parameter                                      :: lkind = FMS_AU_KIND_
 425
 426     n = size(grid1(:))
 427     m = size(grid2(:))
 428
 429     do i = 2, n
 430        if (grid1(i) <= grid1(i-1)) call mpp_error(FATAL, 'grid1 not monotonic')
 431     enddo
 432
 433     do i = 2, m
 434        if (grid2(i) <= grid2(i-1)) call mpp_error(FATAL, 'grid2 not monotonic')
 435     enddo
 436
 437     if (grid1(1) > grid2(1) ) call mpp_error(FATAL, 'grid2 lies outside grid1')
 438     if (grid1(n) < grid2(m) ) call mpp_error(FATAL, 'grid2 lies outside grid1')
 439
 440 if (yp1>0.99e30_lkind) then
 441        y2(1) = 0.0_lkind
 442        u(1)  = 0.0_lkind
 443     else
 444        y2(1) = -0.5_lkind
 445        u(1)  = (3.0_lkind)/(grid1(2)-grid1(1))*((data1(2)-data1(1))/(grid1(2)-grid1(1))-yp1)
 446     endif
 447
 448     do i = 2, n-1
 449        sig   = (grid1(i)-grid1(i-1))/(grid1(i+1)-grid1(i-1))
 450        p     = sig*y2(i-1) + 2.0_lkind
 451        y2(i) = (sig-1.0_lkind)/p
 452        u(i)  = (6.0_lkind*((data1(i+1)-data1(i))/(grid1(i+1)-grid1(i))-(data1(i)-data1(i-1)) &
 453                /(grid1(i)-grid1(i-1)))/(grid1(i+1)-grid1(i-1))-sig*u(i-1))/p
 454     enddo
 455
 456     if (ypn>0.99e30_lkind) then
 457        qn = 0.0_lkind
 458        un = 0.0_lkind
 459     else
 460        qn = 0.5_lkind
 461        un = (3.0_lkind)/(grid1(n)-grid1(n-1))*(ypn-(data1(n)-data1(n-1))/ &
 462             (grid1(n)-grid1(n-1)))
 463     endif
 464
 465     y2(n) = (un-qn*u(n-1))/(qn*y2(n-1)+1.0_lkind)
 466
 467     do  k = n-1,1,-1
 468        y2(k) = y2(k)*y2(k+1)+u(k)
 469     enddo
 470
 471     do k = 1, m
 472        n = nearest_index(grid2(k),grid1)
 473        if (grid1(n) < grid2(k)) then
 474           klo = n
 475        else
 476           if(n==1) then
 477             klo = n
 478           else
 479             klo = n -1
 480           endif
 481        endif
 482
 483        khi      = klo+1
 484        h        = grid1(khi)-grid1(klo)
 485        a        = (grid1(khi) - grid2(k))/h
 486        b        = (grid2(k) - grid1(klo))/h
 487        data2(k) = a*data1(klo) + b*data1(khi)+ ((a**3-a)*y2(klo) + (b**3-b)*y2(khi))*(h**2) &
 488                  /6.0_lkind
 489     enddo
 490
 491   end subroutine INTERP_1D_CUBIC_SPLINE_
 492
 493   !###################################################################
 494
 495   subroutine INTERP_1D_1D_(grid1,grid2,data1,data2, method, yp1, yp2)
 496
 497     real(kind=FMS_AU_KIND_), dimension(:),         intent(in)    :: grid1, data1, grid2
 498     real(kind=FMS_AU_KIND_), dimension(:),         intent(inout) :: data2
 499     character(len=*),                    optional, intent(in)    :: method
 500     real(kind=FMS_AU_KIND_),             optional, intent(in)    :: yp1, yp2
 501
 502     real(kind=FMS_AU_KIND_)                                      :: y1, y2
 503     character(len=32)                                            :: interp_method
 504     integer                                                      :: k2, ks, ke
 505     integer, parameter                                           :: lkind = FMS_AU_KIND_
 506
 507     k2 = size(grid2(:))
 508
 509     interp_method = "linear"
 510     if(present(method)) interp_method = method
 511     y1 = 1.0e30_lkind
 512
 513     if(present(yp1)) y1 = yp1
 514     y2 = 1.0e30_lkind
 515
 516     if(present(yp2)) y2 = yp2
 517     call find_index(grid1, grid2(1), grid2(k2), ks, ke)
 518     select case(trim(interp_method))
 519     case("linear")
 520        call interp_1d_linear(grid1(ks:ke),grid2,data1(ks:ke),data2)
 521     case("cubic_spline")
 522        call interp_1d_cubic_spline(grid1(ks:ke),grid2,data1(ks:ke),data2, y1, y2)
 523     case default
 524        call mpp_error(FATAL,"axis_utils: interp_method should be linear or cubic_spline")
 525     end select
 526
 527     return
 528
 529   end subroutine INTERP_1D_1D_
 530
 531  !###################################################################
 532
 533
 534   subroutine INTERP_1D_2D_(grid1,grid2,data1,data2)
 535
 536     real(kind=FMS_AU_KIND_), dimension(:,:),    intent(in) :: grid1, data1, grid2
 537     real(kind=FMS_AU_KIND_), dimension(:,:), intent(inout) :: data2
 538
 539     integer :: n1, n2, n, k2, ks, ke
 540
 541     n1 = size(grid1,1)
 542     n2 = size(grid2,1)
 543     k2 = size(grid2,2)
 544
 545     if (n1 /= n2) call mpp_error(FATAL,'grid size mismatch')
 546
 547     do n = 1, n1
 548        call find_index(grid1(n,:), grid2(n,1), grid2(n,k2), ks, ke)
 549        call interp_1d_linear(grid1(n,ks:ke),grid2(n,:),data1(n,ks:ke),data2(n,:))
 550     enddo
 551
 552     return
 553
 554   end subroutine INTERP_1D_2D_
 555
 556   !###################################################################
 557
 558   subroutine INTERP_1D_3D_(grid1,grid2,data1,data2, method, yp1, yp2)
 559
 560     real(FMS_AU_KIND_), dimension(:,:,:), intent(in)    :: grid1, data1, grid2
 561     real(FMS_AU_KIND_), dimension(:,:,:), intent(inout) :: data2
 562     character(len=*),           optional, intent(in)    :: method
 563     real(kind=FMS_AU_KIND_),    optional, intent(in)    :: yp1, yp2
 564
 565     integer                                             :: n1, n2, m1, m2, k2, n, m
 566     real(kind=FMS_AU_KIND_)                             :: y1, y2
 567     character(len=32)                                   :: interp_method
 568     integer                                             :: ks, ke
 569     integer, parameter                                  :: lkind = FMS_AU_KIND_
 570
 571     n1 = size(grid1,1)
 572     n2 = size(grid2,1)
 573     m1 = size(grid1,2)
 574     m2 = size(grid2,2)
 575     k2 = size(grid2,3)
 576
 577     interp_method = "linear"
 578     if(present(method)) interp_method = method
 579     y1 = 1.0e30_lkind
 580
 581     if(present(yp1)) y1 = yp1
 582     y2 = 1.0e30_lkind
 583     if(present(yp2)) y2 = yp2
 584
 585     if (n1 /= n2 .or. m1 /= m2) call mpp_error(FATAL,'grid size mismatch')
 586
 587     select case(trim(interp_method))
 588     case("linear")
 589        do m = 1, m1
 590           do n = 1, n1
 591             call find_index(grid1(n,m,:), grid2(n,m,1), grid2(n,m,k2), ks, ke)
 592             call interp_1d_linear(grid1(n,m,ks:ke),grid2(n,m,:),data1(n,m,ks:ke),data2(n,m,:))
 593           enddo
 594        enddo
 595
 596     case("cubic_spline")
 597        do m = 1, m1
 598           do n = 1, n1
 599             call find_index(grid1(n,m,:), grid2(n,m,1), grid2(n,m,k2), ks, ke)
 600             call interp_1d_cubic_spline(grid1(n,m,ks:ke),grid2(n,m,:), data1(n,m,ks:ke),data2(n,m,:), y1, y2)
 601           enddo
 602        enddo
 603
 604     case default
 605        call mpp_error(FATAL,"axis_utils: interp_method should be linear or cubic_spline")
 606     end select
 607
 608     return
 609
 610   end subroutine INTERP_1D_3D_
 611
 612
 613   !#####################################################################
 614   subroutine FIND_INDEX_(grid1, xs, xe, ks, ke)
 615     real(kind=FMS_AU_KIND_), dimension(:), intent(in)  :: grid1
 616     real(kind=FMS_AU_KIND_),               intent(in)  :: xs, xe
 617     integer,                               intent(out) :: ks, ke
 618
 619     integer :: k, nk
 620
 621     nk = size(grid1(:))
 622
 623     ks = 0; ke = 0
 624     do k = 1, nk-1
 625        if(grid1(k) <= xs .and. grid1(k+1) > xs ) then
 626           ks = k
 627           exit
 628        endif
 629     enddo
 630
 631     do k = nk, 2, -1
 632        if(grid1(k) >= xe .and. grid1(k-1) < xe ) then
 633           ke = k
 634           exit
 635        endif
 636     enddo
 637
 638     if(ks == 0 ) call mpp_error(FATAL,' xs locate outside of grid1')
 639     if(ke == 0 ) call mpp_error(FATAL,' xe locate outside of grid1')
 640
 641   end subroutine FIND_INDEX_
 642  !> @}
 643  ! close documentation grouping