subroutine mpicalo(sst,rmpi,ibasin)
c     This routine calculates the maximum potential intensity (kt)
c     from the sst (C) using empirical relationships.
c
c     This version was modifed March 2006 to
c     use the original formulas from DK 1994 and
c     WH 1997, rather than the formulas adjusted
c     for the average storm translational speed.
c
c     Modified June 2013 to add IO and SH functions
c
c     Input: 
c       sst:    SST in deg C
c       tspeed: Storm translational speed in knots
c       ibasin: Basin indicator
c               ibasin=1 for Atlantic
c               ibasin=2 for East Pacific
c               ibasin=3 for West Pacific
c               ibasin=4 for the Indian Ocean
c               ibasin=5 for the Southern hemisphere
c      
c     Output:
c       rmpi:  Maximum potential intensity (kt)
c              Note: rmpi is set to 999.9 for missing SST
c
c     Check for illegal sst values 
      if (sst .gt. 35.0 .or. sst .lt. 0.0) then
         rmpi=999.9
         return
      endif
c
      if (ibasin .eq. 1) then
c        Atlantic function (DeMaria and Kaplan 1994)
c        vcold = 34.2
         vcold = 28.2
         vadd  = 55.8
         a     = 0.1813
         tmax  = 30.00
         vcap  = 165.0
c
         rmpi = vcold + vadd*exp(-a*(tmax-sst))
         rmpi = rmpi*1.944
      elseif (ibasin .eq. 2) then
c        East Pacific function (Whitney and Hobgood 1997)
         a = -79.2
         b = 5.362
         c = 0.0
         vcap  = 165.0
C
         tmin = 20.0
         sstt = sst
         if (sstt .lt. tmin) sstt=tmin
C
         rmpi = a + b*sstt + c
         rmpi = rmpi*1.944
      elseif (ibasin .eq. 3) then
C        West Pacific function
         vcold = 19.7
         vadd  = 88.0
         a     = .1909
         tmax  = 30.00
         vcap  = 165.0
C
         rmpi  = vcold + vadd*exp(-a*(tmax-sst))
         rmpi  = rmpi*1.944
      elseif (ibasin .eq. 4) then
C        N. Indian Ocean function
         vcold = 25.0
         vadd  = 70.0
         a     = .35
         tmax  = 30.00
         vcap  = 165.0
C
         rmpi  = vcold + vadd*exp(-a*(tmax-sst))
         rmpi  = rmpi*1.944
      elseif (ibasin .eq. 5) then
C        S. Hemisphere function (same as W. Pac)
         vcold = 19.7
         vadd  = 88.0
         a     = .1909
         tmax  = 30.00
         vcap  = 165.0
C
         rmpi  = vcold + vadd*exp(-a*(tmax-sst))
         rmpi  = rmpi*1.944
      else
         rmpi = 999.9
         return
      endif
C
      if (rmpi .gt. vcap) rmpi=vcap
c
      return
      end