!======================================================================
! 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 University (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 QCORRECTION && !defined ANA_SST  && defined TEMPERATURE


      subroutine get_sst
!
!  Read in point or grided sea surface temperature surface net heat
!  flux sensitivity to sea surface temperature at the appropriate
!  time from forcing NetCDF file.
!
!  These forcing fields are used when flux correction is activated:
!
!       Q_model ~ Q + dQdSST * (T_model - SST)
!
# ifdef OPENACC
      use exchange_device_host
# endif
# define SST_DATA
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "netcdf.inc"
# include "ncscrum.h"
      real cff
      integer 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 (itsst.eq.0 .or. iic.eq.0) then
        lstr=lenstr(frcname)
c*        call opencdf (frcname,N)
c*        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
! SST is a field or scalar value.
!
        if (ncidfrc.eq.-1) then
          ierr=nf_open(frcname(1:lstr), nf_nowrite, ncidfrc)
          if (ierr. ne. nf_noerr) goto 4                  !--> ERROR
        endif

        ierr=nf_inq_varid (ncidfrc, 'sst_time', sst_tid)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) 'sst_time', frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif

        lvar=lenstr(vname(1,indxSST))
        ierr=nf_inq_varid (ncidfrc, vname(1,indxSST)(1:lvar), sst_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, sst_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              lsstgrd=1
            else
              lsstgrd=0
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxSST)(1:lvar), frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif

        lvar=lenstr(vname(1,indxdQdSST))
        ierr=nf_inq_varid (ncidfrc, vname(1,indxdQdSST)(1:lvar),
     &                                                    dqdt_id)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxdQdSST)(1:lvar),frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
! Determine whether there is cycling to reuse the input data and
! find cycling period "sst_cycle", set initial cycling index
! "sst_ncycle" and record index "sst_rec".
! Set initial value for time index "itsst" and both time record
! bounds to large negative artificial values, so that it will
! trigger the logic in reading part below.
! Also set scale factor to convert input dQdSST from Watts/m2/Celsius
! to meter/second.
!
        call set_cycle (ncidfrc,   sst_tid,    ntsst,
     &                  sst_cycle, sst_ncycle, sst_rec)
        if (may_day_flag.ne.0) return    !-->  EXIT
        itsst=2
        sst_time(1)=-1.E+20
        sst_time(2)=-1.E+20
        scldqdt=1./(rho0*Cp)
      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-itsst
       cff=time+0.5*dt
        if (sst_time(i).le.cff .and. cff.lt.sst_time(itsst)) return
        ierr=advance_cycle (sst_cycle, ntsst, sst_ncycle, sst_rec)
        if (ierr .ne. 0) goto 7                           !--> ERROR

        ierr=nf_get_var1_FTYPE (ncidfrc, sst_tid, sst_rec, cff)
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'sst_time', sst_rec
          goto 99                                         !--> ERROR
        endif
# ifdef USE_CALENDAR
        call tool_origindate(ncidfrc,sst_tid,
     &                       sst_origin_date_in_sec)
        cff=cff+sst_origin_date_in_sec*sec2day
# endif

        sst_time(i)=cff*day2sec+sst_cycle*sst_ncycle
        if (sst_time(itsst).eq.-1.E+20) sst_time(itsst)=sst_time(i)

        if (lsstgrd.eq.1) then
          ierr=nf_fread (sstg(START_2D_ARRAY,i), ncidfrc, sst_id,
     &                                           sst_rec, r2dvar)
# ifdef OPENACC
          call update_device(sstg(:,:,i))
# endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc, sst_id, sst_rec, sstp(i))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'SST', sst_rec
          goto 99                                         !--> ERROR
        endif

        if (lsstgrd.eq.1) then
          ierr=nf_fread (dqdtg(START_2D_ARRAY,i), ncidfrc, dqdt_id,
     &                                            sst_rec,  r2dvar)
# ifdef OPENACC
          call update_device(dqdtg(:,:,i))
# endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc,dqdt_id,sst_rec,dqdtp(i))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'dQdSST', sst_rec
          goto 99                                         !--> ERROR
        endif
       itsst=i
       MPI_master_only write(stdout,'(6x,A,1x,A,1x,g12.4,1x,I4)')
     &              'GET_SST    --',
     &              'Read SST and dQdSST fields     for time =', cff
# ifdef USE_CALENDAR
     &  -sst_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                                      , mynode
# endif
      if (ntsst.gt.1) goto 1
      if (ntsst.eq.1) return
!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_SST - ERROR: unable to find forcing variable',
     &                   ': ',a,/,11x,'in forcing NetCDF file: ',a)
  4   write(stdout,5) frcname(1:lstr)
  5   format(/,' GET_SST - ERROR: unable to open forcing NetCDF ',
     &                                                  'file: ',a)
      goto 99
  6   format(/,' GET_SST - ERROR while reading variable: ',a,2x,
     &                                      ' at TIME index = ',i4)
  7   write(stdout,8) sst_rec, ntsst, frcname(1:lstr), tdays,
     &                                      sst_time(itsst)*sec2day
# ifdef USE_CALENDAR
     &  -sst_origin_date_in_sec*sec2day
# endif
  8   format(/,' GET_SST - ERROR: requested time record ',I4,
     &       1x,'exeeds the last available',  /,  11x,'record ',I4,
     &       1x,'in forcing NetCDF file: ', a, /,   11x,'TDAYS = ',
     &                g12.4,2x,'last available SST_TIME = ',g12.4)
  99  may_day_flag=2
      return
      end



      subroutine set_sst_tile (Istr,Iend,Jstr,Jend)
!
! Set-up sea surface temperature data for current tile.
!
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2
      real cff, cff1,cff2, cff3,cff4, val1,val2
# define SST_DATA
# include "param.h"
# include "forces.h"
# include "scalars.h"
!
# include "compute_extended_bounds.h"
!
      it1=3-itsst
      it2=itsst
      cff=time+0.5*dt
      cff1=sst_time(it2)-cff
      cff2=cff-sst_time(it1)
!
!  Load time invariant SST and dQdSST data.
!
      if (sst_cycle.lt.0.) then
        if (FIRST_RST_TIME_STEP) then
          if (lsstgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                sst(i,j)=sstg(i,j,itsst)
                dqdt(i,j)=scldqdt*dqdtg(i,j,itsst)
              enddo
            enddo
          else
            val1=sstp(itsst)
            val2=scldqdt*dqdtp(itsst)
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                sst(i,j)=val1
                dqdt(i,j)=val2
              enddo
            enddo
          endif
        endif
!
!  Time-interpolate SST and dQdSST from grided or point data.
!  Check that for the next time step [when time=time+dt] time+dt
!  is still between sst_time(it1) and sst_time(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=scldqdt/(cff1+cff2)
        cff3=cff1*cff
        cff4=cff2*cff
        cff=1./(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
        if (lsstgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=IstrR,IendR
              sst(i,j)=cff1*sstg(i,j,it1)+cff2*sstg(i,j,it2)
              dqdt(i,j)=cff3*dqdtg(i,j,it1)+cff4*dqdtg(i,j,it2)
            enddo
          enddo
        else
          val1=cff1*sstp(it1)+cff2*sstp(it2)
          val2=cff3*dqdtp(it1)+cff4*dqdtp(it2)
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=IstrR,IendR
              sst(i,j)=val1
              dqdt(i,j)=val2
            enddo
          enddo
        endif
!
!  Unable to set-up SST and dQdSST:
!  Complain about the error and signal to quit.
!
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'sst_time', tdays, sst_time(it2)*sec2day
# ifdef USE_CALENDAR
     &  -sst_origin_date_in_sec*sec2day
# endif
  1       format(/,' SET_SST - current model time exceeds ending',
     &           1x,'value for variable: ',a,/,11x,'TDAYS = ',g12.4,
     &           2x,'TEND = ',g12.4)
          may_day_flag=2
        endif
      endif
      return
      end
#else
      subroutine get_sst_empty
      return
      end
#endif /* QCORRECTION && !ANA_SST */
