#include "options.inc"

c
c-----------------------------------------------------------------------
c     Module which handles IO for forcing functions
c-----------------------------------------------------------------------
c

      module obc_io_module
      implicit none
c
c-----------------------------------------------------------------------
c     read in and time interpolate open boundary conditions
c     from a netcdf or binary file, which have to be prepared
c     beforehand. 
c
c     linked in the code in subr. setup and tracer only
c
c                                       c.eden
c-----------------------------------------------------------------------
c
      type :: obc_type
        integer index         ! the index for the time interpolator
        logical active     
        logical ts_active     
        logical psi_active     
        real, pointer :: ts_cache(:,:,:,:) ! cache for 2 time steps 
                                           ! of T/S and time scale
        integer cache_pointer(2)   ! pointer for the cache
        real, pointer :: psi_cache(:,:) ! cache for 2 time steps 
        integer length             ! number of time steps of the data
        character (len=80) :: file
        character (len=1) :: orient
        integer :: fdim
        integer :: ntr ! number of actual tracers in input files
      end type obc_type

      type (obc_type) :: obc(4)
      logical  :: netcdf_obc_verbose = .true.

      end module obc_io_module


      subroutine obc_initialize
c-----------------------------------------------------------------------
c     initialize the open boundary time interpolation
c-----------------------------------------------------------------------
      use spflame_module
      use obc_io_module
      implicit none
      integer :: n,io,k
      type( time_type) :: start
      type( time_type), allocatable :: interval(:)
      real, allocatable :: aprec(:)
      character (len=32) :: stamp

      if (my_pe==0) print*,' '
      if (my_pe==0) 
     &   print*,' Initialising of open boundary forcing data '
      if (my_pe==0) print*,' '

      if (my_pe==0 ) then 
       netcdf_obc_verbose=.true.
      else
       netcdf_obc_verbose=.false.
      endif

      obc(1)%orient='n'
      obc(2)%orient='s'
      obc(3)%orient='e'
      obc(4)%orient='w'
      obc(1)%file='obc_north.mom'
      obc(2)%file='obc_south.mom'
      obc(3)%file='obc_east.mom'
      obc(4)%file='obc_west.mom'

      obc(:)%active=.false.
      obc(:)%ts_active=.false.
      obc(:)%psi_active=.false.

      n=1
      if (enable_obc_north.and.my_blk_j==n_pes_j) then
        obc(n)%active=.true.
        if (restore_TS_obc_north) obc(n)%ts_active=.true.
        if (prescribe_psi_obc_north) obc(n)%psi_active=.true.
        obc(n)%fdim=imt
      endif

      n=2
      if (enable_obc_south.and.my_blk_j==1) then
        obc(n)%active=.true.
        if (restore_TS_obc_south) obc(n)%ts_active=.true.
        if (prescribe_psi_obc_south) obc(n)%psi_active=.true.
        obc(n)%fdim=imt
      endif

      n=3
      if (enable_obc_east.and.my_blk_i==n_pes_i) then
        obc(n)%active=.true.
        if (restore_TS_obc_east) obc(n)%ts_active=.true.
        if (prescribe_psi_obc_east) obc(n)%psi_active=.true.
        obc(n)%fdim=jmt
      endif

      n=4
      if (enable_obc_west.and.my_blk_i==1) then
        obc(n)%active=.true.
        if (restore_TS_obc_west) obc(n)%ts_active=.true.
        if (prescribe_psi_obc_west) obc(n)%psi_active=.true.
        obc(n)%fdim=jmt
      endif

      do n=1,4
       if (obc(n)%active) then
        if (obc(n)%ts_active .or. obc(n)%psi_active ) then
         if (my_pe==0) print*,' activating OBC ',obc(n)%orient
         allocate( obc(n)%ts_cache(obc(n)%fdim,km,nt,2) )
         obc(n)%ts_cache(obc(n)%fdim,km,nt,2) = 0.
         allocate( obc(n)%psi_cache(obc(n)%fdim,2) )
         obc(n)%psi_cache(obc(n)%fdim,2) = 0.
         obc(n)%cache_pointer=-1 ! this is important
         if (my_pe==0) print*,' reading file :',obc(n)%file
         call getunit(io,obc(n)%file,'usr ieee')
         read (io) obc(n)%length
         if (my_pe==0) print*,' time steps :',obc(n)%length
         allocate(interval(obc(n)%length))
         allocate(aprec(obc(n)%length))
         rewind(io)
         read (io)  k,stamp,aprec
         if (my_pe==0) print*,' starting at ',stamp
         start= get_stamp(stamp)
         do k=1,obc(n)%length
           call set_time(interval(k),
     &               int((aprec(k)-int(aprec(k)))*24*60*60),
     &               int(aprec(k)))
         enddo
         call init_forcing_interp(obc(n)%length,.true.,
     &                            start,interval,obc(n)%index)
         obc(n)%ntr=nt
 20      rewind(io)
         read(io)
         read (io,err=30) obc(n)%ts_cache(:,:,1:obc(n)%ntr,1), 
     &                    obc(n)%psi_cache(:,1) 
         goto 40
 30      continue
         if (obc(n)%ntr > 2 ) then
          obc(n)%ntr=obc(n)%ntr-1
          goto 20
         else
          if (my_pe==0) print*,' cannot read input file'
          call halt_stop(' in OBC initialization' )
         endif
 40      close(io)
         if (my_pe==0) print*,' found ',obc(n)%ntr,' tracers in files'
         deallocate(interval,aprec)
        endif
       endif
      enddo


      if (my_pe==0) print*,' '
      if (my_pe==0) 
     &   print*,' Initialisation of open boundary forcing data done'
      if (my_pe==0) print*,' '

      end subroutine obc_initialize



      subroutine obc_read
      use spflame_module
      use obc_io_module
      implicit none
      integer :: k,p1,p2,kk,i,j
      real    :: f1,f2
      logical load(2)
      character(len=80) :: text
      integer :: io,n
      real :: dummy
      integer is,ie,js,je
#ifdef SR8000_host
      real, allocatable :: buf(:,:,:), buf2(:)
#endif

      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      do k=1,4
       if (obc(k)%psi_active.or.obc(k)%ts_active) then

        call forcing_interp(obc(k)%index,p1,p2,f1,f2)
        write(text,'(" obc data orientation : ",a1)') obc(k)%orient
        call update_forcing_cache(obc(k)%index,p1,p2,
     &                obc(k)%cache_pointer(:),load,text)
        do kk=1,2
         if (load(kk)) then

#ifdef SR8000_host
c there is apparently a bug in fortran types and precexpio option
c for the compiler on the SR8000: use the following work around

           if (obc(k)%ts_active.or.obc(k)%psi_active) then
            call getunit(io,obc(k)%file,'usr ieee')
            allocate( buf(obc(k)%fdim,km,obc(k)%ntr) )
            allocate( buf2(obc(k)%fdim) )
            read (io)  ! length
            do n=1,obc(k)%cache_pointer(kk)-1
             read(io)
            enddo
            read(io)
     &         (((buf(i,j,n),i=1,obc(k)%fdim),j=1,km),n=1,obc(k)%ntr),
     &               (buf2(i),i=1,obc(k)%fdim)
            if (obc(k)%ts_active) then
              obc(k)%ts_cache(:,:,1:obc(k)%ntr,kk) =buf
            endif
            if (obc(k)%psi_active) then
              obc(k)%psi_cache(1:obc(k)%fdim,kk) =buf2(1:obc(k)%fdim)
            endif
            close(io)
            deallocate(buf, buf2)
           endif

#else
           call getunit(io,obc(k)%file,'usr ieee')
           if (obc(k)%ts_active.or.obc(k)%psi_active) then
            read (io)  ! length
            do n=1,obc(k)%cache_pointer(kk)-1
             read(io)
            enddo
            if (obc(k)%ts_active.and.obc(k)%psi_active) then
             read (io)obc(k)%ts_cache(:,:,1:obc(k)%ntr,kk), 
     &                obc(k)%psi_cache(:,kk) 
            elseif (obc(k)%ts_active) then
             read (io)obc(k)%ts_cache(:,:,1:obc(k)%ntr,kk) 
            elseif (obc(k)%psi_active) then
             read (io)(dummy,i=1,obc(k)%fdim*km*obc(k)%ntr),
     &                 obc(k)%psi_cache(:,kk) 
            endif
           endif
           close(io)
#endif

         endif ! load(kk)
        enddo

        if (obc(k)%ts_active.and.obc(k)%orient=='n')
     &   ts_obc_north(is:ie,:,:)=obc(k)%ts_cache(is:ie,:,:,p1)*f1+
     &                           obc(k)%ts_cache(is:ie,:,:,p2)*f2
        if (obc(k)%ts_active.and.obc(k)%orient=='s')
     &   ts_obc_south(is:ie,:,:)=obc(k)%ts_cache(is:ie,:,:,p1)*f1+
     &                           obc(k)%ts_cache(is:ie,:,:,p2)*f2
        if (obc(k)%ts_active.and.obc(k)%orient=='w')
     &   ts_obc_west(js:je,:,:)=obc(k)%ts_cache(js:je,:,:,p1)*f1+
     &                          obc(k)%ts_cache(js:je,:,:,p2)*f2
        if (obc(k)%ts_active.and.obc(k)%orient=='e')
     &   ts_obc_east(js:je,:,:)=obc(k)%ts_cache(js:je,:,:,p1)*f1+
     &                          obc(k)%ts_cache(js:je,:,:,p2)*f2

        if (obc(k)%psi_active.and.obc(k)%orient=='n') 
     &   psi_wall_north(is:ie) = obc(k)%psi_cache(is:ie,p1)*f1+
     &                           obc(k)%psi_cache(is:ie,p2)*f2
        if (obc(k)%psi_active.and.obc(k)%orient=='s') 
     &   psi_wall_south(is:ie) = obc(k)%psi_cache(is:ie,p1)*f1+
     &                           obc(k)%psi_cache(is:ie,p2)*f2
        if (obc(k)%psi_active.and.obc(k)%orient=='w')
     &   psi_wall_west(js:je) = obc(k)%psi_cache(js:je,p1)*f1+
     &                          obc(k)%psi_cache(js:je,p2)*f2
        if (obc(k)%psi_active.and.obc(k)%orient=='e') 
     &   psi_wall_east(js:je) = obc(k)%psi_cache(js:je,p1)*f1+
     &                          obc(k)%psi_cache(js:je,p2)*f2

       endif
      enddo
      end subroutine obc_read



      module sponge_module
c
c-----------------------------------------------------------------------c
c     read in and time interpolate data for sponge layers
c     from a netcdf or binary file, which have to be prepared
c     beforehand. 
c
c     linked in the code in subr. setup and tracer only
c
c                                       c.eden
c-----------------------------------------------------------------------
c
      implicit none

      type :: spg_type
        integer is,ie,js,je   ! start and end indicees for the sponge layer
        integer isb,ieb,jsb,jeb ! same but restricted to the PE domain
        integer index         ! the index for the time interpolator
        logical active        ! true if sponge layer is in PE domain
        real, pointer :: cache(:,:,:,:,:) ! cache for 2 time steps 
                                          ! of T/S and time scale
        integer cache_pointer(2)   ! pointer for the cache
        integer length             ! number of time steps of the data
        integer ntr                ! number of tracer in input file
      end type spg_type
      integer :: number_of_spg ! total number of sponge layers
                                       ! use not more than 99
      type(spg_type), allocatable :: spg(:) ! the sponge layers
      logical :: netcdf_spg_verbose = .true. ! be verbose
      end module sponge_module


      subroutine sponge_initialize
c
c-----------------------------------------------------------------------
c     initialize the sponge layer time interpolation
c-----------------------------------------------------------------------
c
      use spflame_module
      use sponge_module
      implicit none
      integer k,n,io,i
      character (len=80) :: name
      character (len=32) :: stamp
      type( time_type) :: start
      type( time_type), allocatable :: interval(:)
      real, allocatable :: aprec(:)
      real, allocatable :: trbuf(:,:,:)

      if (my_pe==0) print*,' '
      if (my_pe==0) 
     &   print*,' Initialising of sponge layer data '
      if (my_pe==0) print*,' '

      netcdf_spg_verbose=.false.
      if (my_pe==0 ) netcdf_spg_verbose=.true.

      if (my_pe==0) then
       print*,' reading sponge informations from file sponge_01.mom'
      endif
      call getunit(io,'sponge_01.mom','usr ieee')
      read (io) number_of_spg
      close(io)

      if (my_pe==0) then
        print*,' found ',number_of_spg,' sponge layers'
      endif

      allocate(spg(number_of_spg));

      do n=1,number_of_spg
       if (my_pe==0) print*,' initializing layer #',n
c      read time dimension and horizontal location for sponge
c      layer number n
       write(name, '("sponge_",i2,".mom")') n
       do i=1,len_trim(name)
         if (name(i:i) == ' ') name(i:i)='0'
       enddo
       if (my_pe==0) print*,' reading from file ',name(1:len_trim(name))
       call getunit(io,name,'usr ieee')
       read (io) i,spg(n)%length,
     &             spg(n)%is,spg(n)%ie,spg(n)%js,spg(n)%je
       if (my_pe==0) print*,' time steps :',spg(n)%length
       if (my_pe==0) print*,' range of sponge :',
     &             spg(n)%is,spg(n)%ie,spg(n)%js,spg(n)%je
c      see how many tracers are actually in file
       if (my_pe==0) then
 
        allocate(trbuf(spg(n)%is:spg(n)%ie,spg(n)%js:spg(n)%je,nt*2))
        spg(n)%ntr = nt
20      rewind(io)
         read(io)
         do i=1,km
           read(io, err=30) trbuf(:,:,1:spg(n)%ntr*2)
         enddo
         goto 40
 30      continue
         if (spg(n)%ntr > 2 ) then
          spg(n)%ntr=spg(n)%ntr-1
          goto 20
         else
          if (my_pe==0) print*,' cannot read input file'
          call halt_stop(' in sponge initialization' )
         endif
 40     continue
        deallocate( trbuf )
       endif ! my_pe == 0

       call bcast_integer(spg(n)%ntr,1,0)
       rewind(io)

       if (my_pe==0) print*,' found ',spg(n)%ntr,' tracers in file' 
       call sub_flush(6)
c
c      restrict array to PE domain
       spg(n)%jsb=max(js_pe,spg(n)%js)
       spg(n)%jeb=min(je_pe,spg(n)%je)
       spg(n)%isb=max(is_pe,spg(n)%is)
       spg(n)%ieb=min(ie_pe,spg(n)%ie)

       if (spg(n)%jeb >= spg(n)%jsb .and.
     &     spg(n)%ieb >= spg(n)%isb ) then

        spg(n)%active=.true.
        allocate( spg(n)%cache(spg(n)%isb:spg(n)%ieb,km,
     &                   spg(n)%jsb:spg(n)%jeb,spg(n)%ntr*2,2) )
        spg(n)%cache=0.
        spg(n)%cache_pointer=-1

        allocate(interval(spg(n)%length))
        allocate(aprec(spg(n)%length))
c       read the time variables for sponge layer # n
        read (io)  k,k,i,i,i,i,stamp,aprec
        start= get_stamp (stamp)
        do k=1,spg(n)%length
         call set_time(interval(k),
     &                 int((aprec(k)-int(aprec(k)))*24*60*60),
     &                 int(aprec(k)))
        enddo
        call init_forcing_interp(spg(n)%length,.true.,
     &                           start,interval,spg(n)%index)
        deallocate(interval,aprec)
       else
        spg(n)%active=.false.
       endif
       close(io)
       call sub_flush(6)
      enddo

      if (my_pe==0) print*,' '
      if (my_pe==0) 
     &   print*,' Initialisation of sponge layer data done'
      if (my_pe==0) print*,' '
       call sub_flush(6)

      end subroutine sponge_initialize




      subroutine sponge_read(n)
c
c-----------------------------------------------------------------------
c     read data from netcdf file and compute source term
c     for tracer number n
c-----------------------------------------------------------------------
c
      use spflame_module
      use sponge_module
      implicit none
      integer, intent(in) :: n ! number of the tracer
      integer :: k,p1,p2,kk,i,j,m
      real    :: f1,f2,tstar,tscl
      logical load(2)
      character(len=80) :: text
      integer :: io
      character (len=80) :: name
      real, allocatable :: trbuf(:,:,:),tsbuf(:,:,:)
      character(len=2) :: number
      
      do k=1,number_of_spg
       write(number,'(i2)') k
       if (number(1:1) == ' ') number(1:1)='0'
       if (spg(k)%active) then

        call forcing_interp(spg(k)%index,p1,p2,f1,f2)
        write(text,'(" sponge layer #",i3)') k
        call update_forcing_cache(spg(k)%index,p1,p2,
     &                spg(k)%cache_pointer(:),load,text)
        do kk=1,2
         if (load(kk)) then
           write(name, '("sponge_",i2,".mom")') k
           do i=1,len_trim(name)
            if (name(i:i) == ' ') name(i:i)='0'
           enddo
           call getunit(io,name,'usr ieee')
           read (io)
           do i=1,(spg(k)%cache_pointer(kk)-1)*km
            read(io)
           enddo
           allocate(trbuf(spg(k)%is:spg(k)%ie,spg(k)%js:spg(k)%je,
     &              spg(k)%ntr))
           allocate(tsbuf(spg(k)%is:spg(k)%ie,spg(k)%js:spg(k)%je,
     &              spg(k)%ntr))
           do i=1,km
            read(io) trbuf,tsbuf
            do m=1,spg(k)%ntr
             spg(k)%cache(:,i,:,m,kk)=trbuf(spg(k)%isb:spg(k)%ieb,
     &                                 spg(k)%jsb:spg(k)%jeb,m)
             spg(k)%cache(:,i,:,spg(k)%ntr+m,kk)=
     &        tsbuf(spg(k)%isb:spg(k)%ieb,spg(k)%jsb:spg(k)%jeb,m)
            enddo
           enddo
           close(io)
           deallocate(trbuf,tsbuf) 
         endif
        enddo

        if (n<= spg(k)%ntr ) then

         do j=spg(k)%jsb,spg(k)%jeb
          do kk=1,km
           do i=spg(k)%isb,spg(k)%ieb
            tstar= spg(k)%cache(i,kk,j,n,p1)*f1+
     &             spg(k)%cache(i,kk,j,n,p2)*f2
            tscl = spg(k)%cache(i,kk,j,spg(k)%ntr+n,p1)*f1+
     &             spg(k)%cache(i,kk,j,spg(k)%ntr+n,p2)*f2
c           restore to tracer at tau-1
            source(i,kk,j)=source(i,kk,j)+
     &        tscl*(tstar-t(i,kk,j,n,taum1))*tmask(i,kk,j)
           enddo
          enddo
         enddo
        endif

       endif
      enddo
      end subroutine sponge_read


      module sbc_io_module
c
c-----------------------------------------------------------------------c
c     read in and time interpolate surface boundary conditions
c     from a netcdf or binary file, which have to be prepared
c     beforehand. 
c
c     linked in the code in subr. setup and vmixc only
c
c                                       c.eden
c-----------------------------------------------------------------------
c
      implicit none
      integer                     :: number_of_sbc
      integer                     :: number_of_ustar
      integer                     :: number_of_qsol
      integer,allocatable         :: sbc_index(:)
      real,allocatable            :: sbc_cache(:,:,:,:)
      integer,allocatable         :: sbc_cache_pointer(:,:) 
      integer,allocatable         :: sbc_length(:) 
      character(len=32),allocatable         :: sbc_name(:) 
      character(len=80),allocatable         :: filename(:) 
      end module sbc_io_module


      subroutine sbc_initialize
c-----------------------------------------------------------------------
c     initialize the SBC time interpolation
c-----------------------------------------------------------------------
      use spflame_module
      use sbc_io_module
      implicit none
      integer k,n,io,m
      type( time_type) :: start
      type( time_type), allocatable :: interval(:)
      real, allocatable :: aprec(:)
      character (len=32) :: stamp
      character (len=80) :: name
      real :: pen(0:km),swarg1,swarg2
      
      if (my_pe==0) then
       print*,' '
       print*,' Initialising of surface boundary forcing data '
       print*,' '
       if  (enable_salt_flux_sbc) then
        print*,' using a fixed salt flux b.c. '
       endif
       if (enable_shortwave_sbc) then
        print*,' using penetrative solar radiation '
       endif
      endif
c
c     Number of files to be read
c
      number_of_sbc = 2+nt*2  ! normal setup: taux,tauy,sst/sss_clim/rest

      if (enable_ktmix.or.enable_tkemix) then
        number_of_sbc   = number_of_sbc+1
        number_of_ustar = number_of_sbc
      endif

      if (enable_shortwave_sbc) then
        number_of_sbc   = number_of_sbc+1
        number_of_qsol  = number_of_sbc
c
c     Initialize penetration profile for solar radiation
c     and store divergence in divpen
c     note that pen(0) is set 0.0 instead of 1.0 to compensate for the
c     shortwave part of the total surface flux in "stf(i,1)"
c
        pen(0) = 0.0
        do k=1,km
         swarg1 = -min(zw(k)/efold1_shortwave,70.0)
         swarg2 = -min(zw(k)/efold2_shortwave,70.0)
         pen(k) = rpart_shortwave*exp(swarg1) 
     &         + (1.0-rpart_shortwave)*exp(swarg2)
  	 divpen_shortwave(k) = (pen(k-1) - pen(k))/dzt(k)
        enddo
        if (my_pe==0) then
          print*,' A part of ',rpart_shortwave,' of the solar radiation'
          print*,' penetrates exponentially with a length scale of'
          print*,efold1_shortwave,' cm the rest with a length scale of'
          print*,efold2_shortwave,' cm'
          print*,' Resulting exponential penetration profile: '
          do k=0,km
            print*,' k=',k,' zw(k)=', zw(k),' pen(k)=', pen(k)
          enddo
        endif
      endif

      allocate(sbc_index(number_of_sbc));sbc_index=9999999
      allocate(sbc_cache(number_of_sbc,is_pe:ie_pe,js_pe:je_pe,2))
      sbc_cache=0.
      allocate(sbc_cache_pointer(number_of_sbc,2))
      sbc_cache_pointer = -1
      allocate(sbc_length(number_of_sbc)); sbc_length=0
      allocate(sbc_name(number_of_sbc)); sbc_name='noname'
      allocate(filename(number_of_sbc)); filename='noname.cdf'

      sbc_name(1)='taux'
      sbc_name(2)='tauy'
      sbc_name(3)='sst_clim'
      sbc_name(4)='sst_rest'
      sbc_name(5)='sss_clim'
      if  (enable_salt_flux_sbc) then
       sbc_name(6)='sss_flux'
      else
       sbc_name(6)='sss_rest'
      endif
      if (enable_ktmix.or.enable_tkemix) then
          sbc_name(number_of_ustar)='ustar'
      endif
      if (enable_shortwave_sbc) then
          sbc_name(number_of_qsol)='qsol'
      endif

      if (my_pe==0) print*,' using binary files '
      filename(1)='taux.mom'
      filename(2)='tauy.mom'
      filename(3)='sst_clim.mom'
      filename(4)='sst_rest.mom'
      filename(5)='sss_clim.mom'
      if  (enable_salt_flux_sbc) then
       filename(6)='sss_flux.mom'
      else
       filename(6)='sss_rest.mom'
      endif
      if (enable_ktmix.or.enable_tkemix) then
         filename(number_of_ustar)='ustar.mom'
      endif
      if (enable_shortwave_sbc) then
         filename(number_of_qsol)='qsol.mom'
      endif

      do m=3,nt

       write(name,'("str_",i2,"_clim")') m
       call replace_space_zero(name)
       sbc_name(6+(m-3)*2+1)=name(1:32)
       write(name,'("str_",i2,"_clim.mom")') m
       call replace_space_zero(name)
       filename(6+(m-3)*2+1)=name

       write(name,'("str_",i2,"_rest")') m
       call replace_space_zero(name)
       sbc_name(6+(m-3)*2+2)=name(1:32)
       write(name,'("str_",i2,"_rest.mom")') m
       call replace_space_zero(name)
       filename(6+(m-3)*2+2)=name

      enddo

      do n=1,number_of_sbc
       call getunit(io,filename(n),'usr ieee')
       read (io) sbc_length(n)
       rewind(io)
       allocate(interval(sbc_length(n)))
       allocate(aprec(sbc_length(n)))
       read (io)  k,stamp,aprec
       close(io)
       start= get_stamp (stamp)
       do k=1,sbc_length(n)
        call set_time(interval(k),
     &               int((aprec(k)-int(aprec(k)))*24*60*60),
     &               int(aprec(k)))
       enddo
       call init_forcing_interp(sbc_length(n),.true.,
     &                          start,interval,sbc_index(n))
       deallocate(interval,aprec)
      enddo

      if (my_pe==0) print*,' '
      if (my_pe==0) 
     &   print*,' Initialising of surface boundary forcing data done'
      if (my_pe==0) print*,' '
      end subroutine sbc_initialize


      subroutine sbc_read
c
c-----------------------------------------------------------------------
c     interpolate in time between adjacent months,
c     read in new data from netcdf files into cache if necessary.
c-----------------------------------------------------------------------
c
      use spflame_module
      use sbc_io_module
      implicit none
      real :: buf(imt,jmt),f1,f2
      integer :: p1,p2,i,j,k,kk,io,m
      logical load(2)

      do k=1,number_of_sbc
       call forcing_interp(sbc_index(k),p1,p2,f1,f2)
       call update_forcing_cache(sbc_index(k),p1,p2,
     &                sbc_cache_pointer(k,:),load,sbc_name(k))
       do kk=1,2
        if (load(kk)) then
         if (my_pe==0) then 
          call getunit(io,filename(k),'usr ieee')
          read (io) 
          do i=1,sbc_cache_pointer(k,kk)-1 
           read (io)  
          enddo
          read (io) buf
          close(io)
         endif
         call bcast_real(buf,imt*jmt,0)
         sbc_cache(k,:,:,kk)=buf(is_pe:ie_pe,js_pe:je_pe)
        endif
       enddo

       if (k==1.or.k==2) then

        smf(is_pe:ie_pe,js_pe:je_pe,k) = 
     &     sbc_cache(k,:,:,p1)*f1+sbc_cache(k,:,:,p2)*f2

       elseif ( k <= nt*2+2 ) then

        m=floor((k-3.)/2.) +1  ! number of tracer for this tracer

        if (ceiling( (k-3.)/2.+1 )==m) then
         stf_clim(is_pe:ie_pe,js_pe:je_pe,m) = 
     &     sbc_cache(k,:,:,p1)*f1+sbc_cache(k,:,:,p2)*f2

         if  (enable_salt_flux_sbc) then
c
c        here we are using a fixed salt flux
c        sft = rest*(clim-t(...taum1)
c          with clim = t(...,taum1)+1 and rest = flux
c        stf = flux*(t(...,taum1)+1-t(...taum1)= flux
c
          if (m==2)  stf_clim(is_pe:ie_pe,js_pe:je_pe,m)=
     &      t(is_pe:ie_pe,1,js_pe:je_pe,2,taum1)+1.
c
c        note that this could be also made nicer with
c        a logical :: restoring_sbc(1:nt)
c
         endif
        else
         stf_rest(is_pe:ie_pe,js_pe:je_pe,m) = 
     &     sbc_cache(k,:,:,p1)*f1+sbc_cache(k,:,:,p2)*f2
        endif


       elseif (k==number_of_ustar .and. 
     &          (enable_ktmix.or.enable_tkemix) ) then
        ustar(is_pe:ie_pe,js_pe:je_pe) = 
     &     sbc_cache(k,:,:,p1)*f1+sbc_cache(k,:,:,p2)*f2
       elseif (k==number_of_qsol .and. enable_shortwave_sbc)  then
        qsol(is_pe:ie_pe,js_pe:je_pe) = 
     &     sbc_cache(k,:,:,p1)*f1+sbc_cache(k,:,:,p2)*f2
       else
          print*,' error in sbc_read '
          call halt_stop(' in sbc_read')
       endif
      enddo

      end subroutine sbc_read

