!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers Universit (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#if defined SOLVE3D && ( defined BHFLUX || defined BWFLUX )
      subroutine get_btflux (itrc)
!
!  Read in point or grided bottom hydrothermal heat flux at the appropriate
!  time from forcing NetCDF file.
!
!
# define BTFLUX_DATA
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "ncscrum.h"
# include "netcdf.inc"
      real cff
      integer itrc, i,ierr, lstr,lvar,lenstr, nf_fread, advance_cycle
!
! Initialization: Inquire about the contents of forcing NetCDF file:
!================  variables and dimensions. Check for consistency.
!
      if (may_day_flag.ne.0) return      !-->  EXIT
      if (itbtf(itrc).eq.0 .or. iic.eq.0) then
        lstr=lenstr(btfname)
c**        call opencdf (btfname,N)
        if (may_day_flag.ne.0) return    !-->  EXIT
!
! If not opened yet, open forcing NetCDF file for reading.
! Find and save IDs for relevant variables, determine whether
! surface flux for tracer itrc is a field or scalar value.
!
        if (ncidbtf.eq.-1) then
          ierr=nf_open (btfname(1:lstr), nf_nowrite, ncidbtf)
          if (ierr .ne. nf_noerr) goto 4
        endif

        if (itrc.eq.itemp) then
          ierr=nf_inq_varid (ncidbtf, 'bhf_time', btf_tid(itrc))
          if (ierr .ne. nf_noerr) then
            write(stdout,3) 'bhf_time', btfname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# ifdef SALINITY
        elseif (itrc.eq.isalt) then
          ierr=nf_inq_varid (ncidbtf, 'bwf_time', btf_tid(itrc))
          if (ierr .ne. nf_noerr) then
            write(stdout,3) 'bwf_time', btfname(1:lstr)
            goto 99                                       !--> ERROR
          endif
# endif
        endif
        lvar=lenstr(vname(1,indxBhflx+itrc-1))
        ierr=nf_inq_varid (ncidbtf, vname(1,indxBhflx+itrc-1),
     &                                          btf_id(itrc))
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidbtf, btf_id(itrc), i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              lbtfgrd(itrc)=1
            else
              lbtfgrd(itrc)=0
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxBhflx+itrc-1), btfname(1:lstr)
          goto 99                                         !--> ERROR
        endif

!
! Determine whether there is cycling to reuse the input data and
! find cycling period "btf_cycle", set initial cycling index
! "btf_ncycle" and record index "btf_rec".
! Set initial value for time index "itbtf" and both time record
! bounds to large negative artificial values, so that it will
! trigger the logic in reading part below.
!
        call set_cycle (ncidbtf,  btf_tid(itrc),    ntbtf(itrc),
     &           btf_cycle(itrc), btf_ncycle(itrc), btf_rec(itrc))
        if (may_day_flag.ne.0) return    !-->  EXIT

        itbtf(itrc)=2
        btf_time(1,itrc)=-1.E+20
        btf_time(2,itrc)=-1.E+20
!
! Set scale factors to convert input fluxes to model units:
!
!     Heat  flux - convert from [Watts/m^2] to [Celsius m/s]
!     Fresh Water flux - convert from [cm/day] to [PSU m/s]
!
        btf_scale(itemp)=1./(rho0*Cp)
# ifdef SALINITY
        btf_scale(isalt)=0.01/86400.
# endif
      endif
!
! Reading data from the forcing file:  Get out, if model time is
!======== ==== ==== === ======= =====  already within the interval
! set by the past and future data times. Otherwise flip the time
! index, increment record and cyclin indices and read a new portion
! of data. Repeat it until model time is between the two times from
! data.
!
  1   i=3-itbtf(itrc)
       cff=time+0.5*dt
        if (btf_time(i,itrc).le.cff .and.
     &      cff.lt.btf_time(itbtf(itrc),itrc)) return
        ierr=advance_cycle (btf_cycle(itrc),  ntbtf(itrc),
     &                      btf_ncycle(itrc), btf_rec(itrc))
        if (ierr .ne. 0) goto 7                           !--> ERROR

        ierr=nf_get_var1_FTYPE(ncidbtf, btf_tid(itrc),
     &                             btf_rec(itrc), cff)
# ifdef USE_CALENDAR
        call tool_origindate(ncidbtf,btf_tid,
     &                       btf_origin_date_in_sec)
        cff=cff+btf_origin_date_in_sec*sec2day
# endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'btf_time', btf_rec(i)
          goto 99                                         !--> ERROR
        endif
        btf_time(i,itrc)=cff*day2sec+btf_cycle(itrc)*btf_ncycle(itrc)
        if (btf_time(itbtf(itrc),itrc).eq.-1.E+20)
     &      btf_time(itbtf(itrc),itrc)=btf_time(i,itrc)

        if (lbtfgrd(itrc).eq.1) then
          ierr=nf_fread(btflxg(START_2D_ARRAY,i,itrc), ncidbtf,
     &                     btf_id(itrc), btf_rec(itrc), r2dvar)
        else
          ierr=nf_get_var1_FTYPE(ncidbtf, btf_id(itrc),
     &                     btf_rec(itrc), btflxp(i,itrc))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'btflux', btf_rec
          goto 99                                         !--> ERROR
        endif
       itbtf(itrc)=i
       MPI_master_only write(stdout,'(6x,A,1x,I2,1x,A,1x,g12.4,1x,I4)')
     &            'GET_BTFLUX -- Read surface flux of tracer', itrc,
     &                                             'for time =', cff
# ifdef USE_CALENDAR
     &  -btf_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                                      , mynode
# endif

      if (ntbtf(itrc).gt.1) goto 1
      if (ntbtf(itrc).eq.1) return

!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_BTFLUX - ERROR: unable to find forcing ',
     &              'variable: ',a,/,14x,'in forcing file: ',a)
  4   write(stdout,5) btfname(1:lstr)
  5   format(/,' GET_BTFLUX - ERROR: unable to open file: ', a)
      goto 99
  6   format(/,' GET_BTFLUX - ERROR while reading variable: ',
     &                                a,' at TIME index = ',i4)
  7   write(stdout,8) btf_rec(itrc),ntbtf(itrc),btfname(1:lstr),
     &                tdays, btf_time(itbtf(itrc),itrc)*sec2day
# ifdef USE_CALENDAR
     &  -btf_origin_date_in_sec*sec2day
# endif
  8   format(/,' GET_BTFLUX - ERROR: requested time record ',I4,
     &       1x,'exeeds the last available', /,14x,'record ',I4,
     &       1x,'in forcing file: ', a,/,14x,'TDAYS = ',
     &             g12.4,2x,'last available BTF_TIME = ',g12.4)
  99  may_day_flag=2
      return
      end

      subroutine set_btflux_tile (Istr,Iend,Jstr,Jend, itrc)
!
!--------------------------------------------------------------------
! Set-up surface tracer flux for current tile.
!--------------------------------------------------------------------
!
# define BTFLUX_DATA
      implicit none
      integer Istr,Iend,Jstr,Jend, itrc, i,j, it1,it2
      real cff, cff1, cff2
# include "param.h"
# include "forces.h"
# include "scalars.h"
!
# include "compute_extended_bounds.h"
!
!      write(*,*)'Enter set_btflux'
      it1=3-itbtf(itrc)
      it2=itbtf(itrc)
      cff=time+0.5*dt
      cff1=btf_time(it2,itrc)-cff
      cff2=cff-btf_time(it1,itrc)
!
! Load time invariant surface tracer flux.
! Time interpolation is not performed in this case.
!
      if (btf_cycle(itrc).lt.0.) then
        if (FIRST_RST_TIME_STEP) then
          if (lbtfgrd(itrc).eq.1) then
            cff=btf_scale(itrc)
            do j=JstrR,JendR
              do i=IstrR,IendR
                btflx(i,j,itrc)=cff*btflxg(i,j,itbtf(itrc),itrc)
!                if ( i==10 .and. j==10) then
!                  write(*,*)'I=',i,' J=',j
!                  write(*,*)'btflx(10,10,1)=',btflx(10,10,1)
!                endif
              enddo
            enddo
          else
            cff=btf_scale(itrc)*btflxp(itbtf(itrc),itrc)
            do j=JstrR,JendR
              do i=IstrR,IendR
                btflx(i,j,itrc)=cff
              enddo
            enddo
          endif
        endif
!
! Time-interpolate surface tracer flux from grided or point data.
! Check that for the next time step [when time=time+dt] time+dt is
! still between srf_tintrp(it1) and srf_tintrp(it2); and if not,
! set synchro_flag top signal that the new forcing data should be
! read from the netCDF input file (master thread only).
!
      elseif (cff1.ge.0. .and. cff2.ge.0.) then
        if (ZEROTH_TILE .and. cff1.lt.dt) synchro_flag=.TRUE.
        cff=btf_scale(itrc)/(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
        if (lbtfgrd(itrc).eq.1) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              btflx(i,j,itrc)=cff1*btflxg(i,j,it1,itrc)
     &                       +cff2*btflxg(i,j,it2,itrc)
            enddo
          enddo
        else
          cff=cff1*btflxp(it1,itrc)+cff2*btflxp(it2,itrc)
          do j=JstrR,JendR
            do i=IstrR,IendR
              btflx(i,j,itrc)=cff
            enddo
          enddo
        endif
        !if ( i==10 .and. j==10) then
        !  write(*,*)'I=',i,' J=',j
        !  write(*,*)'btflx(10,10,1)=',btflx(10,10,1)
        !endif
!
!     Unable to set-up surface tracer flux:
!     Complain about the error and signal to quit (ONE THRED ONLY).
!
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'btf_time',tdays,btf_time(it2,itrc)*sec2day
# ifdef USE_CALENDAR
     &  -btf_origin_date_in_sec*sec2day
# endif
  1       format(/,' SET_BTFLUX_TILE - current model time exceeds',
     &          1x,'ending value for variable: ',A8,/,14x,'TDAYS = ',
     &                                    g12.4, 2x, 'TBTF = ',g12.4)
          may_day_flag=2
        endif
      endif
      return
      end
#else
      subroutine get_btflux_empty
      return
      end
#endif /* defined SOLVE3D && (!defined ANA_BTFLUX || !defined ANA_BSFLUX) */
