    SUBROUTINE SSTHIRES(SST, SM, GLAT1, GLON1, IDAT, IM, JM, NM, ND)
!--------------------------------------------------------------------------------------------------
    IMPLICIT NONE   
!
    INTEGER(KIND=4)                                                                             ::&
    & IM      , JM      , NM      , ND      , T       , INSST   , INDXST
!
    REAL   (KIND=4)                                                                             ::&
    & ILAT1   , ILON1   ,                                                                         &
    & ILAT2   , ILON2
!
    REAL   (KIND=4)     , PARAMETER :: INCR_LAT = 0.0833333
    REAL   (KIND=4)     , PARAMETER :: INCR_LON = 0.0833333
!
    REAL   (KIND=4)     , PARAMETER :: H90  =  90.0
    REAL   (KIND=4)     , PARAMETER :: H360 = 360.0
    REAL   (KIND=4)     , PARAMETER :: D5   =   5.E-1
    REAL   (KIND=4)     , PARAMETER :: D00  =   0.0
    REAL   (KIND=4)     , PARAMETER :: H1   =   1.0
!
    REAL   (KIND=4)     , PARAMETER :: DTR  =   3.1415926 / 180.
!
    INTEGER(KIND=4)     , DIMENSION(4)                                                          ::&
    & IDATE
!
    INTEGER(KIND=4)     , DIMENSION(3)                                                          ::&
    & IDAT
!
    INTEGER(KIND=4)     , DIMENSION(12)                                                         ::&
    & MONTH
!
    REAL   (KIND=4)     , DIMENSION(4321,2160)                                                  ::&
    & SSTLL
!
    REAL   (KIND=4)     , DIMENSION(12)                                                         ::&
    & SALTLK
!------------------------------------------
! DRAGAN, ATUALIZADO POR: GSM EM 2017-02-10
!------------------------------------------
    REAL   (KIND=4)     , DIMENSION(2)                                                          ::&
    & SALTLA  , SALTLO  
!
    REAL   (KIND=4)     , DIMENSION(0:IM+1, 0:JM+1, NM)                                         ::&
    & SST     , SM
!
    REAL   (KIND=4)     , DIMENSION(0:IM+1, 0:JM+1, NM)                   , INTENT(IN)          ::&
    & GLAT1   , GLON1
!
    REAL   (KIND=4)     , DIMENSION(0:IM+1, 0:JM+1, NM)                                         ::&
    & GLAT    , GLON
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
    INTEGER(KIND=4)                                                                             ::&
    & IERR    , N       , I       , J       , ID1     , ID2     , MARG0   , MNTH0   , MNTH1   ,   &
    & NUMER   , DENOM   , IARG1   , IARG2   , LON1INDX, LON2INDX, LAT1INDX, LAT2INDX
!
    REAL   (KIND=4)                                                                             ::&
    & ELAT    , ELON    , DIF     , W1      , W2      , AR1     , AR2     , AR3     , AR4     ,   &
    & FRAC
!
    DATA   INSST  /39/
    DATA   INDXST / 0/
!
    DATA MONTH    /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
!
    DATA SALTLK   /273.38, 274.27, 278.50, 283.01, 287.33, 293.41,                                &
    &              297.13, 297.73, 294.97, 289.58, 282.31, 275.67/
!---------------------------------------------------------
! CORNERS OF SALT LAKE LATITUDE / LONGITUDE BOX IN DEGREES 
!                     40.0      42.0
!---------------------------------------------------------
    DATA SALTLA   /0.698132, 0.733038/
!---------------------------------------------------------
!                     111.0    114.0
!---------------------------------------------------------
    DATA SALTLO   /1.937315, 1.989675/
!
    GLON = -GLON1
    GLAT =  GLAT1
!
    WRITE(0,*) "CALL GRIDST FOR READING SST"
!
    CALL GRIDST(INSST,SSTLL)
!						
      ierr = 0
      write(0,*)'IERR1=',ierr

!cjfp ****** final desta modificacao                         !dragan

      write(0,*) sstll(1,1)  
      write(0,*) sstll(100,100)                           !dragan
!GSM      write(0,*) 'msst values' 
      
!!!!!!!! print*,"sstll",sstll
    
							
      if (ierr.ne.0) goto 4500
!
!----  interpolate 1-deg global sst to eta grid  -------
!
!-cp note:  this subroutine and interpolation algorithm assume
!-cp a 1-deg global sst field in the following format:  
!-cp
!-cp  i=1 at 0.5 e,  i=2 at 1.5 e, ... , i=360 at 0.5w
!-cp  j=1 at 89.5s, j=2 at 88.5 s, ..., j=180 at 89.5n
!
!
!        new 0.5 degree data
!
!        i=1 at 0.25 e, i=720 at 0.25 w, i=721 at 0.25 e
!        j=1 at 89.75 s, j=360 at 39.75 n 
!
!
!-cp  
!-cp in the interpolation algorithm below, glon is positive westward,
!-cp from 0 to 360, with 0 at the greenwich meridian.  elon is positive 
!-cp eastward, thus the need to subtract glon from 360 to get the index
!-cp of the correct oisst point.  if your input 1 deg sst field is in
!-cp a different indexing scheme, you will need to change the algorithm
!-cp below - see "grdeta.oldoi"
!-cp

!*** write the sst monthly data n a new file                !dragan, august 2011

     
      do n=1,nm
      do j=0,jm+1
      do i=0,im+1
    
        elat=h90+glat(i,j,n)/dtr
        elon=h360-glon(i,j,n)/dtr
        if(elon.gt.h360)elon=elon-h360

        dif=elon-int(elon)

        if (dif .lt. incr_lon) ilon1=int(elon)+(0.5*incr_lon)
	
        if (dif .ge. incr_lon .and. dif .lt. (2*incr_lon)) then 
            ilon1=int(elon)+(1.5*incr_lon)
        endif 
	
        if (dif .ge. (2*incr_lon) .and. dif .lt. (3*incr_lon)) then 
            ilon1=int(elon)+(2.5*incr_lon)
        endif
								
        if (dif .ge. (3*incr_lon) .and. dif .lt. (4*incr_lon)) then 
            ilon1=int(elon)+(3.5*incr_lon)
        endif 
	
        if (dif .ge. (4*incr_lon) .and. dif .lt. (5*incr_lon)) then 
            ilon1=int(elon)+(4.5*incr_lon)
        endif
	
        if (dif .ge. (5*incr_lon) .and. dif .lt. (6*incr_lon)) then 
            ilon1=int(elon)+(5.5*incr_lon)
        endif

        if (dif .ge. (6*incr_lon) .and. dif .lt. (7*incr_lon)) then 
            ilon1=int(elon)+(6.5*incr_lon)
        endif

        if (dif .ge. (7*incr_lon) .and. dif .lt. (8*incr_lon)) then 
            ilon1=int(elon)+(7.5*incr_lon)
        endif
	
        if (dif .ge. (8*incr_lon) .and. dif .lt. (9*incr_lon)) then 
            ilon1=int(elon)+(8.5*incr_lon)
        endif
	
        if (dif .ge. (9*incr_lon) .and. dif .lt. (10*incr_lon)) then 
            ilon1=int(elon)+(9.5*incr_lon)
        endif
	
        if (dif .ge. (10*incr_lon) .and. dif .lt. (11*incr_lon)) then 
            ilon1=int(elon)+(10.5*incr_lon)
        endif
	
        if (dif .ge. (11*incr_lon)) then 
            ilon1=int(elon)+(11.5*incr_lon)
        endif


        if(ilon1.le.d00)ilon1=360.
        ilon2=ilon1+incr_lon

!
!-mp        new approach sets ilat1, ilon1 to point on sst grid that is
!-mp        sw of the eta grid point.
!
!
        dif=elat-int(elat)
!GSM        write(6,*) 'dif_lat: ',dif
!GSM        if (dif .ge. 0.75) ilat1=int(elat)+0.75
!GSM        if (dif .lt. 0.25) ilat1=int(elat)-0.25
!GSM        if (dif .ge. 0.25 .and. dif .lt. 0.75) ilat1=int(elat)+0.25


	      IF (dif .lt. incr_lat) ilat1=int(elat)+(0.5*incr_lat)

	      IF (dif .ge. incr_lat .and. dif .lt. (2*incr_lat)) then 
	         ilat1=int(elat)+(1.5*incr_lat)
	      endif 

	      IF (dif .ge. (2*incr_lat) .and. dif .lt. (3*incr_lat)) then 
	         ilat1=int(elat)+(2.5*incr_lat)
	      endif

	      IF (dif .ge. (3*incr_lat) .and. dif .lt. (4*incr_lat)) then 
	         ilat1=int(elat)+(3.5*incr_lat)
	      endif 

	      IF (dif .ge. (4*incr_lat) .and. dif .lt. (5*incr_lat)) then 
	         ilat1=int(elat)+(4.5*incr_lat)
	      endif

	      IF (dif .ge. (5*incr_lat) .and. dif .lt. (6*incr_lat)) then 
	         ilat1=int(elat)+(5.5*incr_lat)
	      endif

	      IF (dif .ge. (6*incr_lat) .and. dif .lt. (7*incr_lat)) then 
	         ilat1=int(elat)+(6.5*incr_lat)
      	      endif

	      IF (dif .ge. (7*incr_lat) .and. dif .lt. (8*incr_lat)) then 
	         ilat1=int(elat)+(7.5*incr_lat)
	      endif

	      IF (dif .ge. (8*incr_lat) .and. dif .lt. (9*incr_lat)) then 
	         ilat1=int(elat)+(8.5*incr_lat)
	      endif

	      IF (dif .ge. (9*incr_lat) .and. dif .lt. (10*incr_lat)) then 
	         ilat1=int(elat)+(9.5*incr_lat)
   	      endif

	      IF (dif .ge. (10*incr_lat) .and. dif .lt. (11*incr_lat)) then 
	         ilat1=int(elat)+(10.5*incr_lat)
	      endif

	      IF (dif .ge. (11*incr_lat)) then 
	         ilat1=int(elat)+(11.5*incr_lat)
	      endif

!      dif=elat-ilat1
!hires      if(dif.gt.d5)ilat1=min(ilat1+1,179)
!tst      if(dif.gt.incr/2.)ilat1=amin1((ilat1+incr),179.5)


        ilat2=ilat1+incr_lat

        w1=elon-ilon1+incr_lon/2.
        if(w1.lt.d00)w1=w1+h360
        w2=elat-ilat1+incr_lat/2.
        ar1=w1*w2
        ar2=w1*(h1-w2)
        ar3=(h1-w1)*(h1-w2)
        ar4=(h1-w1)*w2
!        lon1indx=2*ilon1+1
!        lon2indx=2*ilon2+1
!        lat1indx=2*ilat1+1
!        lat2indx=2*ilat2+1
        lon1indx=12*(ilon1)
        lon2indx=12*(ilon2)
        lat1indx=12*(ilat1)
        lat2indx=12*(ilat2)
!      sst(i,j) = ar1*sstll(ilon2,ilat2)+ar2*sstll(ilon2,ilat1)+
!     1            ar3*sstll(ilon1,ilat1)+ar4*sstll(ilon1,ilat2)
        
        if (lon1indx .lt. 1 .or. lon1indx .gt. 4321) then
        write(0,*) 'out of bounds on index!!', lon1indx
        endif

!***********************************************************************************
!if somebody wants to use some other way to get sst, from some binary or grib file,          !dragan, april, 2012
!then, uncomment this, comment call to sststp, and modify gribst.f90
!***********************************************************************************

!!!!!!        sst(i,j,n)=ar1*sstll(lon2indx,lat2indx)+  &
!!!!!!                        ar2*sstll(lon2indx,lat1indx)+  &
!!!!!!                        ar3*sstll(lon1indx,lat1indx)+  &
!!!!!!                        ar4*sstll(lon1indx,lat2indx)
!!!!!!        allocate(msstt(0:im+1,0:jm+1,nm,t))         
!***********************************************************************************

        
               sst2(i,j,n)=ar1*sstll(lon2indx,lat2indx)+ar2*sstll(lon2indx,lat1indx)+ar3*sstll(lon1indx,lat1indx)+ar4*sstll(lon1indx,lat2indx)
	       sst(i,j,n)=sstll(lon2indx,lat2indx)
               write(1004,*) sst(i,j,n)     
      enddo         
      enddo
      enddo
        write(1003,*) sst        

!***
!***  insert temperatures for the great salt lake
!***
      id1=idat(1)
      id2=idat(2)+nd
      marg0=id1-1
      if(marg0.lt.1)marg0=12
      mnth0=month(marg0)
      mnth1=month(id1)
      if(id2.lt.15)then
        numer=id2+mnth0-15
        denom=mnth0
        iarg1=marg0
        iarg2=id1
      else
        numer=id2-15
        denom=mnth1
        iarg1=id1
        iarg2=id1+1
        if(iarg2.gt.12)iarg2=1
      endif
      frac=numer/denom
      do n=1,nm
      do j=1,jm
      do i=1,im
        if(glat(i,j,n).gt.saltla(1).and.glat(i,j,n).lt.saltla(2))then
          if(glon(i,j,n).gt.saltlo(1).and.glon(i,j,n).lt.saltlo(2))then
            if(sm(i,j,n).gt.0.5)  &
              sst(i,j,n)=saltlk(iarg1)+  &
                      (saltlk(iarg2)-saltlk(iarg1))*frac
          endif
        endif
      enddo
      enddo
      enddo
!
      return
!
 4500 continue
!              error occurred when inputing sst from grib.
      write (0, 4550) insst
 4550 format ('0', 'error occurred when reading in sst        ',  &
                   'on unit', i3, ' grib ' /  &
              ' ', 'execution terminating.')
!
      stop 222
!
      end subroutine ssthires
