! source file: /sfs/fs1/work-geomar6/smomw113/UVic_ESCM/2.9_dk_0814/updates/07_aouassim/source/common/sstdata.F
      subroutine sstdata(locali,localj,sstreturn)

!=======================================================================
!     routine to read and interpolate three dimensional SST data for ASE experiments
! wkoeve@geomar.de, 22.06.16, experimental code so far, under development
! This routine shall be called from gasbc.F, similar to cfcdata.f.
!=======================================================================

      implicit none

      character(120) :: fname, name, new_file_name, text

      integer iou, n, ln, nfound,ib(10), ic(10)
      integer locali, localj, c,i,j
      integer :: AllocateStatus, DeAllocateStatus

      logical inqvardef, exists
      logical show_this

      real data_time, wt1, wt3
      real sstreturn

      real, allocatable :: time(:),tmpij(:,:)
      real, allocatable :: datasst(:,:,:), tim(:,:,:)
      real, allocatable :: dat(:,:,:)

!      real datasst(100,100,5000)
!      real tim(100,100,3)
!      real dat(100,100,3)
      save dat, tim, ln, time
      save datasst

!     All variable values stored via save are available each time sstdata()
!     is being entered, kind of internal global variables.

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "cembm.h"
      include "switch.h"
      include "tmngr.h"

      include "csbc.h"

      if (.not. allocated (time)) then
!       variable time has not yet been alocated, hence this section of the code is only entered
!       once during the model run
        name = "O_sstase.nc"
        fname = new_file_name (name)
        inquire (file=trim(fname), exist=exists)
        if (.not. exists) then
          print*, "==> Fatal Error: ", trim(fname), " does not exist."
          print*, "stopping UVic in sstdata.F"
          stop

        else
          call openfile (fname, iou)
          call getdimlen ('MYTIME', iou, ln)
          write(20,*) 'sstdata: dim of Mytime', ln
          if (ln.gt.5000) then
            write(*,*) 'sstdata(): ln is too large, update code'
            stop
          endif
          allocate ( time(ln) )
          allocate ( tmpij(imtm2,jmtm2), STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
          allocate ( datasst(imtm2,jmtm2,ln), STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
          allocate ( tim(imtm2,jmtm2,3), STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
          allocate ( dat(imtm2,jmtm2,3), STAT = AllocateStatus )
          IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"

          datasst(1:imtm2,1:jmtm2,1:ln) = c0
          tmpij(1:imtm2,1:jmtm2) = c0

!         note: the file has been written with ferret and has an offset of 1700 yrs
          call getvara ('MYTIME', iou, ln, 1, ln, time, c1, c0)
          text = 'years'
          call getatttext (iou, 'MYTIME', 'units', text)
          if (trim(text) .eq. "days since 1-1-1")
     &      time(:) = time(:)/yrlen - 1.
          if (trim(text) .eq. "days since 0-1-1")
     &       time(:) = time(:)/yrlen
          if (trim(text) .eq. "years since 1-1-1")
     &      time(:) = time(:) - 1.
          do c=1,ln
             time(c)=time(c)+1700
          enddo

          exists = inqvardef('O_SSTFORC', iou)
          if (.not. exists) then
            print*, "==>  Fatal Error: O_sstforc data does not exist."
            stop
          else
            ib(:) = 1
            ic(:) = 1
            do c=1,ln
             ib(1)=1
             ib(2)=1
             ib(3)=1
             ib(4)=c
             ic(1)=imtm2
             ic(2)=jmtm2
             ic(3)=1
             ic(4)=1
             call getvara ('O_SSTFORC', iou, imtm2*jmtm2, ib, ic
     &,       tmpij, c1, c0)
             do i=1,imtm2
               do j=1,jmtm2
                  datasst(i,j,c) =
     &            tmpij(i,j)
                enddo
              enddo
            enddo

           DEALLOCATE (tmpij, STAT = DeAllocateStatus)

            write(20,*)
     &        'Note: O_sstforc data have been loaded into memory.'
          endif
        endif

        tim(:,:,:) = time(1)

        show_this=.false.
        do j=1,jmtm2
          if (show_this)
     &       write(20,*) 'j = ',j, '  -------------------------'
          do i=1,imtm2
             if (show_this) write(20,*) i, datasst(i,j,1)
             dat(i,j,:) = datasst(i,j,1)
          enddo
        enddo
        if (.not. show_this)
     &    write(20,*) 'c=0 initial probing truned off'

        show_this=.false.
        if (show_this) then
          do c=1,ln
           write(20,*) c, time(c), datasst(50,50,c)
          enddo
        else
          write(20,*) 'i=50,j=50 initial probing turned off'
        endif
        show_this=.false.
!     next line: endif of (.not. allocated (time))
      endif

!     The following code is used each time sstdata.F is called

      data_time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
      tim(locali,localj,2) = min(time(ln), max(time(1), data_time))

!     The following part is the time-interpolation part.
!     For time values outside the range of the sstdata-file, I leave it
!     as in the cfc case. Note, that  this causes the sst-for-ase to be
!     not very well defined outside this range!

       show_this=.false.
       nfound=0.

       if (tim(locali,localj,2) .le. time(1)) then
         dat(locali,localj,2) = datasst(locali,localj,1)

         if  ((show_this).and.(locali.eq.50).and.
     &            (localj.eq.50)) then
            write(20,*) 'tim(2) le time(1) --------'
            write (20,*) time(1), ' ', tim(locali,localj,2)
         endif

       elseif (tim(locali,localj,2) .ge. time(ln)) then
         dat(locali,localj,2) = datasst(locali,localj,ln)
         if  ((show_this).and.(locali.eq.50).and.
     &            (localj.eq.50)) then
            write(20,*) 'tim(2) ge time(ln) --------'
            write (20,*) time(ln), ' ',tim(locali,localj,2)
         endif

       else
         if (tim(locali,localj,2) .gt. tim(locali,localj,3)) then
          do n=2,ln
            if (time(n-1).le.tim(locali,localj,2)
     &         .and. time(n).ge.tim(locali,localj,2)) then
              tim(locali,localj,1) = time(n-1)
              dat(locali,localj,1) = datasst(locali,localj,n-1)
              tim(locali,localj,3) = time(n)
              dat(locali,localj,3) = datasst(locali,localj,n)
              if ((show_this) .and.(locali.eq.50).and.
     &            (localj.eq.50)) then
                write(20,*) ' next time step: '
     &,           locali,' ',localj,' -----------'
              endif
              nfound=n
            endif
            if (nfound.gt.0.) exit
          enddo
         endif
         if ((show_this).and.(locali.eq.50).and.(localj.eq.50)) then
           write(20,*) 'data_t: ', data_time
           write(20,*) 'n:      ', nfound
           write(20,*) 'tim(1): ', tim(locali,localj,1)
           write(20,*) 'tim(2): ', tim(locali,localj,2)
           write(20,*) 'tim(3): ', tim(locali,localj,3)
           write(20,*) 'dat(1): ', dat(locali,localj,1)
           write(20,*) 'dat(2): ', dat(locali,localj,2)
           write(20,*) 'dat(3): ', dat(locali,localj,3)
         endif

        wt1 = 1.
        if (tim(locali,localj,3) .ne. tim(locali,localj,1))
     &      wt1 = (tim(locali,localj,3)-tim(locali,localj,2))
     &         /(tim(locali,localj,3)-tim(locali,localj,1))
        wt1 = max(0., min(1., wt1))
        wt3 = 1. - wt1
        dat(locali,localj,2) = dat(locali,localj,1)*wt1
     &      + dat(locali,localj,3)*wt3
       endif
       sstreturn = dat(locali,localj,2)

       if ((show_this) .and. (locali.eq.50) .and. (localj.eq.50))
     &    write(20,*) tim(locali,localj,2)
     &,      locali, localj, wt1, sstreturn

       show_this=.false.

      return
      end