subroutine sststp(out) ! taken from regional eta and modified, dragan, april 2012
                             ! used to get sst as initial condition from the monthly sst archive 1981(dec)-2011(july)
                             ! out=sst from monthly sst archive
      include'param_o.h'
      include 'const.h'

      real,dimension(0:im+1,0:jm+1,nm),intent(out):: out
      integer                   :: day
      integer                   :: mon
      integer                   :: yr
      integer                   :: utc
      integer                   :: ma
      integer                   :: mb
      integer                   :: imes
      integer                   :: imes2
      integer                   :: nmessst
      real*4                    :: yday
      real*4                    :: add
      real*4                    :: fa
      real*4                    :: fb
      integer, dimension(12)    :: monl
      character(len=10)         :: chardate

      character(len=256):: m_before, m_after
      character(len=4)::cn

      real,dimension(0:im+1,0:jm+1,nm):: monbefore, monafter

      integer::x,ntsdday,calc_mon,i,j,n

      data monl /31,28,31,30,31,30,31,31,30,31,30,31/

      iyr=idat(3)  !
      imo=idat(1)  ! initial date, day, month, year
      idy=idat(2)  !

      print*,"idy=",idy,"imo=",imo,"iyr=",iyr

      iutc=0        !when used to get initial conditions for the model

      call getdate(iyr,imo,idy,iutc,yr,mon,day,utc)

      if (mod(yr,4) .eq. 0) then
        monl(2)=29
      endif

      yday=day
      ma=mon-1

!      if (yday.ge.(float(monl(mon))/2.0)) then
      if (yday.ge.(1.0+float(monl(mon))/2.0)) then  ! changing the month of sst from the archive,
                                                    ! each 16th of the month, for the months with 30 days,
                                                    ! each 17th of the month, for the months with 31 days,
                                                    ! 15th of the month each february,
                                                    ! except for the case when it has 29 days, then it changes at the day 16th

        ma=mon
      endif

      mb=ma+1

      if (ma .lt. 1) then
        ma=12
      endif

      if (mb .gt. 12)then
        mb=1
      endif


      print*,"ma",ma,"day",yday, "yr",yr


!=====================================================
!+++ new part of the code - dragan, september 2011 +++
!=====================================================



                                                   ! the archive consists of the monthly sst data for the period december 1981 - july 2011

      initsst=(iyr-1982)*12+1+imo               ! this line calculates the initial sst
                                                   ! if the day is in the first half of each month, it will be recalculated below, in calc_mon,
                                                   ! it will use ma=mon-1, and calc_mon will start from -1


      if (ma.eq.12.and.yday.lt.(float(monl(mon))/2.0+1.0)) then       ! defines the problem when the year is changing, for the first 15 days of each january
         calc_mon=(yr-iyr)*12+(ma-imo)-12                        ! calculates the current month of sst taken from archive, for the period 1-15th of january each year
      else
         calc_mon=(yr-iyr)*12+(ma-imo)                           ! calculates the current month of sst taken from archive for all the other dates
      endif

      imes=initsst+calc_mon                      ! "previous" month from the archive, for making interpolation of monthly sst for specific date
      imes2=imes+1                               ! "future" month from the archive, for making interpolation of monthly sst for specific date


      if (imes.lt.1) then
        imes=12
      endif
!in the case that the initial conditions are from the period 1-15th of january, and if yr=iyr
      if (imes2.gt.12)then
        mb=1
      endif


!++++++++++++++++++++++++++++++++++++++++

!       ntsdday=86400/dt

!     if (idy.ge.(1.0+float(monl(imo))/2.0)) then    ! a part of the code added to fix the old version used in regional eta, in the case that each month has 30 days
!       x=(float(monl(imo))-idy)*ntsdday
!     else
!       x=-(float(monl(imo))-idy)*ntsdday
!     endif




!+++++++++++++++++++++++++++++++++++++++++
!         initsst=297
!         ntsd1=(ntsd-x+(nmessst/2))

!        imes=int(ntsd1/nmessst)+initsst             ! old part of the code taken from regional eta, which will work correctly only in the case
!        imes2=imes+1                                ! that each month has 30 days, but together with the definition of x, declared in the code above

!        print*,'x=',x
!++++++++++++++++++++++++++++++++++++++++++

       add=float(monl(ma))/2.0-1.0

       if (ma.eq.mon) then
         add=-add-2.
       endif

       fb=2.0*(yday+add)/float(monl(ma)+monl(mb))
       fa=1.0-fb


       print*,"imes=",imes,"imes2=",imes2
       print*,"fa=",fa,"fb=",fb

       write(cn,'(i4.4)')imes
       m_before='/scratchout/grupos/grpeta/projetos/tempo/oper/gef_v1.0.0/initdata/monthly_sst'//cn//'.bin'
       open(11,file=m_before,form='unformatted',status='old',access='sequential')
       read(11) monbefore
       close(11)

       write(cn,'(i4.4)')imes2
       m_after='/scratchout/grupos/grpeta/projetos/tempo/oper/gef_v1.0.0/initdata/monthly_sst'//cn//'.bin'
       open(22,file=m_after,form='unformatted',status='old',access='sequential')
       read(22) monafter
       close(22)

       out(:,:,:)=fa*monbefore(:,:,:)+fb*monafter(:,:,:)

       return
       end