!======================================================================
! 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 SALINITY && defined SFLX_CORR && !defined ANA_SSS

      subroutine get_sss
!
!  Read in point or grided sea surface salinity used to correct
!  the salinity fluxes at the appropriate
!  time from forcing NetCDF file.
!
!  These forcing fields are used when salt flux correction is activated:
!
!    SSSFLX_model ~ SSS*(E-P) + CST * (SSS_model - SSS)
!
!  we use DQDSST for CST....
!
# ifdef OPENACC
      use exchange_device_host
# endif
# define SSS_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 (itsss.eq.0 .or. iic.eq.0) then
        lstr=lenstr(frcname)
!
! If not opened yet, open forcing NetCDF file for reading.
! Find and save IDs for relevant variables, determine whether
! SSS 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, 'sss_time', sss_tid)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) 'sss_time', frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif

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

# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
        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
# endif
!
! Determine whether there is cycling to reuse the input data and
! find cycling period "sss_cycle", set initial cycling index
! "sss_ncycle" and record index "sss_rec".
! Set initial value for time index "itsss" 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,   sss_tid,    ntsss,
     &                  sss_cycle, sss_ncycle, sss_rec)
        if (may_day_flag.ne.0) return    !-->  EXIT
        itsss=2
        sss_time(1)=-1.E+20
        sss_time(2)=-1.E+20
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
        scldqdt=1./(rho0*Cp)
# 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-itsss
       cff=time+0.5*dt
        if (sss_time(i).le.cff .and. cff.lt.sss_time(itsss)) return
        ierr=advance_cycle (sss_cycle, ntsss, sss_ncycle, sss_rec)
        if (ierr .ne. 0) goto 7                           !--> ERROR

        ierr=nf_get_var1_FTYPE (ncidfrc, sss_tid, sss_rec, cff)
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'sss_time', sss_rec
          goto 99                                         !--> ERROR
        endif
# ifdef USE_CALENDAR
            call tool_origindate(ncidfrc,sss_tid,
     &                           sss_origin_date_in_sec)
              cff=cff+sss_origin_date_in_sec*sec2day
# endif
        sss_time(i)=cff*day2sec+sss_cycle*sss_ncycle
        if (sss_time(itsss).eq.-1.E+20) sss_time(itsss)=sss_time(i)

        if (lsssgrd.eq.1) then
          ierr=nf_fread (sssg(START_2D_ARRAY,i), ncidfrc, sss_id,
     &                                           sss_rec, r2dvar)
# ifdef OPENACC
          call update_device(sssg(:,:,i))
# endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc, sss_id, sss_rec, sssp(i))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'SSS', sss_rec
          goto 99                                         !--> ERROR
        endif

# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
        if (lsssgrd.eq.1) then
          ierr=nf_fread (dqdtg(START_2D_ARRAY,i), ncidfrc, dqdt_id,
     &                                            sss_rec,  r2dvar)
#  ifdef OPENACC
          call update_device(dqdtg(:,:,i))
#  endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc,dqdt_id,sss_rec,dqdtp(i))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'dQdSST', sss_rec
          goto 99                                         !--> ERROR
        endif
# endif
       itsss=i
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
      MPI_master_only write(stdout,'(6x,A,1x,A,1x,g12.4,1x,I4)')
     &              'GET_SSS    --',
     &              'Read SSS and dQdSST fields     for time =', cff
# else
      MPI_master_only write(stdout,'(6x,A,16x,A,1x,g12.4,1x,I4)')
     &                  'GET_SSS    -- Read SSS fields',
     &                                     'for time =', cff
# endif
# ifdef USE_CALENDAR
     &  -sss_origin_date_in_sec*sec2day
# endif

# ifdef MPI
     &                                                 , mynode
# endif
      if (ntsss.gt.1) goto 1
      if (ntsss.eq.1) return
!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_SSS - ERROR: unable to find forcing variable',
     &                   ': ',a,/,11x,'in forcing NetCDF file: ',a)
  4   write(stdout,5) frcname(1:lstr)
  5   format(/,' GET_SSS - ERROR: unable to open forcing NetCDF ',
     &                                                  'file: ',a)
      goto 99
  6   format(/,' GET_SSS - ERROR while reading variable: ',a,2x,
     &                                      ' at TIME index = ',i4)
  7   write(stdout,8) sss_rec, ntsss, frcname(1:lstr), tdays,
     &                                      sss_time(itsss)*sec2day
  8   format(/,' GET_SSS - 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 SSS_TIME = ',g12.4)
  99  may_day_flag=2
      return
      end



      subroutine set_sss_tile (Istr,Iend,Jstr,Jend)
!
! Set-up sea surface salinity data for current tile.
!
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2
      real cff, cff1,cff2, cff3,cff4, val1, val2
# define SSS_DATA
# define SST_DATA /* PIERRICK */
# include "param.h"
# include "forces.h"
# include "scalars.h"
!
# include "compute_extended_bounds.h"
!
      it1=3-itsss
      it2=itsss
      cff=time+0.5*dt
      cff1=sss_time(it2)-cff
      cff2=cff-sss_time(it1)
!
!  Load time invariant SSS data.
!
      if (sss_cycle.lt.0.) then
        if (FIRST_RST_TIME_STEP) then
          if (lsssgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                sss(i,j)=sssg(i,j,itsss)
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
                dqdt(i,j)=scldqdt*dqdtg(i,j,itsss)
# endif
              enddo
            enddo
          else
            val1=sssp(itsss)
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
            val2=scldqdt*dqdtp(itsss)
# endif
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                sss(i,j)=val1
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
                dqdt(i,j)=val2
# endif
              enddo
            enddo
          endif
        endif
!
!  Time-interpolate SSS from grided or point data.
!  Check that for the next time step [when time=time+dt] time+dt
!  is still between sss_time(it1) and sss_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.
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
        cff=scldqdt/(cff1+cff2)
        cff3=cff1*cff
        cff4=cff2*cff
# endif
        cff=1./(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
        if (lsssgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=IstrR,IendR
              sss(i,j)=cff1*sssg(i,j,it1)+cff2*sssg(i,j,it2)
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
              dqdt(i,j)=cff3*dqdtg(i,j,it1)+cff4*dqdtg(i,j,it2)
# endif
            enddo
          enddo
        else
          val1=cff1*sssp(it1)+cff2*sssp(it2)
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
          val2=cff3*dqdtp(it1)+cff4*dqdtp(it2)
# endif
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=IstrR,IendR
              sss(i,j)=val1
# if defined SFLX_CORR && !defined SFLX_CORR_COEF && !defined QCORRECTION
              dqdt(i,j)=val2
# endif
            enddo
          enddo
        endif
!
!  Unable to set-up SSS:
!  Complain about the error and signal to quit.
!
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'sss_time', tdays, sss_time(it2)*sec2day
  1       format(/,' SET_SSS - 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_sss_empty
      return
      end
#endif /* SALINITY && SFLX_CORR && !defined ANA_SSS */
