!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
      SUBROUTINE SSTHIRES2 (SST,SM,GLAT1,GLON1,IDAT,im,jm,nm,nd)
!
      IMPLICIT REAL (A-H, O-Z)
!
!

!hires
      integer::im,jm,nm,nd
	REAL INCR,ILAT1,ILON1,ILAT2,ILON2
	PARAMETER(INCR=0.25)
!hires

      PARAMETER  (H90=90.0,H360=360.0,D5=5.E-1,D00=0.0,H1=1.0)
      real,parameter::dtr=3.1415926/180.
!
      INTEGER IDATE(4),IDAT(3),MONTH(12)
      DIMENSION SSTLL(1441,720),SALTLK(12),SALTLA(2),SALTLO(2)
!
      DIMENSION  SST(0:IM+1,0:JM+1,nm), SM(0:IM+1,0:JM+1,nm)  
      real,intent(in)::GLAT1(0:IM+1,0:JM+1,nm), GLON1(0:IM+1,0:JM+1,nm)
      real::GLAT(0:IM+1,0:JM+1,nm), GLON(0:IM+1,0:JM+1,nm)
!
      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            111.0    114.0
      DATA SALTLA/0.698132,0.733038/,SALTLO/1.937315,1.989675/
!
      character(len=150)::filename
      glon=-glon1
      glat=glat1
      filename='/scratchin/grupos/grpeta/projetos/tempo/oper/gef_v1.0.0/gef_trunk/PRP/data_in/init/tmisst'
      CALL read_rss_sst(filename,sstll(1:1440,:))
      sstll(1441,:)=sstll(1,:)
      do i=1,1441
      do j=1,720
         if(sstll(i,j).ne.0)then
!	   print *,i,j,sstll(i,j)
         endif
      enddo
      enddo
      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.25 degree data
!
!	I=1 at 0.125 E, I=1440 at 0.125 W, I=1441 at 0.125 E
!	J=1 at 89.875 S, J=360 at 89.875 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
      do n=1,nm
      DO J=1,JM
      DO I=1,IM
      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 .ge. 0.875) ILON1=INT(ELON)+0.875
	IF (DIF .lt. 0.125) ILON1=INT(ELON)-0.125
	IF (DIF .ge. 0.125 .and. DIF .lt. 0.875) ILON1=INT(ELON)+0.125

!old      IF(ILON1.EQ.D00)ILON1=360.
      IF(ILON1.LE.D00)ILON1=360.
      ILON2=ILON1+INCR

!
!-MP	New approach sets ILAT1, ILON1 to point on SST grid that is
!-MP	SW of the Eta Grid point.
!
!
        DIF=ELAT-INT(ELAT)
        IF (DIF .ge. 0.875) ILAT1=INT(ELAT)+0.875
        IF (DIF .lt. 0.125) ILAT1=INT(ELAT)-0.125
        IF (DIF .ge. 0.125 .and. DIF .lt. 0.875) ILAT1=INT(ELAT)+0.125


!      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)


!     IF(ILAT1.EQ.180.OR.ILAT1.EQ.0)THEN
!       WRITE(6,6788)I,J,GLAT(I,J),GLON(I,J),ELAT,ELON
!6788   FORMAT(' I,J=',2I4,' GLAT=',E12.5,' GLON=',E12.5,
!    1   ' ELAT=',E12.5,' ELON=',E12.5)
!       STOP 333
!     ENDIF

      ILAT2=ILAT1+INCR

!hires,notsure      W1=ELON-ILON1+D5
      W1=ELON-ILON1+INCR/4.
      IF(W1.LT.D00)W1=W1+H360
!hires,notsure      W2=ELAT-ILAT1+D5
      W2=ELAT-ILAT1+INCR/4.
      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=4*(ILON1+INCR/4.)
	LON2INDX=4*(ILON2+INCR/4.)
	LAT1INDX=4*(ILAT1+INCR/4.)
	LAT2INDX=4*(ILAT2+INCR/4.)
	if (mod (I,20) .eq. 0 .and. mod(J,20) .eq. 0) then
!	write(6,*) 'weights: ',AR1,AR2,AR3,AR4
!	write(6,*) 'ILAT1,ILON1,ELAT,ELON: ', ILAT1,ILON1,ELAT,ELON
!	write(6,*) '------------------------------------------'
!	write(6,*) 'corresponding indices: ', LAT1INDX,LON1INDX,
!     +					      LAT2INDX,LON2INDX
	endif
!      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. 1441) then
!	write(6,*) 'out of bounds on index!!', LON1INDX,ilon1,incr
	endif

        if(sstll(LON2INDX,LAT2INDX).lt.500.and.   &
	   SSTLL(LON2INDX,LAT1INDX).lt.500.and.  &
           SSTLL(LON1INDX,LAT1INDX).lt.500.and.  &
           SSTLL(LON1INDX,LAT2INDX).lt.500)then
	SST(I,J,n)=AR1*SSTLL(LON2INDX,LAT2INDX)+  &
      	 	 AR2*SSTLL(LON2INDX,LAT1INDX)+  &
      	 	 AR3*SSTLL(LON1INDX,LAT1INDX)+  &
      	 	 AR4*SSTLL(LON1INDX,LAT2INDX)
        else
	  sst(i,j,n)=min(sstll(LON2INDX,LAT2INDX),  &
      	 	 SSTLL(LON2INDX,LAT1INDX),  &
      	 	 SSTLL(LON1INDX,LAT1INDX),  &
      	 	 SSTLL(LON1INDX,LAT2INDX))
        endif
      ENDDO
      ENDDO
      ENDDO
!***
!***  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

      END