axis_utils/axis_utils2.F90

   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 module axis_utils2_mod
  28   use mpp_mod,      only: mpp_error, FATAL, stdout
  29   use fms_mod,      only: lowercase, uppercase, string_array_index, fms_error_handler
  30   use fms2_io_mod,  only: FmsNetcdfDomainFile_t, variable_att_exists, FmsNetcdfFile_t, &
  31                           get_variable_num_dimensions, get_variable_attribute,  &
  32                           get_variable_size, read_data, variable_exists
  33   use platform_mod, only: r4_kind, r8_kind
  34
  35   implicit none
  36
  37   public get_axis_cart, get_axis_modulo, lon_in_range, &
  38          tranlon, frac_index, nearest_index, interp_1d, get_axis_modulo_times, axis_edges
  39
  40   private
  41
  42   integer, parameter :: maxatts = 100
  43   real(r8_kind), parameter    :: epsln = 1.e-10_r8_kind
  44   real(r8_kind), parameter    :: fp5 = 0.5_r8_kind, f360 = 360.0_r8_kind
  45
  46 !> @}
  47 ! Include variable "version" to be written to log file.
  48 #include<file_version.h>
  49
  50   !> Perform 1D interpolation between grids.
  51   !!
  52   !> Data and grids can have 1, 2, or 3 dimensions.
  53   !! @param grid1 grid for data1
  54   !! @param grid2 grid for data2
  55   !! @param data1 Data to interpolate
  56   !! @param [inout] data2 Interpolated data
  57   !! @param method Either "linear" or "cubic_spline" interpolation method, default="linear"
  58   !! @ingroup axis_utils2_mod
  59
  60   interface axis_edges
  61     module procedure axis_edges_r4, axis_edges_r8
  62   end interface axis_edges
  63
  64   interface lon_in_range
  65     module procedure lon_in_range_r4, lon_in_range_r8
  66   end interface lon_in_range
  67
  68   interface frac_index
  69     module procedure frac_index_r4, frac_index_r8
  70   end interface frac_index
  71
  72   interface nearest_index
  73       module procedure nearest_index_r4, nearest_index_r8
  74   end interface nearest_index
  75
  76   interface tranlon
  77       module procedure tranlon_r4, tranlon_r8
  78   end interface tranlon
  79
  80   interface interp_1d_linear
  81       module procedure interp_1d_linear_r4, interp_1d_linear_r8
  82   end interface interp_1d_linear
  83
  84   interface interp_1d_cubic_spline
  85         module procedure interp_1d_cubic_spline_r4, interp_1d_cubic_spline_r8
  86   end interface interp_1d_cubic_spline
  87
  88   interface interp_1d
  89      module procedure interp_1d_1d_r4, interp_1d_1d_r8
  90      module procedure interp_1d_2d_r4, interp_1d_2d_r8
  91      module procedure interp_1d_3d_r4, interp_1d_3d_r8
  92   end interface interp_1d
  93
  94   interface find_index
  95       module procedure find_index_r4, find_index_r8
  96   end interface find_index
  97
  98 !> @addtogroup axis_utils2_mod
  99 !> @{
 100
 101 contains
 102
 103   !> @brief Returns X,Y,Z or T cartesian attribute
 104   subroutine get_axis_cart(fileobj, axisname, cart)
 105     type(FmsNetcdfFile_t), intent(in) :: fileobj !< file object to read from
 106     character(len=*), intent(in) :: axisname !< name of axis to retrieve
 107     character(len=1), intent(out) :: cart !< Returned attribute axis
 108
 109     character(len=1) :: axis_cart
 110     character(len=16), dimension(2) :: lon_names, lat_names
 111     character(len=16), dimension(3) :: z_names
 112     character(len=16), dimension(2) :: t_names
 113     character(len=16), dimension(3) :: lon_units, lat_units
 114     character(len=8) , dimension(4) :: z_units
 115     character(len=3) , dimension(6) :: t_units
 116     character(len=32) :: name
 117     integer :: i
 118
 119     lon_names = (/'lon','x  '/)
 120     lat_names = (/'lat','y  '/)
 121     z_names = (/'depth ','height','z     '/)
 122     t_names = (/'time','t   '/)
 123     lon_units = (/'degrees_e   ', 'degrees_east', 'degreese    '/)
 124     lat_units = (/'degrees_n    ', 'degrees_north', 'degreesn     '/)
 125     z_units = (/'cm ','m  ','pa ','hpa'/)
 126     t_units = (/'sec', 'min','hou','day','mon','yea'/)
 127
 128     cart = "N"
 129     if (variable_exists(fileobj, axisname)) then
 130       if (variable_att_exists(fileobj, axisname, "cartesian_axis")) then
 131         call get_variable_attribute(fileobj, axisname, "cartesian_axis", cart(1:1))
 132       elseif (variable_att_exists(fileobj, axisname, "axis")) then
 133         call get_variable_attribute(fileobj, axisname, "axis", cart(1:1))
 134       endif
 135       axis_cart = uppercase(cart)
 136       if (axis_cart .eq. 'X' .or. axis_cart .eq. 'Y' .or. axis_cart .eq. 'Z' &
 137           .or. axis_cart .eq. 'T') then
 138         cart = axis_cart
 139         return
 140       endif
 141     endif
 142
 143     if (cart /= 'X' .and. cart /= 'Y' .and. cart /= 'Z' .and. cart /= 'T') then
 144        name = lowercase(axisname)
 145        do i=1,size(lon_names(:))
 146           if (trim(name(1:3)) == trim(lon_names(i))) cart = 'X'
 147        enddo
 148        do i=1,size(lat_names(:))
 149           if (trim(name(1:3)) == trim(lat_names(i))) cart = 'Y'
 150        enddo
 151        do i=1,size(z_names(:))
 152           if (trim(name) == trim(z_names(i))) cart = 'Z'
 153        enddo
 154        do i=1,size(t_names(:))
 155           if (trim(name) == t_names(i)) cart = 'T'
 156        enddo
 157     end if
 158
 159     if (cart /= 'X' .and. cart /= 'Y' .and. cart /= 'Z' .and. cart /= 'T') then
 160        name = lowercase(axisname)
 161        do i=1,size(lon_units(:))
 162           if (trim(name) == trim(lon_units(i))) cart = 'X'
 163        enddo
 164        do i=1,size(lat_units(:))
 165           if (trim(name) == trim(lat_units(i))) cart = 'Y'
 166        enddo
 167        do i=1,size(z_units(:))
 168           if (trim(name) == trim(z_units(i))) cart = 'Z'
 169        enddo
 170        do i=1,size(t_units(:))
 171           if (name(1:3) == trim(t_units(i))) cart = 'T'
 172        enddo
 173     end if
 174   end subroutine get_axis_cart
 175
 176   !> @brief Checks if 'modulo' variable exists for a given axis.
 177   !!
 178   !> @return true if modulo variable exists in fileobj for the given axis name.
 179   function get_axis_modulo(fileobj, axisname)
 180     type(FmsNetcdfFile_t), intent(in) :: fileobj
 181     character(len=*), intent(in) :: axisname
 182     logical :: get_axis_modulo
 183
 184     get_axis_modulo = variable_att_exists(fileobj, axisname, "modulo")
 185   end function get_axis_modulo
 186
 187   !> @return true if modulo_beg and modulo_end exist in fileobj with the given
 188   !! axis, and returns their values in tbeg and tend.
 189   function get_axis_modulo_times(fileobj, axisname, tbeg, tend)
 190     type(FmsNetcdfFile_t), intent(in) :: fileobj
 191     character(len=*), intent(in) :: axisname
 192     character(len=*), intent(out) :: tbeg, tend
 193     logical :: get_axis_modulo_times
 194     logical :: found_tbeg, found_tend
 195
 196     found_tbeg = variable_att_exists(fileobj, axisname, "modulo_beg")
 197     found_tend = variable_att_exists(fileobj, axisname, "modulo_end")
 198
 199     if (found_tbeg .and. .not. found_tend) then
 200       call mpp_error(FATAL,'error in get: Found modulo_beg but not modulo_end')
 201     endif
 202     if (.not. found_tbeg .and. found_tend) then
 203       call mpp_error(FATAL,'error in get: Found modulo_end but not modulo_beg')
 204     endif
 205
 206     if (found_tbeg) then
 207       call get_variable_attribute(fileobj, axisname, "modulo_beg", tbeg)
 208       call get_variable_attribute(fileobj, axisname, "modulo_end", tend)
 209     else
 210       tbeg = ""
 211       tend = ""
 212     endif
 213     get_axis_modulo_times = found_tbeg
 214   end function get_axis_modulo_times
 215
 216 #include "axis_utils2_r4.fh"
 217 #include "axis_utils2_r8.fh"
 218
 219 end module axis_utils2_mod
 220 !> @}
 221 ! close documentation grouping