subroutine rsstclim(ityr,itmon,itday,tlat,tlon,rmiss,
     +                 crsst,ierr)
c
c     Version 2.1.0
c     Last updated: 13 May 2020
c
c     This routine calculates the climatological values of  
c     sea surface temperature.
c
c
c     Input:  ityr  - 4-digit year
c             itmon - month (1-12)
c             itday - day (1-31)
c             tlat  - TC latitude (deg N)
c             tlon  - TC longitude (0-360degE)
c             rmiss - error value (real)
c
c     Output: 
c             crsst  - Climatological SST value
c             ierr - error flag =0 for normal completion
c                               =2 for illegal itmon
c                               =3 for illegal itday
c                               =4 for problem with reading input file
c                               =5 clim sst value missing (likely over land)
c
c     11 Oct 2011 - Modified, increase mx to 360 and my to 180 to
c                   accomodate global 1deg/1deg climo array
c                                
c     15 Aug 2014 - ioper=1 option added to add path to data input file for WCOSS
c
c     10 Mar 2020 - removed ioper flag
c     
c     13 May 2020 - updated rsstclim to 1982-2018 version
c
c-----------------------------------------------------------------------
c
c     Array for day of month represented by each monthly mean 
c     (used for time interpolation 0=dec,1=jan,...,12=dec,13=jan)
      dimension ndmo(0:13),jdmo(0:13),kdmo(12)
c
c     Inputs: Climatological SST
      parameter (mx=2000,my=2000)      
      real rlonl,rlonr,dlon
      real rlatt,rlatb,dlat
      dimension climsst(mx,my,12)
      character *3 label

      data luinp /77/
      character*256 coef_location, fninp
c
c     Month information
      dimension ngm(12)
c
      data kdmo /   31,28,31,30,31,30,31,31,30,31,30,31   /
      data ndmo /16,16,14,16,15,16,15,16,16,15,16,15,16,16/
c
c     Flag for file reading
      data irfread /1/
c
      save irfread,ndmo,jdmo,kdmo,ngm
      save rlonl,rlonr,dlon,nlon,
     +     rlatt,rlatb,dlat,nlat 
      save climsst
c
      ierr = 0
      crsst = rmiss
c
c     Calculate Julian Days of center of each month
c     for later use by time interpolation routines
      ityrt=1999
      do k=1,12
         call jday(k,ndmo(k),ityrt,jdmo(k))
      enddo
c
      jdmo(13) = 365 + ndmo(13)
      jdmo (0) =     - ndmo( 0)
c
      if (irfread .eq. 1) then
         irfread = 0
c        Open and read climatological SST file
c
         call getenv( "SHIPS_COEF", coef_location )
         fninp = trim( coef_location )//'/'//'rsstclim_1982_2018.dat'
c
         open(file=fninp,
     +        unit=luinp,form='formatted',status='old',err=900)
c
c        Read the header line
         read(luinp,505,err=900) rlonl,rlonr,dlon,nlon,
     +                           rlatb,rlatt,dlat,nlat,
     +                           rmissc
c
 505     format(8x,2(3(f9.2),i6),f9.1)	 
c
c        Check that the size of the climo grid in the .dat file
c        doesnt exceed the size of the working grid, climohc
         if (nlon .gt. mx .or. nlat .gt. my) then
	    print*, 'nlon/nlat exceeds mx/my: ',nlon,nlat,mx,my
	    stop
	 endif
c	 
c        Read the clim RSST data for each month
         do m=1,12
            read(luinp,100,err=900) label
  100       format(a3)
c
            do j=1,nlat
               read(luinp,110,err=900) (climsst(i,j,m),i=1,nlon)
  110          format(500(f6.1))
            enddo
         enddo
c
         close(luinp)
      endif
c
c     Determine the index of the climatological SST value
c        closest to the TC lat/lon:
      itc = nint((tlon-rlonl)/dlon)+1
      jtc = nint((tlat-rlatb)/dlat)+1 
c
c     Check that lat/lon are within the climo grid domain
      if ((itc .le. 0) .or. (itc .gt. nlon)) then
         ierr = 1
	 return
      endif
      
      if ((jtc .le. 0) .or. (jtc .gt. nlat)) then
         ierr = 1
	 return
      endif       
c
c     Check input dates
      if (itmon .gt. 12 .or. itmon .lt. 1) then
         ierr=2
         return
      endif
c
      ndmax = 2*ndmo(itmon)
      if (ndmax .gt. 31) ndmax=31
      if (itmon .eq. 2 .and. mod(ityr,4) .eq. 0) ndmax=29
c
      if (itday .gt. ndmax .or. itday .lt. 1) then
         ierr=3
         return
      endif
c
c     Calculate Julian day of requested input date
c     but under the assumption that it is not a leap year
      ityrt=1999
      itdaytm=itday
      if (itmon .eq. 2 .and. itday .eq. 29) itdaytm=28
      call jday(itmon,itdaytm,ityrt,ijday)
c
c     Calculate the weights for time interpolation
      if (itday .le. ndmo(itmon)) then
         im = itmon-1
         ip = itmon
      else
         im = itmon
         ip = itmon+1
      endif
      jm = jdmo(im)
      jp = jdmo(ip)
c
      dday = float(jp-jm)
      wtm = float(jp-ijday)/dday
      wtp = float(ijday-jm)/dday
c
c     Adjust indices for Jan and Dec
      if (im .eq.  0) im = 12
      if (ip .eq. 13) ip= 1
c
c     Perform time interpolation
c
      if (climsst(itc,jtc,im) .ne. rmissc .and. 
     +    climsst(itc,jtc,ip) .ne. rmissc) then    
         crsst = wtm*climsst(itc,jtc,im) + wtp*climsst(itc,jtc,ip)      
      else
         ierr = 5
      endif
c
      return
c
  900 ierr=4
      return
c
      end