!======================================================================
! 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"
#ifndef ANA_SMFLUX
                              ! Read point or grided surface momentum
      subroutine get_smflux   ! flux (wind stress) at the appropriate
                              ! time from forcing netCDF file.
# ifdef OPENACC
      use exchange_device_host
# endif
      implicit none
# include "param.h"
# include "scalars.h"
# include "ncscrum.h"
# include "forces.h"
# include "netcdf.inc"
      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 (itsms.eq.0 .or. iic.eq.0) then
        lstr=lenstr(frcname)
c**        call opencdf (frcname,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 momentum stress
! components exist as fields or scalars.
!
        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, 'sms_time', sms_tid)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) 'sms_time', frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif

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

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

! Determine whether there is cycling to reuse the input data and
! what is cycling period "sms_cycle", find initial cycling index
! "sms_ncycle" and record index "sms_rec". Set initial value for
! time index "itsms" 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 stress to
! model units: from [Newton/m^2] to [m^2/sec^2].
!
        call set_cycle (ncidfrc,   sms_tid,    ntsms,
     &                  sms_cycle, sms_ncycle, sms_rec)
        if (may_day_flag.ne.0) return    !-->  EXIT
        itsms=2
        sms_time(1)=-1.E+20
        sms_time(2)=-1.E+20
        sms_scale=1./rho0
      endif                   !<-- itsms.eq.0 .or. iic.eq.0
!
! 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-itsms
       cff=time+0.5*dt
        if (sms_time(i).le.cff .and. cff.lt.sms_time(itsms)) return
        ierr=advance_cycle (sms_cycle, ntsms, sms_ncycle, sms_rec)
        if (ierr.ne.0) goto 7                           !--> ERROR

        ierr=nf_get_var1_FTYPE(ncidfrc, sms_tid, sms_rec, cff)
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'sms_time', sms_rec
          goto 99                                         !--> ERROR
        endif
# ifdef USE_CALENDAR
        call tool_origindate(ncidfrc,sms_tid,
     &                       sms_origin_date_in_sec)
        cff=cff+sms_origin_date_in_sec*sec2day
# endif
        sms_time(i)=cff*day2sec+sms_cycle*sms_ncycle
        if (sms_time(itsms).eq.-1.E+20) sms_time(itsms)=sms_time(i)

        if (lsusgrd.eq.1) then
          ierr=nf_fread(sustrg(START_2D_ARRAY,i), ncidfrc, susid,
# ifdef STRESS_AT_RHO_POINTS
     &                                              sms_rec, r2dvar)
# else
     &                                              sms_rec, u2dvar)
# endif
# ifdef OPENACC
          call update_device(sustrg(:,:,i))
# endif        
        else
          ierr=nf_get_var1_FTYPE(ncidfrc,susid,sms_rec,sustrp(i))
#ifdef OPENACC
          call update_device(sustrp)
#endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'sustr', sms_rec
          goto 99                                         !--> ERROR
        endif

        if (lsvsgrd.eq.1) then
          ierr=nf_fread(svstrg(START_2D_ARRAY,i), ncidfrc, svsid,
# ifdef STRESS_AT_RHO_POINTS
     &                                              sms_rec, r2dvar)
# else
     &                                              sms_rec, v2dvar)
# endif
# ifdef OPENACC
          call update_device(svstrg(:,:,i))
# endif
        else
          ierr=nf_get_var1_FTYPE(ncidfrc,svsid,sms_rec,svstrp(i))
#ifdef OPENACC
          call update_device(svstrp)
#endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'svstr', sms_rec
          goto 99                                         !--> ERROR
        endif

       itsms=i
       MPI_master_only write(stdout,'(6x,A,1x,A,1x,g12.4,1x,I4)')
     &              'GET_SMFLUX --',
     &              'Read surface momentum stresses for time =', cff
# ifdef USE_CALENDAR
     &  -sms_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                                      , mynode
# endif
      if (ntsms.gt.1) goto 1
      if (ntsms.eq.1) return
!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_SMFLUX - unable to find forcing variable: ',A,
     &                          /,14x,'in forcing netCDF file: ',A)
      goto 99
  4   write(stdout,5) frcname(1:lstr)
  5   format(/,' GET_SMFLUX - unable to open forcing netCDF file:',
     &                                                        1x,A)
  6   format(/,' GET_SMFLUX - error while reading variable: ',A,2x,
     &                                      ' at TIME index = ',i4)
  7   write(stdout,8) sms_rec, ntsms, frcname(1:lstr), tdays,
     &                                sms_time(itsms)*sec2day
# ifdef USE_CALENDAR
     &  -sms_origin_date_in_sec*sec2day
# endif
  8   format(/,' GET_SMFLUX - ERROR: requested time record ',I4,
     &       1x,'exeeds the last available', /,14x,'record ',I4,
     &       1x,'in forcing netCDF file: ',a, /,14x,'TDAYS = ',
     &       g12.4,2x,'last available SMS_TIME = ',g12.4)
 99   may_day_flag=2
      return
      end

      subroutine set_smflux_tile (Istr,Iend,Jstr,Jend)
# define SMFLUX_DATA
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2
      real cff, cff1, cff2, ustress, vstress
# include "param.h"
# ifdef SMFLUX_CURVGRID
#  include "grid.h"
# endif
# include "forces.h"
# include "scalars.h"
# ifdef SFLUX_CFB
      real stau
      parameter (stau=-0.01)
#  include "ocean3d.h"
# endif
!
# include "compute_extended_bounds.h"
!
      it1=3-itsms
      it2=itsms
      cff=time+0.5*dt
      cff1=sms_time(it2)-cff
      cff2=cff-sms_time(it1)
!
! Load time invariant surface stresses.
!
      if (sms_cycle.lt.0.) then
        if (FIRST_RST_TIME_STEP) then
          if (lsusgrd.eq.1 .and. lsvsgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=Istr,IendR
                sustr(i,j)=sms_scale*sustrg(i,j,itsms)
              enddo
            enddo
!$acc parallel loop if(compute_on_device) default(present)
            do j=Jstr,JendR
              do i=IstrR,IendR
                svstr(i,j)=sms_scale*svstrg(i,j,itsms)
              enddo
            enddo
          else
            ustress=sms_scale*sustrp(itsms)
            vstress=sms_scale*svstrp(itsms)
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=Istr,IendR
# ifdef SMFLUX_CURVGRID
                cff=0.5*(angler(i,j)+angler(i-1,j))
                sustr(i,j)=ustress*cos(cff)+vstress*sin(cff)
# else
                sustr(i,j)=ustress
# endif
              enddo
            enddo
!$acc parallel loop if(compute_on_device) default(present)
            do j=Jstr,JendR
              do i=IstrR,IendR
# ifdef SMFLUX_CURVGRID
                cff=0.5*(angler(i,j)+angler(i,j-1))
                svstr(i,j)=-ustress*sin(cff)+vstress*cos(cff)
# else
                svstr(i,j)=vstress
# endif
              enddo
            enddo
          endif
        endif
!
! Interpolate surface momentum stresses in time.
!-----------------------------------------------
!
      elseif (cff1.ge.0. .and. cff2.ge.0.) then
        if (ZEROTH_TILE .and. cff1.lt.dt) synchro_flag=.TRUE.
        cff=sms_scale/(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
!
! Process grided data. Grided stresses are assumed to be already
! rotated by the appropriate angle in curvilinear configurations.
!
        if (lsusgrd.eq.1 .and. lsvsgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=Istr,IendR
              sustr(i,j)=cff1*sustrg(i,j,it1)+cff2*sustrg(i,j,it2)
            enddo
          enddo
!$acc parallel loop if(compute_on_device) default(present)
          do j=Jstr,JendR
            do i=IstrR,IendR
              svstr(i,j)=cff1*svstrg(i,j,it1)+cff2*svstrg(i,j,it2)
            enddo
          enddo
!
! Process point data. If appropriate, rotate stress components.
!
        else
          ustress=cff1*sustrp(it1)+cff2*sustrp(it2)
          vstress=cff1*svstrp(it1)+cff2*svstrp(it2)
# ifdef SMFLUX_CURVGRID
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=Istr,IendR
              cff=0.5*(angler(i,j)+angler(i-1,j))
              sustr(i,j)=ustress*cos(cff)+vstress*sin(cff)
            enddo
          enddo
!$acc parallel loop if(compute_on_device) default(present)
          do j=Jstr,JendR
            do i=IstrR,IendR
              cff=0.5*(angler(i,j)+angler(i,j-1))
              svstr(i,j)=-ustress*sin(cff)+vstress*cos(cff)
            enddo
          enddo
# else
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=Istr,IendR
              sustr(i,j)=ustress
            enddo
          enddo
!$acc parallel loop if(compute_on_device) default(present)
          do j=Jstr,JendR
            do i=IstrR,IendR
              svstr(i,j)=vstress
            enddo
          enddo
# endif /* SMFLUX_CURVGRID */
        endif
# if defined SFLUX_CFB && !defined BULK_FLUX
!
!$acc parallel loop if(compute_on_device) default(present)
        do j=JstrR,JendR
          do i=Istr,IendR
             sustr(i,j)=sustr(i,j) + sms_scale*stau*u(i,j,N,nstp)
          enddo
        enddo
!$acc parallel loop if(compute_on_device) default(present)
        do j=Jstr,JendR
          do i=IstrR,IendR
             svstr(i,j)=svstr(i,j) + sms_scale*stau*v(i,j,N,nstp)
          enddo
        enddo
!
# endif
!
! Unable to set-up surface momentum stresses: Complain about the
! error and signal to quit (ONE THREAD ONLY).
!
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'sms_time', tdays, sms_time(it1)*sec2day
# ifdef USE_CALENDAR
     &  -sms_origin_date_in_sec*sec2day
# endif
     &                                       ,sms_time(it2)*sec2day
# ifdef USE_CALENDAR
     &  -sms_origin_date_in_sec*sec2day
# endif
  1       format(/,' SET_SMFLUX_TILE - current model time',
     &           1x,'exceeds ending value for variable: ',a,/,
     &           14x,'TDAYS = ',g12.4,2x,'SMS_TSTART = ',g12.4,2x,
     &                                        'SMS_TEND = ',g12.4)
          may_day_flag=2
        endif
      endif
      return
      end
#else
      subroutine get_smflux_empty
      return
      end
#endif /* !ANA_SMFLUX */
