axis_utils2.inc

!***********************************************************************
!*                             Apache License 2.0
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* Licensed under the Apache License, Version 2.0 (the "License");
!* you may not use this file except in compliance with the License.
!* You may obtain a copy of the License at
!*
!*     http://www.apache.org/licenses/LICENSE-2.0
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
!* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
!* PARTICULAR PURPOSE. See the License for the specific language
!* governing permissions and limitations under the License.
!***********************************************************************
!> @defgroup axis_utils2_mod axis_utils2_mod
!> @ingroup axis_utils
!> @brief A set of utilities for manipulating axes and extracting axis attributes.
!! FMS2_IO equivalent version of @ref axis_utils_mod.
!> @author M.J. Harrison
 
!> @addtogroup axis_utils2_mod
!> @{
 
  !> get axis edge data from a given file
  subroutine axis_edges_(fileobj, name, edge_data, reproduce_null_char_bug_flag)
 
  class(fmsnetcdffile_t), intent(in)            :: fileobj  !< File object to read from
  character(len=*), intent(in)                  :: name  !< Name of a given axis
  real(FMS_AU_KIND_), dimension(:), intent(out) :: edge_data  !< Returned edge data from given axis name
  logical, intent(in), optional                 :: reproduce_null_char_bug_flag  !< Flag indicating to reproduce
                                     !! the mpp_io bug where the null characters were not removed
                                     !! after reading a string attribute
 
  integer                                              :: ndims
  character(len=128)                                   :: buffer
  integer,                 dimension(:),   allocatable :: dim_sizes
  real(kind=fms_au_kind_), dimension(:),   allocatable :: r_var
  real(kind=fms_au_kind_), dimension(:,:), allocatable :: r2d
  integer                                              :: i
  integer                                              :: n
  logical                                              :: reproduce_null_char_bug  !< Local flag
                                              !! indicating to reproduce the mpp_io bug where
                                              !! the null characters were not removed after reading a string attribute
  integer, parameter                                   :: lkind = fms_au_kind_
  integer :: edge_index(2) !< Index to use when reading the edges from the file
                           !! (/1, 2/) if the axis data is monotonically increasing
                           !! (/2, 1/) if the axis data is monotonically decreasing
 
  ndims = get_variable_num_dimensions(fileobj, name)
  allocate(dim_sizes(ndims))
 
  call get_variable_size(fileobj, name, dim_sizes)
 
  n = dim_sizes(1)
  if (size(edge_data) .ne. n+1) then
    call mpp_error(fatal, "axis_edge: incorrect size of edge_data array.")
  endif
  deallocate(dim_sizes)
 
  reproduce_null_char_bug = .false.
  if (present(reproduce_null_char_bug_flag)) reproduce_null_char_bug = reproduce_null_char_bug_flag
 
  buffer = ""
  if (variable_att_exists(fileobj, name, "edges")) then
   !! If the reproduce_null_char_bug flag is turned on fms2io will not remove the null character
    call get_variable_attribute(fileobj, name, "edges", buffer(1:128), &
        reproduce_null_char_bug_flag=reproduce_null_char_bug)
 
   !! Check for a null character here, if it exists *_bnds will be calculated instead of read in
    if (reproduce_null_char_bug) then
        i = 0
        i = index(buffer, char(0))
        if (i > 0) buffer = ""
    endif
  elseif (variable_att_exists(fileobj, name, "bounds")) then
   !! If the reproduce_null_char_bug flag is turned on fms2io will not remove the null character
    call get_variable_attribute(fileobj, name, "bounds", buffer(1:128), &
        reproduce_null_char_bug_flag=reproduce_null_char_bug)
 
    !! Check for a null character here, if it exists *_bnds will be calculated instead of read in
    if (reproduce_null_char_bug) then
        i = 0
        i = index(buffer, char(0))
        if (i > 0) buffer = ""
    endif
  endif
  if (trim(buffer) .ne. "") then
    ndims = get_variable_num_dimensions(fileobj, buffer)
    allocate(dim_sizes(ndims))
 
    call get_variable_size(fileobj, buffer, dim_sizes)
 
    if (size(dim_sizes) .eq. 1) then
      if (dim_sizes(1) .ne. n+1) then
        call mpp_error(fatal, "axis_edges: incorrect size of edge data.")
      endif
 
      call read_data(fileobj, buffer, edge_data)
 
    elseif (size(dim_sizes) .eq. 2) then
      if (dim_sizes(1) .ne. 2) then
        call mpp_error(fatal, "axis_edges: first dimension of edge must be of size 2")
      endif
      if (dim_sizes(2) .ne. n) then
        call mpp_error(fatal, "axis_edges: incorrect size of edge data.")
      endif
 
      allocate(r2d(dim_sizes(1), dim_sizes(2)))
      call read_data(fileobj, buffer, r2d)
      edge_index = (/1, 2/)
      if (r2d(1,1) .gt. r2d(1,2)) edge_index = (/2, 1 /)
      edge_data(1:dim_sizes(2)) = r2d(edge_index(1),:)
      edge_data(dim_sizes(2)+1) = r2d(edge_index(2),dim_sizes(2))
      deallocate(r2d)
    endif
    deallocate(dim_sizes)
  else
      allocate(r_var(n))
 
      call read_data(fileobj, name, r_var)
 
      do i = 2, n
         edge_data(i) = r_var(i-1) + 0.5_lkind*(r_var(i) - r_var(i-1))
      enddo
      edge_data(1)    = r_var(1) - 0.5_lkind*(r_var(2) - r_var(1))
      if (abs(edge_data(1)) .lt. 1.e-10_lkind) then
         edge_data(1) = 0.0_lkind
      endif
      edge_data(n+1)  = r_var(n) + 0.5_lkind*(r_var(n) - r_var(n-1))
      deallocate(r_var)
  endif
  end subroutine axis_edges_
 
  !> @brief Returns lon_strt <= longitude <= lon_strt+360
  !! @return real lon_in_range */
 
  function lon_in_range_(lon, l_strt)
    real(kind=fms_au_kind_), intent(in) :: lon, l_strt
    real(kind=fms_au_kind_)             :: lon_in_range_
    real(kind=fms_au_kind_)             :: l_end
    integer, parameter                  :: lkind = fms_au_kind_
 
    lon_in_range_ = lon
    l_end = l_strt + 360.0_lkind
 
    if (abs(lon_in_range_ - l_strt) < 1.e-4_lkind) then
      lon_in_range_ = l_strt
      return
    endif
 
    if (abs(lon_in_range_ - l_end) < 1.e-4_lkind) then
      lon_in_range_ = l_strt
      return
    endif
 
    do
      if (lon_in_range_ < l_strt) then
        lon_in_range_ = real(lon_in_range_, fms_au_kind_) + real(f360, fms_au_kind_)
      else if (lon_in_range_  >  l_end) then
        lon_in_range_ = real(lon_in_range_, fms_au_kind_) - real(f360, fms_au_kind_)
      else
        exit
      end if
    end do
 
  end function lon_in_range_
 
  !> @brief Returns monotonic array of longitudes s.t., lon_strt <= lon(:) < lon_strt+360.
  !!
  !! This may require that entries be moved from the beginning of the array to
  !! the end. If no entries are moved (i.e., if lon(:) is already monotonic in
  !! the range from lon_start to lon_start + 360), then istrt is set to 0. If
  !! any entries are moved, then istrt is set to the original index of the entry
  !! which becomes lon(1).
  !!
  !! e.g.,
  !!
  !!        lon = 0 1 2 3 4 5  ...  358 359; lon_strt = 3
  !!    ==> lon = 3 4 5 6 7 8  ...  359 360 361 362; istrt = 4
  !!
  subroutine tranlon_(lon, lon_start, istrt)
    real(kind=fms_au_kind_), intent(inout), dimension(:) :: lon
    real(kind=fms_au_kind_), intent(in)                  :: lon_start
    integer, intent(out)                                 :: istrt
    integer                 :: len, i
    real(kind=fms_au_kind_) :: lon_strt, tmp(size(lon(:))-1)
 
    len = size(lon(:))
 
    do i = 1, len
       lon(i) = lon_in_range(lon(i),lon_start)
    enddo
 
    istrt = 0
    do i = 1,len-1
       if (lon(i+1) < lon(i)) then
          istrt = i+1
          exit
       endif
    enddo
 
    if (istrt>1) then ! grid is not monotonic
       if (abs(lon(len)-lon(1)) < real(epsln, fms_au_kind_)) then
          tmp = cshift(lon(1:len-1),istrt-1)
          lon(1:len-1) = tmp
          lon(len)     = lon(1)
       else
          lon = cshift(lon,istrt-1)
       endif
 
       lon_strt = lon(1)
       do i=2,len
          lon(i)   = lon_in_range(lon(i),lon_strt)
          lon_strt = lon(i)
       enddo
    endif
 
    return
  end subroutine tranlon_
 
 
  function frac_index_(rval, array)
 
    integer                               :: ia, i, ii, iunit
    real(kind=fms_au_kind_)               :: rval !< arbitrary data...same units as elements in "array"
    real(kind=fms_au_kind_)               :: frac_index_
    real(kind=fms_au_kind_), dimension(:) :: array !< array of data points  (must be monotonically increasing)
    logical                               :: keep_going
    integer, parameter                    :: lkind = fms_au_kind_
    ia = size(array(:))
 
    do i = 2, ia
       if (array(i) < array(i-1)) then
          iunit = stdout()
          write (iunit,*) '=> Error: "frac_index" array must be monotonically' &
          & // 'increasing when searching for nearest value to ', rval
          write (iunit,*) '          array(i) < array(i-1) for i=',i
          write (iunit,*) '          array(i) for i=1..ia follows:'
          do ii = 1, ia
             write (iunit,*) 'i=',ii, ' array(i)=',array(ii)
          enddo
          call mpp_error(fatal,' "frac_index" array must be monotonically increasing.')
       endif
    enddo
 
    if (rval < array(1) .or. rval > array(ia)) then
        frac_index_ = -1.0_lkind
    else
       i = 1
       keep_going = .true.
       do while (i <= ia .and. keep_going)
          i = i+1
          if (rval <= array(i)) then
             frac_index_ = real((i-1), lkind) + (rval-array(i-1)) / (array(i) - array(i-1))
             keep_going  = .false.
          endif
       enddo
    endif
  end function frac_index_
 
  !> @brief Return index of nearest point along axis
  !!
  !>     nearest_index = index of nearest data point within "array" corresponding to
  !!            "value".
  !!
  !!     inputs:
  !!
  !!     rval   = arbitrary data...same units as elements in "array"
  !!     array  = array of data points  (must be monotonically)
  !!     ia     = dimension of "array"
  !!
  !!     output:
  !!
  !!     nearest_index =  index of nearest data point to "value"
  !!             if "value" is outside the domain of "array" then nearest_index = 1
  !!             or "ia" depending on whether array(1) or array(ia) is
  !!             closest to "value"
  !!
  !!             note: if "array" is dimensioned array(0:ia) in the calling
  !!                   program, then the returned index should be reduced
  !!                   by one to account for the zero base.
  !!
  !!     example:
  !!
  !!     let model depths be defined by the following:
  !!     parameter (km=5)
  !!     dimension z(km)
  !!     data z /5.0, 10.0, 50.0, 100.0, 250.0/
  !!
  !!     k1 = nearest_index (12.5, z, km)
  !!     k2 = nearest_index (0.0, z, km)
  !!
  !!     k1 would be set to 2, and k2 would be set to 1 so that
  !!     z(k1) would be the nearest data point to 12.5 and z(k2) would
  !!     be the nearest data point to 0.0
  !! @return integer nearest_index
  function nearest_index_(rval, array)
    real(kind=fms_au_kind_), intent(in)               :: rval !< arbitrary data...same units as elements in "array"
    real(kind=fms_au_kind_), intent(in), dimension(:) :: array !< array of data points  (must be monotonic)
 
    integer                               :: nearest_index_
    integer                               :: ia !< dimension of "array"
    integer                               :: i !< For looping through "array"
 
    logical :: increasing !< .True. if the array is increasing
 
    ia = SIZE(array(:))
 
    ! check if array is increasing
    increasing = .true.
    DO i = 2, ia-1
       IF( array(i) .lt. array(i-1)) then
         increasing = .false.
         exit
       endif
    END DO
 
    if (.not. increasing) then
      ! if not increasing, check that it is decreasing
      DO i = 2, ia-1
         IF( array(i) .gt. array(i-1)) &
            call mpp_error(fatal, 'axis_utils2::nearest_index array is NOT monotonously ordered')
      END DO
    endif
 
    array_is_increasing: if (increasing) then
      !< Check if the rval is outside the range of the array
      if (rval .le. array(1)) then
        nearest_index_ = 1
        return
      elseif (rval .ge. array(ia)) then
        nearest_index_ = ia
        return
      endif
 
      DO i = 2, ia
        if (rval .le. array(i)) then
          nearest_index_ = i
          if (array(i) -rval .gt. rval - array(i-1)) nearest_index_ = i - 1
          return
        endif
      END DO
    else !array_is_decreasing
      !< Check if the rval is outside the range of the array
      if (rval .le. array(ia)) then
        nearest_index_ = ia
        return
      elseif (rval .gt. array(1)) then
        nearest_index_ = 1
        return
      endif
 
      DO i = 2, ia
        if (rval .ge. array(i)) then
          nearest_index_ = i
          if (rval - array(i) .gt. array(i-1) -rval ) nearest_index_ = i - 1
            return
        endif
      END DO
    endif array_is_increasing
  end function nearest_index_
 
  !#############################################################################
 
  subroutine interp_1d_linear_(grid1,grid2,data1,data2)
 
    real(kind=fms_au_kind_), dimension(:),    intent(in) :: grid1, data1, grid2
    real(kind=fms_au_kind_), dimension(:), intent(inout) :: data2
 
    integer                 :: n1, n2, i, n
    real(kind=fms_au_kind_) :: w
    integer, parameter      :: lkind = fms_au_kind_
 
    n1 = size(grid1(:))
    n2 = size(grid2(:))
 
 
    do i = 2, n1
       if (grid1(i) <= grid1(i-1)) call mpp_error(fatal, 'grid1 not monotonic')
    enddo
 
    do i = 2, n2
       if (grid2(i) <= grid2(i-1)) call mpp_error(fatal, 'grid2 not monotonic')
    enddo
 
    if (grid1(1) > grid2(1) ) call mpp_error(fatal, 'grid2 lies outside grid1')
    if (grid1(n1) < grid2(n2) ) call mpp_error(fatal, 'grid2 lies outside grid1')
 
    do i = 1, n2
       n = nearest_index(grid2(i),grid1)
 
       if (grid1(n) < grid2(i)) then
          w        = (grid2(i)-grid1(n))/(grid1(n+1)-grid1(n))
          data2(i) = (1.0_lkind-w)*data1(n) + w*data1(n+1)
       else
          if(n==1) then
             data2(i) = data1(n)
          else
             w        = (grid2(i)-grid1(n-1))/(grid1(n)-grid1(n-1))
             data2(i) = (1.0_lkind-w)*data1(n-1) + w*data1(n)
          endif
       endif
    enddo
 
 
    return
 
  end subroutine interp_1d_linear_
 
  !###################################################################
  subroutine interp_1d_cubic_spline_(grid1, grid2, data1, data2, yp1, ypn)
 
    real(kind=fms_au_kind_), dimension(:),    intent(in)    :: grid1, grid2, data1
    real(kind=fms_au_kind_), dimension(:),    intent(inout) :: data2
    real(kind=fms_au_kind_),                  intent(in)    :: yp1, ypn
 
    real(kind=fms_au_kind_), dimension(size(grid1))         :: y2, u
    real(kind=fms_au_kind_)                                 :: sig, p, qn, un, h, a ,b
    integer                                                 :: n, m, i, k, klo, khi
    integer, parameter                                      :: lkind = fms_au_kind_
 
    n = size(grid1(:))
    m = size(grid2(:))
 
    do i = 2, n
       if (grid1(i) <= grid1(i-1)) call mpp_error(fatal, 'grid1 not monotonic')
    enddo
 
    do i = 2, m
       if (grid2(i) <= grid2(i-1)) call mpp_error(fatal, 'grid2 not monotonic')
    enddo
 
    if (grid1(1) > grid2(1) ) call mpp_error(fatal, 'grid2 lies outside grid1')
    if (grid1(n) < grid2(m) ) call mpp_error(fatal, 'grid2 lies outside grid1')
 
if (yp1>0.99e30_lkind) then
       y2(1) = 0.0_lkind
       u(1)  = 0.0_lkind
    else
       y2(1) = -0.5_lkind
       u(1)  = (3.0_lkind)/(grid1(2)-grid1(1))*((data1(2)-data1(1))/(grid1(2)-grid1(1))-yp1)
    endif
 
    do i = 2, n-1
       sig   = (grid1(i)-grid1(i-1))/(grid1(i+1)-grid1(i-1))
       p     = sig*y2(i-1) + 2.0_lkind
       y2(i) = (sig-1.0_lkind)/p
       u(i)  = (6.0_lkind*((data1(i+1)-data1(i))/(grid1(i+1)-grid1(i))-(data1(i)-data1(i-1)) &
               /(grid1(i)-grid1(i-1)))/(grid1(i+1)-grid1(i-1))-sig*u(i-1))/p
    enddo
 
    if (ypn>0.99e30_lkind) then
       qn = 0.0_lkind
       un = 0.0_lkind
    else
       qn = 0.5_lkind
       un = (3.0_lkind)/(grid1(n)-grid1(n-1))*(ypn-(data1(n)-data1(n-1))/ &
            (grid1(n)-grid1(n-1)))
    endif
 
    y2(n) = (un-qn*u(n-1))/(qn*y2(n-1)+1.0_lkind)
 
    do  k = n-1,1,-1
       y2(k) = y2(k)*y2(k+1)+u(k)
    enddo
 
    do k = 1, m
       n = nearest_index(grid2(k),grid1)
       if (grid1(n) < grid2(k)) then
          klo = n
       else
          if(n==1) then
            klo = n
          else
            klo = n -1
          endif
       endif
 
       khi      = klo+1
       h        = grid1(khi)-grid1(klo)
       a        = (grid1(khi) - grid2(k))/h
       b        = (grid2(k) - grid1(klo))/h
       data2(k) = a*data1(klo) + b*data1(khi)+ ((a**3-a)*y2(klo) + (b**3-b)*y2(khi))*(h**2) &
                 /6.0_lkind
    enddo
 
  end subroutine interp_1d_cubic_spline_
 
  !###################################################################
 
  subroutine interp_1d_1d_(grid1,grid2,data1,data2, method, yp1, yp2)
 
    real(kind=fms_au_kind_), dimension(:),         intent(in)    :: grid1, data1, grid2
    real(kind=fms_au_kind_), dimension(:),         intent(inout) :: data2
    character(len=*),                    optional, intent(in)    :: method
    real(kind=fms_au_kind_),             optional, intent(in)    :: yp1, yp2
 
    real(kind=fms_au_kind_)                                      :: y1, y2
    character(len=32)                                            :: interp_method
    integer                                                      :: k2, ks, ke
    integer, parameter                                           :: lkind = fms_au_kind_
 
    k2 = size(grid2(:))
 
    interp_method = "linear"
    if(present(method)) interp_method = method
    y1 = 1.0e30_lkind
 
    if(present(yp1)) y1 = yp1
    y2 = 1.0e30_lkind
 
    if(present(yp2)) y2 = yp2
    call find_index(grid1, grid2(1), grid2(k2), ks, ke)
    select case(trim(interp_method))
    case("linear")
       call interp_1d_linear(grid1(ks:ke),grid2,data1(ks:ke),data2)
    case("cubic_spline")
       call interp_1d_cubic_spline(grid1(ks:ke),grid2,data1(ks:ke),data2, y1, y2)
    case default
       call mpp_error(fatal,"axis_utils: interp_method should be linear or cubic_spline")
    end select
 
    return
 
  end subroutine interp_1d_1d_
 
 !###################################################################
 
 
  subroutine interp_1d_2d_(grid1,grid2,data1,data2)
 
    real(kind=fms_au_kind_), dimension(:,:),    intent(in) :: grid1, data1, grid2
    real(kind=fms_au_kind_), dimension(:,:), intent(inout) :: data2
 
    integer :: n1, n2, n, k2, ks, ke
 
    n1 = size(grid1,1)
    n2 = size(grid2,1)
    k2 = size(grid2,2)
 
    if (n1 /= n2) call mpp_error(fatal,'grid size mismatch')
 
    do n = 1, n1
       call find_index(grid1(n,:), grid2(n,1), grid2(n,k2), ks, ke)
       call interp_1d_linear(grid1(n,ks:ke),grid2(n,:),data1(n,ks:ke),data2(n,:))
    enddo
 
    return
 
  end subroutine interp_1d_2d_
 
  !###################################################################
 
  subroutine interp_1d_3d_(grid1,grid2,data1,data2, method, yp1, yp2)
 
    real(fms_au_kind_), dimension(:,:,:), intent(in)    :: grid1, data1, grid2
    real(fms_au_kind_), dimension(:,:,:), intent(inout) :: data2
    character(len=*),           optional, intent(in)    :: method
    real(kind=fms_au_kind_),    optional, intent(in)    :: yp1, yp2
 
    integer                                             :: n1, n2, m1, m2, k2, n, m
    real(kind=fms_au_kind_)                             :: y1, y2
    character(len=32)                                   :: interp_method
    integer                                             :: ks, ke
    integer, parameter                                  :: lkind = fms_au_kind_
 
    n1 = size(grid1,1)
    n2 = size(grid2,1)
    m1 = size(grid1,2)
    m2 = size(grid2,2)
    k2 = size(grid2,3)
 
    interp_method = "linear"
    if(present(method)) interp_method = method
    y1 = 1.0e30_lkind
 
    if(present(yp1)) y1 = yp1
    y2 = 1.0e30_lkind
    if(present(yp2)) y2 = yp2
 
    if (n1 /= n2 .or. m1 /= m2) call mpp_error(fatal,'grid size mismatch')
 
    select case(trim(interp_method))
    case("linear")
       do m = 1, m1
          do n = 1, n1
            call find_index(grid1(n,m,:), grid2(n,m,1), grid2(n,m,k2), ks, ke)
            call interp_1d_linear(grid1(n,m,ks:ke),grid2(n,m,:),data1(n,m,ks:ke),data2(n,m,:))
          enddo
       enddo
 
    case("cubic_spline")
       do m = 1, m1
          do n = 1, n1
            call find_index(grid1(n,m,:), grid2(n,m,1), grid2(n,m,k2), ks, ke)
            call interp_1d_cubic_spline(grid1(n,m,ks:ke),grid2(n,m,:), data1(n,m,ks:ke),data2(n,m,:), y1, y2)
          enddo
       enddo
 
    case default
       call mpp_error(fatal,"axis_utils: interp_method should be linear or cubic_spline")
    end select
 
    return
 
  end subroutine interp_1d_3d_
 
 
  !#####################################################################
  subroutine find_index_(grid1, xs, xe, ks, ke)
    real(kind=fms_au_kind_), dimension(:), intent(in)  :: grid1
    real(kind=fms_au_kind_),               intent(in)  :: xs, xe
    integer,                               intent(out) :: ks, ke
 
    integer :: k, nk
 
    nk = size(grid1(:))
 
    ks = 0; ke = 0
    do k = 1, nk-1
       if(grid1(k) <= xs .and. grid1(k+1) > xs ) then
          ks = k
          exit
       endif
    enddo
 
    do k = nk, 2, -1
       if(grid1(k) >= xe .and. grid1(k-1) < xe ) then
          ke = k
          exit
       endif
    enddo
 
    if(ks == 0 ) call mpp_error(fatal,' xs locate outside of grid1')
    if(ke == 0 ) call mpp_error(fatal,' xe locate outside of grid1')
 
  end subroutine find_index_
 !> @}
 ! close documentation grouping