suydam_kgs.f90

Go to the documentation of this file.

subroutine suydam_kgs(t, t0, cp, cq, ipsi)
!-----------------------------------------------------------------------
! subroutine to evaluate the p and q coefficient, used in the suydam
! routine. 
! note : the routine asumes that only values on grid points are
! requested. no interpolation is done!
!-----------------------------------------------------------------------      

  use itm_types
  use mod_parameter
  use mod_suydam

  implicit none

  integer(itm_i4), intent(in) :: ipsi
  real(r8), intent(in) :: t, t0
  real(r8), intent(out) :: cp, cq

  integer(itm_i4) :: jm, ijm
  real(r8) :: dt, tm, frac
  real(r8) :: cp0d, cp1d, cp2d, cq0d, cq1d

  dt = twopi / real(ncpq - 1)
  tm = mod(mod(t, twopi) + twopi, twopi)

!------ ncpq = 2 * nchi - 1 (ias=0) or ncpq = nchi + 1 (ias = 1)      
  jm = int((ncpq - 1) * tm / twopi) + 1
  ijm = (ipsi - 1) * ncpq + jm
  frac = (tm - (jm - 1) * dt) / dt

  cp0d = cp0(ijm) + (cp0(ijm + 1) - cp0(ijm)) * frac
  cp1d = cp1(ijm) + (cp1(ijm + 1) - cp1(ijm)) * frac
  cp2d = cp2(ijm) + (cp2(ijm + 1) - cp2(ijm)) * frac
  cq0d = cq0(ijm) + (cq0(ijm + 1) - cq0(ijm)) * frac 
  cq1d = cq1(ijm) + (cq1(ijm + 1) - cq1(ijm)) * frac           

  cp = cp0d + cp1d * (t - t0) + cp2d * (t - t0)**2 
  cq = cq0d + cq1d * (t - t0)

  return

end subroutine suydam_kgs