! $Id: benmod.f90 2022-11-02 09:51:28Z ychen $

!****s* fsi/benmod *
!
! NAME
!   benmod  - calculates bending angle for a given
!             latitude and impact parameter based on CIRA+Q model
!
! SYNOPSIS
!    call benmod(imonth, glat, a, b, b1, amin, bmax)
!
! DESCRIPTION
!   This subroutine calculates bending angle for a given
!   latitude and impact parameter based on CIRA+Q model.
!   Interpolation of latitude is linear;
!   interpolation of impact parameter is a cubic spline.
!   The BENMOD input file is used implicitly.
!
! INPUTS
!   INTEGER      :: imonth            month
!   REAL(wp)     :: glat              latitude
!   REAL(wp)     :: a                 input bending angle
!
! OUTPUT
!   REAL(wp)     :: b                 output bending angle
!   REAL(wp)     :: b1                output bending angle derivative
!   REAL(wp)     :: amin              min impact height
!   REAL(wp)     :: bmax              max bending angle
!
! AUTHOR
!   S.V.Sokolovskiy, UCAR
!   Update: Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE benmod(imonth, glat, a, b, b1, amin, bmax)

      USE typesizes, ONLY: wp => EightByteReal
      USE load_ciraq_model, only: am, bm
      USE ropp_pp_spline

      IMPLICIT NONE
      INTEGER,  INTENT(in)     :: imonth           ! month
      REAL(wp), INTENT(in)     :: glat             ! latitude
      REAL(wp), INTENT(in)     :: a                ! input bending angle
      REAL(wp), INTENT(out)    :: b                ! output bending angle
      REAL(wp), INTENT(out)    :: b1               ! output bending angle derivative
      REAL(wp), INTENT(out)    :: amin             ! min impact height
      REAL(wp), INTENT(out)    :: bmax             ! max bending angle


      REAL(wp), DIMENSION(:), ALLOCATABLE :: amd, bmd, d2

      INTEGER  :: i, k, n
      REAL(wp) :: gl, dgl1, dgl

      gl    = glat + 90.0_wp
      k     = floor(gl / 10._wp) + 1
      dgl1  = (gl - 10.0_wp * (k - 1) ) / 10.0_wp
      dgl   = 1.0_wp - dgl1

      n = size(am, 3)
      ALLOCATE(amd(n))
      ALLOCATE(bmd(n))
      ALLOCATE(d2(n))

      DO i = 1, n
        amd(i) = am(imonth, k, i) * dgl + am(imonth, k + 1, i) * dgl1
        bmd(i) = bm(imonth, k, i) * dgl + bm(imonth, k + 1, i) * dgl1
      enddo

      !CALL splinx1(amd, bmd, a, b, b1)

      CALL ropp_pp_init_spline(amd, bmd, d2)
      CALL ropp_pp_interpol_spline(amd, bmd, d2, a, b, b1)

      b = exp(b)
      b1 = b * b1

      amin = amd(1)
      bmax = exp(bmd(1))

      DEALLOCATE(amd)
      DEALLOCATE(bmd)
      DEALLOCATE(d2)

 END SUBROUTINE benmod