!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      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_lat,inc_lon,ilat1,ilon1,ilat2,ilon2

      parameter(incr_lat=0.0833333)
      parameter(incr_lon=0.0833333)

      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(4321,2160),saltlk(12),saltla(2),saltlo(2)   !dragan,  Atualizado Por: GSM em 20170210
    
!!!!!!!!!!!     real,dimension(:,:,:,:),allocatable:: msstt

      dimension  sst(0:im+1,0:jm+1,nm),sm(0:im+1,0:jm+1,nm)  
      dimension  sst2(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)


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

!GSM      character(len=256)::monthlysst
!GSM      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
    
    
!GSM      print*,'IERR=',ierr
      write(0,*)' before gridst in ssthires'
!GSM      call gribst(ierr,msst)        ! dragan-removed sstll, sst is got from monthly sst archive 1981(dec)-2011(july), 
!GSM                                                  ! calling subroutine sstch.f90

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

!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


!GSM      do t=1,356          !dragan, august 2011

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

        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 < 1 .OR. LON2INDX > 4320) WRITE(0,*) "N=",N,"  I=",I,"  J=",J,"  LON1INDX=",LON1INDX,"  LON2INDX=",LON2INDX    ! LUCCI
      IF (LAT1INDX < 1 .OR. LAT2INDX > 2160) WRITE(0,*) "N=",N,"  I=",I,"  J=",J,"  LAT1INDX=",LAT1INDX,"  LAT2INDX=",LAT2INDX    ! LUCCI
!-------------------------
! SOMENTE TESTE DE MBOUNDS
!-------------------------
      IF (LAT2INDX > 2160) LAT2INDX = 2160
      IF (LAT1INDX <    1) LAT1INDX =    1 
        
!        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))         
!***********************************************************************************

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

	      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)
     
      enddo         
      enddo
      enddo
        
!GSM        write(33) msstt               !dragan, august 2011           
!GSM        close(33)            
 !!!!     deallocate(msstt)
!GSM      enddo


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