C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
      SUBROUTINE SST_ECMWF (SST,SM,GLAT,GLON,IDAT,LIST,DTR)
C
      IMPLICIT REAL (A-H, O-Z)
C
      INCLUDE "parmeta"
C

Chires
	REAL INCR,ILAT1,ILON1,ILAT2,ILON2
        PARAMETER(INCR=0.5)
Chires

      PARAMETER  (H90=90.0,H360=360.0,D5=5.E-1,D00=0.0,H1=1.0)
       PARAMETER (NX=201,NY=101,DX=0.5,DY=0.5)        !cspir
       PARAMETER (west=-40,east=60,north=60,south=10) !cspir
C
      INTEGER IDATE(4),IDAT(3),MONTH(12)
      DIMENSION SSTLL(721,361),SALTLK(12),SALTLA(2),SALTLO(2)
      DIMENSION SMLL(721,361)! ioannis
!      integer nx,ny ! ioannis
!      real DX,DY,west,east,south,north ! ioannis
c      DIMENSION SSTTEMP(245,111),SMTEMP(245,111)  ! ioannis
ccc      DIMENSION TEMP(221,141)  ! ioannis
!      REAL,ALLOCATABLE:: SSTTEMP(:,:),SMTEMP(:,:) !ioannis
C
      DIMENSION SSTTEMP(NX,NY),SMTEMP(NX,NY) !cspir
      DIMENSION  SST(IM,JM), SM(IM,JM), GLAT(IM,JM), GLON(IM,JM)
C
      CHARACTER*140 INFN,INFNSM  ! ioannis
C
      DATA   INSST/39/
      DATA   INDXST/0/

      DATA MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
C
      DATA SALTLK/273.38,274.27,278.50,283.01,287.33,293.41
     1,           297.13,297.73,294.97,289.58,282.31,275.67/
C
C     CORNERS OF SALT LAKE LATITUDE/LONGITUDE BOX
C     in degrees---> 40.0     42.0            111.0    114.0
      DATA SALTLA/0.698132,0.733038/,SALTLO/1.937315,1.989675/
C
      IOUTUPRT = LIST
c--------------------------------------------------------------------
cioannis      CALL GRIBST(INSST,INDXST,SSTLL,IERR)
cioannis      IERR=0
c--------------------------------------------------------------------
c      IF (IERR.NE.0) GOTO 4500
c--------------------------------------------------------------------
!       ALLOCATE(SSTTEMP(nx,ny),SMTEMP(nx,ny))
       INFN='/vol1/meteo/ETA/worketa.v5.0/data/prep/dump'
         write(6,*)'ioannisSST',infn
       OPEN(UNIT=INSST,FILE=INFN(1:(INDEX(INFN,' '))-1),
     &STATUS='OLD',FORM='FORMATTED')
        read(INSST,*) mx, my
        write(6,*) mx, my
       DO J=NY,1,-1
        DO I=1,NX
            read(INSST,*) SSTTEMP(I,J)
        ENDDO
       ENDDO 
!      INFN='../../data/prep/06011006000.ETA_SST'
!      print*,'ioannisSST',infn
c       OPEN(UNIT=INSST,FILE=INFN(1:(INDEX(INFN,' '))-1),
c     &STATUS='OLD',FORM='FORMATTED')
c      read(INSST,*)NX,NY,DX,DY,west,east,south,north
c      print*,'SST - NX - NY=',nx,ny ! ioannis
c      print*,'SST - DX - DY=',dx,dy ! ioannis
c     IF (DX.NE.0.5.OR.DY.NE.0.5) THEN
c         STOP 997
C      ENDIF
!      ALLOCATE(SSTTEMP(nx,ny),SMTEMP(nx,ny))
c      read(INSST,*)           ! SICE
c      read(INSST,*) SSTTEMP   ! SST
c      read(INSST,*)           ! SD
c      CLOSE(INSST)

      DO J=1,NY
      DO I=1,NX
c        IF (SSTTEMP(I,J).LT.0.5) THEN
        IF (SSTTEMP(I,J).LT.250.OR.SSTTEMP(I,J).GT.350) THEN !cspir
            SMTEMP(i,j)=0.0   ! land
        ELSE
            SMTEMP(i,j)=1.0   ! sea
        ENDIF
      ENDDO
      ENDDO

      DO j=1,361
      DO i=1,721
       SSTLL(i,j)=-999.9
       SMLL(i,j)=-999.9
      ENDDO
      ENDDO

c --------------------------------------------------------------------
c SSTTEMP(1:nx,1:ny) temporary array with the SSTs in the ECMWF domain
c SSTLL(721,361) global array with resolution 0.5 deg. The array 
c                SSTTEMP will now be fitted inside SSTLL
c in SSTLL i=1 for 0E, i=2 for 0.5E, i=3 for 1E,...,i=720 for 359.5E
c in SSTLL j=1 for -90N, j=2 for -89.5N, ..., j=361 for 90N
c
c     jsouth = index of SSTLL at the southern limit of SSTTEMP
c     jnorth = index of SSTLL at the northern limit of SSTTEMP
c     iwest  = index of SSTLL at the western limit of SSTTEMP
c     ieast  = index of SSTLL at the eastern limit of SSTTEMP
c     that is, (jnorth-jsouth+1)=ny
c     while the total number of points between iwest and ieast is nx
c ---------------------------------------------------------------------

      jsouth=int((south+90.0)/0.5)+1 
      jnorth=int((north+90.0)/0.5)+1      
      if (west.lt.0.0) then
      iwest=int((360.0+west)/0.5)+1
      else
      iwest=int(west/0.5)+1
      endif

      if (east.lt.0.0) then
      ieast=int((360.0+east)/0.5)+1
      else
      ieast=int(east/0.5)+1
      endif

      print*,'SST coordinates'
      print*,'west,east,iwest,ieast=',west,east,iwest,ieast
      print*,'south,north,isouth,inorth=',south,north,jsouth,jnorth

      if ((west.lt.0.0).and.(east.lt.0.0)) then
      if (iwest.ge.ieast) STOP 992 
      print*,'SST domain case 1'      

      DO j=jsouth,jnorth
        do i=iwest,ieast
           SSTLL(i,j)=SSTTEMP(i-iwest+1,j-jsouth+1)
           SMLL(i,j)=SMTEMP(i-iwest+1,j-jsouth+1)
        enddo
      ENDDO

      elseif ((west.ge.0.0).and.(east.ge.0.0)) then
      if (iwest.ge.ieast) STOP 993 
      print*,'SST domain case 2'      

      DO j=jsouth,jnorth
        do i=iwest,ieast
           SSTLL(i,j)=SSTTEMP(i-iwest+1,j-jsouth+1)
           SMLL(i,j)=SMTEMP(i-iwest+1,j-jsouth+1)
        enddo
      ENDDO

      elseif ((west.lt.0.0).and.(east.ge.0.0)) then
      if (iwest.le.ieast) STOP 994 
      print*,'SST domain case 3'      

      DO j=jsouth,jnorth
        do i=iwest,720        ! west to -0.5E
           SSTLL(i,j)=SSTTEMP(i-iwest+1,j-jsouth+1)
           SMLL(i,j)=SMTEMP(i-iwest+1,j-jsouth+1)
        enddo
      ENDDO

      DO j=jsouth,jnorth
        do i=1,ieast          ! 0E to east
           SSTLL(i,j)=SSTTEMP(i+720-iwest+1,j-jsouth+1)
           SMLL(i,j)=SMTEMP(i+720-iwest+1,j-jsouth+1)
        enddo
      ENDDO

      else
      print*,'SST domain case 4'      
      STOP 996
      endif

c      DO j=jnorth,jsouth,-1   ! 321,211,-1 --- 70N to 15N
c        do i=1,141 ! 0E to 70E
c           SSTLL(i,j)=SSTTEMP(i+104,321-j+1)
c           SMLL(i,j)=SMTEMP(i+104,321-j+1)
c        enddo
c      ENDDO

c      DO j=jnorth,jsouth,-1   ! 321,211,-1  --- 70N to 15N
c        do i=617,720 ! -52E to -0.5E
c           SSTLL(i,j)=SSTTEMP(i-616,321-j+1)
c           SMLL(i,j)=SMTEMP(i-616,321-j+1)
c        enddo
c      ENDDO

      do j=1,361  ! Fix greenwich meridian
        SSTLL(721,j)=SSTLL(1,j)
        SMLL(721,j)=SMLL(1,j)
      enddo

!      DEALLOCATE(SSTTEMP) ! ioannis
!      DEALLOCATE(SMTEMP)  ! ioannis

C-CP  
C-CP In the interpolation algorithm below, glon is positive westward,
C-CP from 0 to 360, with 0 at the greenwich meridian.  Elon is positive 
C-CP eastward, thus the need to subtract glon from 360 to get the index
C-CP of the correct point.  
C-CP
      print*,'ioannis SST 0.5 degrees'
      DO J=1,JM
      DO I=1,IM
      ELAT=H90+GLAT(I,J)/DTR
      ELON=H360-GLON(I,J)/DTR
      IF(ELON.GT.H360)ELON=ELON-H360
      print*,"IOANNIS ELON, ELAT=",i,j,ELON,ELAT

       DIF=ELON-INT(ELON)
       IF (DIF .ge. 0.5) ILON1=INT(ELON)+0.5 ! ioannis
       IF (DIF .lt. 0.5) ILON1=INT(ELON)     ! ioannis
cioannis      IF(ILON1.LE.D00) ILON1=360.
      IF(ILON1.LT.D00) THEN  ! ioannis
        print*,'IOANNIS ERROR1'
        STOP 777
      ENDIF
      ILON2=ILON1+INCR
c      print*,"IOANNIS DIF,ILON1,ILON2=",DIF,ILON1,ILON2

C
C-MP	New approach sets ILAT1, ILON1 to point on SST grid that is
C-MP	SW of the Eta Grid point.
C
C
        DIF=ELAT-INT(ELAT)
        IF (DIF .ge. 0.5) ILAT1=INT(ELAT)+0.5 ! ioannis
        IF (DIF .lt. 0.5) ILAT1=INT(ELAT)     ! ioannis
       ILAT2=ILAT1+INCR
c      print*,"IOANNIS DIF,ILAT1,ILAT2=",DIF,ILAT1,ILAT2

c      W1=ELON-ILON1+INCR/2.
c      IF(W1.LT.D00)W1=W1+H360
c      W2=ELAT-ILAT1+INCR/2.
c      AR1=W1*W2
c      AR2=W1*(H1-W2)
c      AR3=(H1-W1)*(H1-W2)
c      AR4=(H1-W1)*W2
c
      W1=(ELON-ILON1)/INCR
      W2=(ELAT-ILAT1)/INCR
      AR1=W1*W2
      AR2=W1*(H1-W2)
      AR3=(H1-W1)*(H1-W2)
      AR4=(H1-W1)*W2
c
      LON1INDX=2*ILON1+1
      LON2INDX=2*ILON2+1
      LAT1INDX=2*ILAT1+1
      LAT2INDX=2*ILAT2+1
      print*,"IOANNIS indices=",LON1INDX,LON2INDX,LAT1INDX,LAT2INDX

      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,*) '------------------------------------------'
C	write(6,*) 'corresponding indices: ', LAT1INDX,LON1INDX,
C     +					      LAT2INDX,LON2INDX
      endif

       if (LON1INDX .lt. 1 .or. LON1INDX .gt. 721) then
        write(6,*) 'out of bounds on index!!', LON1INDX
        STOP 778
       endif
       if (LON2INDX .lt. 1 .or. LON2INDX .gt. 721) then
        write(6,*) 'out of bounds on index!!', LON2INDX
        STOP 779
       endif

       SMCHECK=SMLL(LON1INDX,LAT1INDX)+SMLL(LON2INDX,LAT1INDX)+
     $         SMLL(LON1INDX,LAT2INDX)+SMLL(LON2INDX,LAT2INDX)

       IF (SMCHECK.EQ.4.0) THEN ! all 4 nearest points are sea-points
       SST(I,J)=AR1*SSTLL(LON2INDX,LAT2INDX)+
     +   AR2*SSTLL(LON2INDX,LAT1INDX)+
     +   AR3*SSTLL(LON1INDX,LAT1INDX)+
     +   AR4*SSTLL(LON1INDX,LAT2INDX)
       ELSEIF ((SMCHECK.LT.4.0).AND.(SMCHECK.GT.0.0)) THEN ! at least 1 points is sea 
       SST(I,J)=(SMLL(LON2INDX,LAT2INDX)*SSTLL(LON2INDX,LAT2INDX)+
     +   SMLL(LON2INDX,LAT1INDX)*SSTLL(LON2INDX,LAT1INDX)+
     +   SMLL(LON1INDX,LAT1INDX)*SSTLL(LON1INDX,LAT1INDX)+
     +   SMLL(LON1INDX,LAT2INDX)*SSTLL(LON1INDX,LAT2INDX))/SMCHECK
       ELSEIF (SMCHECK.EQ.0.0) THEN ! all 4 nearest points are land
       SST(I,J)=(COS(glat(i,j))**2.)*30.+269.
c       print*,'!!!!=',SMCHECK,SM(i,j),sst(i,j),elat,glat(i,j)
c     $,cos(glat(i,j))
       ENDIF

      ENDDO
      ENDDO

c--------------------------------------------------------------------
c--------------------------------------------------------------------
C***
C***  INSERT TEMPERATURES FOR THE GREAT SALT LAKE
C***
c      ID1=IDAT(1)
c      ID2=IDAT(2)
c      MARG0=ID1-1
c      IF(MARG0.LT.1)MARG0=12
c      MNTH0=MONTH(MARG0)
c      MNTH1=MONTH(ID1)
c      IF(ID2.LT.15)THEN
c        NUMER=ID2+MNTH0-15
c        DENOM=MNTH0
c        IARG1=MARG0
c        IARG2=ID1
c      ELSE
c        NUMER=ID2-15
c        DENOM=MNTH1
c        IARG1=ID1
c        IARG2=ID1+1
c        IF(IARG2.GT.12)IARG2=1
c      ENDIF
c      FRAC=NUMER/DENOM
c      DO J=1,JM
c      DO I=1,IM
c        IF(GLAT(I,J).GT.SALTLA(1).AND.GLAT(I,J).LT.SALTLA(2))THEN
c          IF(GLON(I,J).GT.SALTLO(1).AND.GLON(I,J).LT.SALTLO(2))THEN
c            IF(SM(I,J).GT.0.5)
c     1        SST(I,J)=SALTLK(IARG1)+
c     2                (SALTLK(IARG2)-SALTLK(IARG1))*FRAC
c          ENDIF
c        ENDIF
c      ENDDO
c      ENDDO
C
      RETURN
C
 4500 CONTINUE
C              ERROR OCCURRED WHEN INPUTING SST FROM GRIB.
      WRITE (IOUTUPRT, 4550) INSST
 4550 FORMAT ('0', 'ERROR OCCURRED WHEN READING IN SST        ',
     1             'ON UNIT', I3, ' GRIB ' /
     2        ' ', 'EXECUTION TERMINATING.')
C
      STOP 222
C
      END
