module glc_elevclass_mod

  !---------------------------------------------------------------------
  !
  ! Purpose:
  !
  ! This module contains data and routines for operating on GLC elevation classes.

#include "shr_assert.h"
  use shr_kind_mod, only : r8 => shr_kind_r8
  use shr_sys_mod
  use seq_comm_mct, only : logunit
  use shr_log_mod, only : errMsg => shr_log_errMsg

  implicit none
  save
  private

  !--------------------------------------------------------------------------
  ! Public interfaces
  !--------------------------------------------------------------------------

  public :: glc_elevclass_init            ! initialize GLC elevation class data
  public :: glc_elevclass_clean           ! deallocate memory allocated here
  public :: glc_get_num_elevation_classes ! get the number of elevation classes
  public :: glc_get_elevation_class       ! get the elevation class index for a given elevation
  public :: glc_get_elevclass_bounds      ! get the boundaries of all elevation classes
  public :: glc_mean_elevation_virtual    ! get the mean elevation of a virtual elevation class
  public :: glc_elevclass_as_string       ! returns a string corresponding to a given elevation class
  public :: glc_all_elevclass_strings     ! returns an array of strings for all elevation classes
  public :: glc_errcode_to_string         ! convert an error code into a string describing the error

  interface glc_elevclass_init
     module procedure glc_elevclass_init_default
     module procedure glc_elevclass_init_override
  end interface glc_elevclass_init


  !--------------------------------------------------------------------------
  ! Public data
  !--------------------------------------------------------------------------

  ! Possible error code values
  integer, parameter, public :: GLC_ELEVCLASS_ERR_NONE = 0      ! err_code indicating no error
  integer, parameter, public :: GLC_ELEVCLASS_ERR_UNDEFINED = 1 ! err_code indicating elevation classes have not been defined
  integer, parameter, public :: GLC_ELEVCLASS_ERR_TOO_LOW = 2   ! err_code indicating topo below lowest elevation class
  integer, parameter, public :: GLC_ELEVCLASS_ERR_TOO_HIGH = 3  ! err_code indicating topo above highest elevation class

  ! String length for glc elevation classes represented as strings
  integer, parameter, public :: GLC_ELEVCLASS_STRLEN = 2

  !--------------------------------------------------------------------------
  ! Private data
  !--------------------------------------------------------------------------

  ! number of elevation classes
  integer :: glc_nec

  ! upper elevation limit of each class (m)
  ! indexing starts at 0, with topomax(0) giving the lower elevation limit of EC 1
  real(r8), allocatable :: topomax(:)


contains

  !-----------------------------------------------------------------------
  subroutine glc_elevclass_init_default(my_glc_nec)
    !
    ! !DESCRIPTION:
    ! Initialize GLC elevation class data to default boundaries, based on given glc_nec
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    integer, intent(in) :: my_glc_nec  ! number of GLC elevation classes
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'glc_elevclass_init'
    !-----------------------------------------------------------------------

    glc_nec = my_glc_nec
    allocate(topomax(0:glc_nec))

    select case (glc_nec)
    case(0)
       ! do nothing
    case(1)
       topomax = [0._r8, 10000._r8]
    case(3)
       topomax = [0._r8,  1000._r8,  2000._r8, 10000._r8]
    case(5)
       topomax = [0._r8,   500._r8,  1000._r8,  1500._r8, 2000._r8, 10000._r8]
    case(10)
       topomax = [0._r8,   200._r8,   400._r8,   700._r8,  1000._r8,  1300._r8,  &
            1600._r8,  2000._r8,  2500._r8,  3000._r8, 10000._r8]
    case(36)
       topomax = [  0._r8,   200._r8,   400._r8,   600._r8,   800._r8,  &
            1000._r8,  1200._r8,  1400._r8,  1600._r8,  1800._r8,  &
            2000._r8,  2200._r8,  2400._r8,  2600._r8,  2800._r8,  &
            3000._r8,  3200._r8,  3400._r8,  3600._r8,  3800._r8,  &
            4000._r8,  4200._r8,  4400._r8,  4600._r8,  4800._r8,  &
            5000._r8,  5200._r8,  5400._r8,  5600._r8,  5800._r8,  &
            6000._r8,  6200._r8,  6400._r8,  6600._r8,  6800._r8,  &
            7000._r8, 10000._r8]
    case default
       write(logunit,*) subname,' ERROR: unknown glc_nec: ', glc_nec
       call shr_sys_abort(subname//' ERROR: unknown glc_nec')
    end select

  end subroutine glc_elevclass_init_default

  !-----------------------------------------------------------------------
  subroutine glc_elevclass_init_override(my_glc_nec, my_topomax)
    !
    ! !DESCRIPTION:
    ! Initialize GLC elevation class data to the given elevation class boundaries.
    !
    ! The input, my_topomax, should have (my_glc_nec + 1) elements.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    integer, intent(in) :: my_glc_nec      ! number of GLC elevation classes
    real(r8), intent(in) :: my_topomax(0:) ! elevation class boundaries (m)
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'glc_elevclass_init_override'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(my_topomax) == (/my_glc_nec/)), __FILE__, __LINE__)

    glc_nec = my_glc_nec
    allocate(topomax(0:glc_nec))
    topomax = my_topomax

  end subroutine glc_elevclass_init_override

  !-----------------------------------------------------------------------
  subroutine glc_elevclass_clean()
    !
    ! !DESCRIPTION:
    ! Deallocate memory allocated in this module

    character(len=*), parameter :: subname = 'glc_elevclass_clean'
    !-----------------------------------------------------------------------

    if (allocated(topomax)) then
       deallocate(topomax)
    end if
    glc_nec = 0

  end subroutine glc_elevclass_clean

  !-----------------------------------------------------------------------
  function glc_get_num_elevation_classes() result(num_elevation_classes)
    !
    ! !DESCRIPTION:
    ! Get the number of GLC elevation classes
    !
    ! !ARGUMENTS:
    integer :: num_elevation_classes  ! function result
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'glc_get_num_elevation_classes'
    !-----------------------------------------------------------------------

    num_elevation_classes = glc_nec

  end function glc_get_num_elevation_classes

  !-----------------------------------------------------------------------
  subroutine glc_get_elevation_class(topo, elevation_class, err_code)
    !
    ! !DESCRIPTION:
    ! Get the elevation class index associated with a given topographic height.
    !
    ! The returned elevation_class will be between 1 and num_elevation_classes, if this
    ! topographic height is contained in an elevation class. In this case, err_code will
    ! be GLC_ELEVCLASS_ERR_NONE (no error).
    !
    ! If there are no elevation classes defined, the returned value will be 0, and
    ! err_code will be GLC_ELEVCLASS_ERR_UNDEFINED
    !
    ! If this topographic height is below the lowest elevation class, the returned value
    ! will be 1, and err_code will be GLC_ELEVCLASS_ERR_TOO_LOW.
    !
    ! If this topographic height is above the highest elevation class, the returned value
    ! will be (num_elevation_classes), and err_code will be GLC_ELEVCLASS_ERR_TOO_HIGH.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8), intent(in) :: topo ! topographic height (m)
    integer, intent(out) :: elevation_class  ! elevation class index
    integer, intent(out) :: err_code ! error code (see above for possible codes)
    !
    ! !LOCAL VARIABLES:
    integer :: ec  ! temporary elevation class

    character(len=*), parameter :: subname = 'glc_get_elevation_class'
    !-----------------------------------------------------------------------

    if (glc_nec < 1) then
       elevation_class = 0
       err_code = GLC_ELEVCLASS_ERR_UNDEFINED
    else if (topo < topomax(0)) then
       elevation_class = 1
       err_code = GLC_ELEVCLASS_ERR_TOO_LOW
    else if (topo >= topomax(glc_nec)) then
       elevation_class = glc_nec
       err_code = GLC_ELEVCLASS_ERR_TOO_HIGH
    else
       err_code = GLC_ELEVCLASS_ERR_NONE
       elevation_class = 0
       do ec = 1, glc_nec
          if (topo >= topomax(ec - 1) .and. topo < topomax(ec)) then
             elevation_class = ec
             exit
          end if
       end do

       SHR_ASSERT(elevation_class > 0, subname//' elevation class was not assigned')
    end if

  end subroutine glc_get_elevation_class

  !-----------------------------------------------------------------------
  function glc_get_elevclass_bounds() result(elevclass_bounds)
    !
    ! !DESCRIPTION:
    ! Get the boundaries of all elevation classes.
    !
    ! This returns an array of size glc_nec+1, since it contains both the lower and upper
    ! bounds of each elevation class.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8) :: elevclass_bounds(0:glc_nec)  ! function result
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'glc_get_elevclass_bounds'
    !-----------------------------------------------------------------------

    elevclass_bounds(:) = topomax(:)

  end function glc_get_elevclass_bounds


  !-----------------------------------------------------------------------
  function glc_mean_elevation_virtual(elevation_class) result(mean_elevation)
    !
    ! !DESCRIPTION:
    ! Returns the mean elevation of a virtual elevation class
    !
    ! !ARGUMENTS:
    real(r8) :: mean_elevation  ! function result
    integer, intent(in) :: elevation_class
    !
    ! !LOCAL VARIABLES:
    integer :: resulting_elevation_class
    integer :: err_code

    character(len=*), parameter :: subname = 'glc_mean_elevation_virtual'
    !-----------------------------------------------------------------------

    if (elevation_class == 0) then
       ! Bare land "elevation class"
       mean_elevation = 0._r8
    else
       if (elevation_class < glc_nec) then
          ! Normal case
          mean_elevation = (topomax(elevation_class - 1) + topomax(elevation_class)) / 2._r8
       else if (elevation_class == glc_nec) then
          ! In the top elevation class; in this case, assignment of a "mean" elevation is
          ! somewhat arbitrary (because we expect the upper bound of the top elevation
          ! class to be very high).

          if (glc_nec > 1) then
             mean_elevation = 2._r8 * topomax(elevation_class - 1) - topomax(elevation_class - 2)
          else
             ! entirely arbitrary
             mean_elevation = 1000._r8
          end if
       else
          write(logunit,*) subname,' ERROR: elevation class out of bounds: ', elevation_class
          call shr_sys_abort(subname // ' ERROR: elevation class out of bounds')
       end if
    end if

    ! Ensure that the resulting elevation is within the given elevation class
    if (elevation_class > 0) then
       call glc_get_elevation_class(mean_elevation, resulting_elevation_class, err_code)
       if (err_code /= GLC_ELEVCLASS_ERR_NONE) then
          write(logunit,*) subname, ' ERROR: generated elevation that results in an error'
          write(logunit,*) 'when trying to determine the resulting elevation class'
          write(logunit,*) glc_errcode_to_string(err_code)
          write(logunit,*) 'elevation_class, mean_elevation = ', elevation_class, mean_elevation
          call shr_sys_abort(subname // ' ERROR: generated elevation that results in an error')
       else if (resulting_elevation_class /= elevation_class) then
          write(logunit,*) subname, ' ERROR: generated elevation outside the given elevation class'
          write(logunit,*) 'elevation_class, mean_elevation, resulting_elevation_class = ', &
               elevation_class, mean_elevation, resulting_elevation_class
          call shr_sys_abort(subname // ' ERROR: generated elevation outside the given elevation class')
       end if
    end if

  end function glc_mean_elevation_virtual


  !-----------------------------------------------------------------------
  function glc_elevclass_as_string(elevation_class) result(ec_string)
    !
    ! !DESCRIPTION:
    ! Returns a string corresponding to a given elevation class.
    !
    ! This string can be used as a suffix for fields in MCT attribute vectors.
    !
    ! ! NOTE(wjs, 2015-01-19) This function doesn't fully belong in this module, since it
    ! doesn't refer to the data stored in this module. However, I can't think of a more
    ! appropriate place for it.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    character(len=GLC_ELEVCLASS_STRLEN) :: ec_string  ! function result
    integer, intent(in) :: elevation_class
    !
    ! !LOCAL VARIABLES:
    character(len=16) :: format_string

    character(len=*), parameter :: subname = 'glc_elevclass_as_string'
    !-----------------------------------------------------------------------

    ! e.g., for GLC_ELEVCLASS_STRLEN = 2, format_string will be '(i2.2)'
    write(format_string,'(a,i0,a,i0,a)') '(i', GLC_ELEVCLASS_STRLEN, '.', GLC_ELEVCLASS_STRLEN, ')'

    write(ec_string,trim(format_string)) elevation_class
  end function glc_elevclass_as_string

  !-----------------------------------------------------------------------
  function glc_all_elevclass_strings(include_zero) result(ec_strings)
    !
    ! !DESCRIPTION:
    ! Returns an array of strings corresponding to all elevation classes from 1 to glc_nec
    !
    ! If include_zero is present and true, then includes elevation class 0 - so goes from
    ! 0 to glc_nec
    !
    ! These strings can be used as suffixes for fields in MCT attribute vectors.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    character(len=GLC_ELEVCLASS_STRLEN), allocatable :: ec_strings(:)  ! function result
    logical, intent(in), optional :: include_zero   ! if present and true, include elevation class 0 (default is false)
    !
    ! !LOCAL VARIABLES:
    logical :: l_include_zero  ! local version of optional include_zero argument
    integer :: lower_bound
    integer :: i

    character(len=*), parameter :: subname = 'glc_all_elevclass_strings'
    !-----------------------------------------------------------------------

    if (present(include_zero)) then
       l_include_zero = include_zero
    else
       l_include_zero = .false.
    end if

    if (l_include_zero) then
       lower_bound = 0
    else
       lower_bound = 1
    end if

    allocate(ec_strings(lower_bound:glc_nec))
    do i = lower_bound, glc_nec
       ec_strings(i) = glc_elevclass_as_string(i)
    end do

  end function glc_all_elevclass_strings


  !-----------------------------------------------------------------------
  function glc_errcode_to_string(err_code) result(err_string)
    !
    ! !DESCRIPTION:
    !
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    character(len=256) :: err_string  ! function result
    integer, intent(in) :: err_code   ! error code (one of the GLC_ELEVCLASS_ERR* values)
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'glc_errcode_to_string'
    !-----------------------------------------------------------------------

    select case (err_code)
    case (GLC_ELEVCLASS_ERR_NONE)
       err_string = '(no error)'
    case (GLC_ELEVCLASS_ERR_UNDEFINED)
       err_string = 'Elevation classes have not yet been defined'
    case (GLC_ELEVCLASS_ERR_TOO_LOW)
       err_string = 'Topographic height below the lower bound of the lowest elevation class'
    case (GLC_ELEVCLASS_ERR_TOO_HIGH)
       err_string = 'Topographic height above the upper bound of the highest elevation class'
    case default
       err_string = 'UNKNOWN ERROR'
    end select

  end function glc_errcode_to_string


end module glc_elevclass_mod
