c     This group of subroutines is for interpolating the
c     SST and Mainelli or NCODA Ocean Heat Content files.
c
c     Version 2.3.0
c
c     Modified Apr 2009: ohcparam and related routines added (MD)
c     Modified Dec 2011: ohcparam, ohclim, rsstclim, jday, bassel  routines removed 
c                        and West Pacific option added to rysst and ryohc (MD)
c     Modified Oct 2013: Arrays larger, basin input omitted. 
c                        SST file is assumed to be global version.
c     Modified Feb 2018: Add DSST and DSTA processing.
c                        Modify ryohc to use new eohc files.
c                        Version number added (v2.2.0)
c     Modified May 2018: Added nahcon subroutine to read all variables from the new    
c                        ocean heat content files (eohc) (MD)                   
c     Modified Dec 2018: The nahcon and nrdata OHC subroutines were removed
c                        for v2.3.0 so they can remain consistent with
c                        the SHIPS developmental code. Also, ryohc is
c                        for the old OHC files and is no longer used 
c                        in ilsin.f. (MD)
c     Modified Feb 2023: NCO specifications (KM)
c             
c     Routines included: rysst
c                        rydsst
c                        ryohc                        
c                        getohc
c
c     Additional routines needed: 
c                        getvar
c                        getgsst
c                        jday
c                        llintsp
c
      subroutine rysst(rlat,rlon,sst,iyr,imon,iday,itime,ifirst,ierr)
c
c     This routine finds the sst from global weekly SST analyses.
c
c     The weekly sst data is assumed to be in an ASCII file called
c     GSSTYYYYMMDD.dat with global fields. 
c
c     Input:
c        rlat:   The latitude in deg N
c        rlon:   The longitude in deg E (deg W negative can also be used)
c
c     Output:
c        sst:    The interpolated SST value
c        iyr:    The year on the SST data file
c        imon:   The month on the SST data file
c        iday:   The day on the SST data file
c        itime:  The time on the SST data file
c        ifirst: Flag for calling getgsst (skipped if ifirst=0)
c        ierr:   =0 for normal return or
c                =3 if there is an error reading the data input file
c                =4 if SST is not found on data input file
c
c     *******************************************************************
c
      parameter (mxlon=1801,mxlat=901)
      dimension sstll(mxlon,mxlat)
c
      character *1 type
c
      character *21 fnsst
c
      save sstll
      save glon1,glon2,dlong,glat1,glat2,dlatg 
      save nlong,nlatg
c
c     Set default value
      sst = 999.9
c
c     Specify unit number for input files
      lusst  = 60
c
      if (rlon .ge. 0.0) then
         rloneast = rlon
      else
         rloneast = rlon + 360.0
      endif
c
      fnsst  = 'GSSTYYYYMMDD.dat     '
c
      if (ifirst .eq. 1) then
c        Open the SST file
         open(unit=lusst,file=fnsst,status='old',form='formatted',
     +        err=920)
c
c        Get the sst array from the data file
         call getgsst(lusst,mxlon,mxlat,iyearg,imong,idayg,ifound,
     +                glon1,glon2,dlong,nlong,
     +                glat1,glat2,dlatg,nlatg,sstll)
         close(lusst)
         ifirst = 0
         if (ifound .ne. 1) then
c           Error reading the data file 
            sst = 999.9
            nday=9999
            ierr=4
            return
         endif
      endif
c
      slon = glon1
      dlon = dlong
      nlon = nlong
c
      slat = glat1
      dlat = dlatg
      nlat = nlatg
c
      iyr = iyearg - (iyearg/100)*100
      imon = imong
      iday = idayg
      itime = 0
c
c     Interpolate to find the sst
      izp  = 0
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sst,rloneast,rlat,izp,ierrll)
c
c     Normal return
      ierr=0
      close(lusst)
      return
c
  920 continue
c     Error opening the data file 
      sst = 999.9
      nday=9999
      ierr=3
      close(lusst)
      return
c
      end 
      subroutine rydsst(rlat,rlon,sst,ssta,distavg,
     +                  iyr,imon,iday,itime,ifirst,ierr)
c
c     This routine finds the sst from daily SST analyses.
c     The area average sst (ssta) is also calculated. 
c
c     The daily sst data is assumed to be in an ASCII file called
c     DSSTYYYYMMDD.dat with global fields. 
c
c     Input:
c        rlat:   The latitude in deg N
c        rlon:   The longitude in deg E (deg W negative can also be used)
c
c     Output:
c        sst:    The interpolated SST value at the storm center
c        ssta:   The SST averaged at the center and 4 points E, N, W and S
c                of the center. 
c        distavg: The distance in km for calculating ssta. 
c        iyr:    The year on the SST data file
c        imon:   The month on the SST data file
c        iday:   The day on the SST data file
c        itime:  The time on the SST data file
c        ifirst: Flag for calling getgsst (skipped if ifirst=0)
c        ierr:   =0 for normal return or
c                =3 if there is an error reading the data input file
c                =4 if SST is not found on data input file
c
c     *******************************************************************
c
      parameter (mxlon=1801,mxlat=901)
      dimension sstll(mxlon,mxlat)
c
      character *1 type
c
      character *21 fnsst
c
      save sstll
      save glon1,glon2,dlong,glat1,glat2,dlatg 
      save nlong,nlatg
c
c     Set default value
      sst = 999.9
c
c     Specify unit number for input files
      lusst  = 62
c
      if (rlon .ge. 0.0) then
         rloneast = rlon
      else
         rloneast = rlon + 360.0
      endif
c
      fnsst  = 'DSSTYYYYMMDD.dat     '
c
c     Get the dsst array from the data file
      if (ifirst .eq. 1) then
c        Open the SST file
         open(unit=lusst,file=fnsst,status='old',form='formatted',
     +        err=920)
c
         call getgsst(lusst,mxlon,mxlat,iyearg,imong,idayg,ifound,
     +                glon1,glon2,dlong,nlong,
     +                glat1,glat2,dlatg,nlatg,sstll)
         ifirst=0
         close(lusst)
         if (ifound .ne. 1) then
c           Error reading the data file 
            sst = 999.9
            nday=9999
            ierr=4
            return
         endif
      endif
c
      slon = glon1
      dlon = dlong
      nlon = nlong
c
      slat = glat1
      dlat = dlatg
      nlat = nlatg
c
      iyr = iyearg - (iyearg/100)*100
      imon = imong
      iday = idayg
      itime = 0
c
c     Interpolate to find the sst
      izp  = 0
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sst,rloneast,rlat,izp,ierrll)
c
c     Calculate area averaged dsst
      dtr = 3.14159/180.0
      dlataa = distavg/111.1
      dlonaa = dlataa/(cos(abs(rlat*dtr)))
      rlatpr = rlat     + dlataa
      rlatmr = rlat     - dlataa
      rlonpr = rloneast + dlonaa
      rlonmr = rloneast - dlonaa
c
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sstp10,rlonpr,rlat,izp,ierrll)
c
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sstm10,rlonmr,rlat,izp,ierrll)
c
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sstp01,rloneast,rlatpr,izp,ierrll)
c
      call llintsp(sstll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             sstm01,rloneast,rlatmr,izp,ierrll)
c
      ssta = 0.2*(sst + sstp01 + sstm01 + sstp10 + sstm10) 
c
c     Normal return
      ierr=0
      close(lusst)
      return
c
  920 continue
c     Error opening the data file 
      sst = 999.9
      nday=9999
      ierr=3
      close(lusst)
      return
c
      end 
      subroutine ryohc(rlat,rlon,ohc,iyear,imon,iday,itime,ifirst,ierr)
c
c     This routine finds the Ocean Heat Content (OHC) from analyses
c     on an evenly spaced lat/lon grid generated in the new eohc format. 
c     The OHC data is assumed to be in an ascii file OHCYYYYMMDD.DAT.
c     The data are interpolated to the requested location (rlat,rlon).
c     If the location is outside of the analyses domain, the OHC value is
c     set to rmiss.
c
c     Input:
c        rlat:   The latitude in deg N
c        rlon:   The longitude in deg E (deg W is negative)
c
c     Output:
c        ohc:    The interpolated OHC value
c        iyear:  The year (4-digit) on the OHC data file
c        imon:   The month on the OHC data file
c        iday:   The day on the OHC data file
c        ifirst: Flag for calling getohc (skipped if ifirst=0)
c        ierr:   =0 for normal return or
c                =3 if there is an error reading the data input file
c                =4 if OHC is not found on data input file
c
c     *******************************************************************
c
      parameter (mxlon=1801,mxlat=901)
      dimension ohcll(mxlon,mxlat)
      save ohcll
      save nlat,nlon,ierrd
      save rlatmn,rlatmx,rlonmn,rlonmx,dlat,dlon
c
      character *20 fnohc
c
c     Specify default value for ohc
      rmiss = 9999.
c
c     Specify unit number for input files
      luohc  = 61
c
      fnohc  = 'OHCYYYYMMDD.DAT'
c
c     Get the OHC array from the data file
      if (ifirst .eq. 1) then
c        Open the OHC file
         open(unit=luohc,file=fnohc,status='old',form='formatted',
     +        err=920)
c
         call getohc(luohc,iyear,imon,iday,
     +               rlatmn,rlatmx,rlonmn,rlonmx,dlat,dlon,
     +               nlat,nlon,mxlat,mxlon,ohcll,ierrd)
         ifirst=0
         close(luohc)
         if (ierrd .ne. 0) then
c           Error reading the data file 
            ohc = rmiss
            nday=9999
            ierr=4
            return
         endif
      endif
c
c
      slon =  rlonmn
      slat =  rlatmn
c
c     Set time to 00 UTC since OHC files are multi-day composites
      itime = 0
c
c     Interpolate to find the ohc
      izp  = 0
      if (rlon .lt. 0) then
         rloneast = rlon + 360.0
      else
         rloneast = rlon
      endif
c
      call llintsp(ohcll,slon,slat,dlon,dlat,mxlon,mxlat,nlon,nlat,
     +             ohc,rloneast,rlat,izp,ierrll)
      if (ohc .lt. 0.0) ohc=rmiss
c
      if (ierrll .ne. 0) then
c        Error reading the data file 
         ohc = rmiss
         nday=9999
         ierr=4
         close(luohc)
         return
      endif
c
c     Normal return
      ierr=0
      close(luohc)
      return
c
  920 continue
c     Error opening the data file or slon,slat point outside of OHC domain
      ohc = rmiss
      nday=9999
      ierr=3
      close(luohc)
      return
c
      end 
      subroutine getohc(luohc,iyear,imon,iday,
     +           rlatmn,rlatmx,rlonmn,rlonmx,dlat,dlon,
     +           nlat,nlon,mxlat,mxlon,ohcll,ierrd)
c
c     This routine reads the file containing the OHC data
c     on an evenly spaced lat/lon grid.
c
      dimension ohcll(mxlon,mxlat)
c
c     Initialize the ohc array to the default
      do j=1,mxlat
         do i=1,mxlon
	    ohcll(i,j) = 9999.
         enddo
      enddo
c
c     Read the header line
      read(luohc,100,err=900,end=900) iyear,imon,iday,itime,
     +                                rlonmn,rlonmx,dlon,
     +                                rlatmn,rlatmx,dlat
  100 format(i4,1x,i2,i2,1x,i2,
     +       f9.3,f9.3,f9.3,
     +       f9.3,f9.3,f9.3)
c
c     Calculate the number of lat,lon points
      if (rlatmx .le. rlatmn .or. dlat .le. 0.0) then
         ierrd=1
         return
      endif
      if (rlonmx .le. rlonmn .or. dlon .le. 0.0) then
         ierrd=1
         return
      endif
c
      eps=0.0001
      nlat = 1 + nint( eps + (rlatmx-rlatmn)/dlat )
      nlon = 1 + nint( eps + (rlonmx-rlonmn)/dlon )
c
c     Read the OHC data
      do j=1,nlat
         do i=1,nlon
c           read(luohc,*,err=900,end=900) 
c    +           dum1,dum2,dum3,dum4,dum5,ohcll(i,j),dum7
c
ccc  110       format(50x,f5.1)
            read(luohc,110,err=900,end=900) ohcll(i,j)
  110       format(168x,f10.3)
         enddo
      enddo
c
      ierrd=0
      return
c
  900 continue
      ierrd=1
      return
c
      end