subroutine ohcparam2(ibasin,ityr,itmon,itday,tlat,tlon,
     +                     rsst,rmiss,pohc,ierr)
c      
c-----------------------------------------------------------------------      
c     This is a routine of the SHIPS model.
c
c     This program calculates the ocean heat content (OHC) at a given 
c     point and time from a linear combination of Reynolds SST and
c     climatological OHC.  Coefficients were statistically derived
c     using a simple linear regression model.
c
c
c     Input:  ibasin - 1=Atlantic, 2=EPac, 3=WPac
c             ityr   - 4-digit year
c             itmon  - month (1-12)
c             itday  - day (1-31)
c             tlat   - TC latitude (deg N)
c             tlon   - TC longitude (deg E, ranging from 0-360E )
c             rsst   - Reynolds SST value (deg C)
c
c     Output: 
c             pohc  - Parameterized OHC value
c             ierr - error flag =0 for normal completion
c                               =1 invalid ibasin
c                               =2 rsst = rmiss (missing value)
c                               =3 ohcclim.f failed
c                               =4 rsstclim.f failed  
c
c     Uses:  ohcclim.f
c            rsstclim.f
c            jday.f
c
c            ohcclim2.atl	! 1995-2008 (source: M. Mainelli, J. Brewster)
c            ohcclim2.pac	! 2000-2008 (source: J. Brewster)
c            ohcclim2.wpn	! 2005-2010 (source: NCODA / B. Sampson)
c            rsstclim2.dat	! 1995-2009 (global pack on pasiphae:/data/andrea/rsst
c
c     11 Oct 2011 - Modified by A. Schumacher to include W. Pacific basin:
c                   included ibasin to passed variables, added W. Pac
c                   means/stdev/coefficients (from NCODA data)
c
c     24 May 2013 - Modified by M. DeMaria to separate rsstclim routine
c                   and replace ohcclim.f routine with global ncodaclim routine
c
c     24 Mar 2017 - Modified by K. Musgrave to add SH and IO (uses WP files, coef)
c    
c----------------------------------------------------------------------- 
c
c     Atlantic arrays
      dimension    meana(3)
      dimension    stdeva(3)
      dimension    coefa(2)
c     E Pacific arrays
      dimension    meane(3)
      dimension    stdeve(3)
      dimension    coefe(2)
c     W Pacific arrays
      dimension    meanw(3)
      dimension    stdevw(3)
      dimension    coefw(2)
c            
      character*3  bas
c
c     Atlantic means, standard deviations and regression coefficients.
c     Model explains  82.4% of variance in developmental dataset.
c                     dSST    COHC    ROHC
      data meana  /  0.201, 43.612, 44.016 /
      data stdeva /  0.482, 28.744, 30.743 /
      data coefa  /  0.215,  0.936         /
c
c     East Pac means, standard deviations and regression coefficients.
c     Model explains  60.9% of variance in developmental dataset.
c                     dSST    COHC    ROHC
      data meane  /  0.126, 26.437, 24.455 /
      data stdeve /  0.587, 15.564, 18.246 /
      data coefe  /  0.302,  0.777         /
c
c     West Pac means, standard deviations and regression coefficients.
c     Model explains  81.3% of variance in developmental dataset.
c                     dSST    COHC    ROHC
      data meanw  /  0.167, 65.212, 63.693 /
      data stdevw /  0.568, 32.659, 35.998 /
      data coefw  /  0.195,  0.933         /
c      
      save meana,stdeva,coefa,
     +     meane,stdeve,coefe
     +     meanw,stdevw,coefw
c
      ierr = 0
      pohc = rmiss
c
c     If RSST is 26.0 C or less, OHC = 0
      if (rsst .le. 26.0) then
         pohc = 0.0
	 return
      endif      	 
c            
      if (ibasin .eq. 1) then
         bas = 'atl'
      elseif (ibasin .eq. 2) then
         bas = 'pac'
      elseif (ibasin .eq. 3) then
         bas = 'wpn'
      elseif (ibasin .eq. 4) then
         bas = 'wpn'
      elseif (ibasin .eq. 5) then
         bas = 'wpn'
      else
         ierr = 1
	 return
      endif
c
      if (rsst .eq. rmiss) then
         ierr = 2
      endif	 
c        
c     call ohcclim(bas,ityr,itmon,itday,tlat,tlon,rmiss,
c    +                 cohc,ierrc)  
      call ncodaclim(tlat,tlon,ityr,itmon,itday,cohc,
     +               d20v,d26v,sstv,ierrc)
c
      if ((ierrc .gt. 0) .or. (cohc .eq. rmiss)) then
         ierr = 30+ierrc
	 return
      endif

      call rsstclim(ityr,itmon,itday,tlat,tlon,rmiss,
     +                 crsst,ierrs)
      if ((ierrs .gt. 0) .or. (crsst .eq. rmiss)) then
         ierr = 40+ierrs
	 return
      endif    

      if (bas .eq. 'atl') then      
c        Compute normalized values of Reynolds SST and clim OHC      
         rrsst = ((rsst-crsst)-meana(1))/stdeva(1)
         rcohc = (cohc-meana(2))/stdeva(2)
c
c        Calculate normalized value of OHC from simple linear regression model  
         pohc  = coefa(1)*rrsst + coefa(2)*rcohc
c
c        Convert normalized OHC to actual OHC      
         pohc  = (pohc*stdeva(3)) + meana(3)
      elseif (bas .eq. 'pac') then
c        Compute normalized values of Reynolds SST and clim OHC      
         rrsst = ((rsst-crsst)-meane(1))/stdeve(1)
         rcohc = (cohc-meane(2))/stdeve(2)
c
c        Calculate normalized value of OHC from simple linear regression model  
         pohc  = coefe(1)*rrsst + coefe(2)*rcohc
c
c        Convert normalized OHC to actual OHC      
         pohc  = (pohc*stdeve(3)) + meane(3)
      elseif (bas .eq. 'wpn') then
c        Compute normalized values of Reynolds SST and clim OHC      
         rrsst = ((rsst-crsst)-meanw(1))/stdevw(1)
         rcohc = (cohc-meanw(2))/stdevw(2)
c
c        Calculate normalized value of OHC from simple linear regression model  
         pohc  = coefw(1)*rrsst + coefw(2)*rcohc
c
c        Convert normalized OHC to actual OHC      
         pohc  = (pohc*stdevw(3)) + meanw(3)
      endif	          
c    
c     Since OHC > 0 by definition, set all negative predicted OHC to 0  
      if (pohc .lt. 0) pohc = 0.0       
c         
      return
      
      end