!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine ssthires (sst,sm,glat1,glon1,idat,im,jm,nm,nd)

      implicit real (a-h, o-z)

      integer::im,jm,nm,nd, t
      real incr,ilat1,ilon1,ilat2,ilon2
      parameter(incr=0.5)

      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(721,360),saltlk(12),saltla(2),saltlo(2)   !dragan
    
!!!!!!!!!!!     real,dimension(:,:,:,:),allocatable:: msstt

      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)


      real:: msst(722,360,356), msstt(0:im+1,0:jm+1,nm)     !dragan

      character(len=256)::monthlysst
      character(len=4)::cn

      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/
!
      glon=-glon1
      glat=glat1
      
!cjfp ****** incluido uma nova matriz na chamada da funcao (msst)  !dragan
      
!!!!      call gribst(insst,indxst,sstll,ierr,nd,msst)                !dragan-added msst
    
    
      print*,'IERR=',ierr
      write(0,*)' before gribst in ssthires'
      call gribst(ierr,msst)        ! dragan-removed sstll, sst is got from monthly sst archive 1981(dec)-2011(july), 
                                                  ! calling subroutine sstch.f90

      write(0,*)'IERR=',ierr
!   call gribst(insst,indxst,sstll,ierr,nd)

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

      write(0,*) msst(1,1,1)  
      write(0,*) msst(100,100,100)                           !dragan
      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 t=1,356          !dragan, august 2011

        write(cn,'(i4.4)') t
        monthlysst='/home/lucci/GEF/GEF-reta/GEF_tempo/scratchout/initdata/monthly_sst'//cn//'.bin'  !dragan, august 2011
        open(33,file=monthlysst,form='unformatted',status='unknown',access='sequential')  
     
      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.75) ilon1=int(elon)+0.75
        if (dif .lt. 0.25) ilon1=int(elon)-0.25
        if (dif .ge. 0.25 .and. dif .lt. 0.75) ilon1=int(elon)+0.25

        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.75) ilat1=int(elat)+0.75
        if (dif .lt. 0.25) ilat1=int(elat)-0.25
        if (dif .ge. 0.25 .and. dif .lt. 0.75) ilat1=int(elat)+0.25


!      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

        w1=elon-ilon1+incr/2.
        if(w1.lt.d00)w1=w1+h360
        w2=elat-ilat1+incr/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=2*(ilon1+incr/2.)
        lon2indx=2*(ilon2+incr/2.)
        lat1indx=2*(ilat1+incr/2.)
        lat2indx=2*(ilat2+incr/2.)
!      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. 721) then
        write(0,*) 'out of bounds on index!!', lon1indx
        endif
        if (lat1indx .lt. 1 .or. lat1indx .gt. 360) then
!d        write(0,*)' lat1indx = ', lat1indx, ilat1, incr
        lat1indx=min(360,max(1,lat1indx))
        end if
        if (lat2indx .lt. 1 .or. lat2indx .gt. 360) then
!d        write(0,*)' lat2indx = ', lat2indx, ilat2, incr
        lat2indx=min(360,max(1,lat2indx))
        end if

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

        
        msstt(i,j,n)=ar1*msst(lon2indx,lat2indx,t)+ar2*msst(lon2indx,lat1indx,t)+  &
                 ar3*msst(lon1indx,lat1indx,t)+ar4*msst(lon1indx,lat2indx,t)                           !dragan, august 2011
        
	if(msstt(i,j,n).gt.350.or.msstt(i,j,n).lt.150)then
	print*,'msstt>350=',msstt(i,j,n), i,j,n
        endif
     
     
      enddo         
      enddo
      enddo
        
        write(33) msstt               !dragan, august 2011           
        close(33)            
 !!!!     deallocate(msstt)
      enddo


!***********************************      
       call sststp(sst)                       !dragan-getting sst from monthly sst archive 1981(dec)-2011(july)  
!***********************************      

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