!======================================================================
! 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
!======================================================================
!
!
!  This is the "online_set_bulk.F" script
!------------------------------------------------------------------------------
!  This file contains the subfunctions enabling the ONLINE extraction of the
!  forcing from an atmospheric dataset. A spatial and time interpolation are 
!  applied to the extracted data in order to adapt these to the considered 
!  simulation domain as well as the associated MPI/OPENMP discretisation (if 
!  defined MPI/OPENMP).
!------------------------------------------------------------------------------
#include "cppdefs.h"

#if defined BULK_FLUX && defined ONLINE

!*******************************************************************************
      subroutine set_bulk_tile_online (Istr,Iend,Jstr,Jend)
!------------------------------------------------------------------------------
!  This subfunction enables the time interpolation of the roms bulk forcings
!  (tair, rhum, prate, radlw, radsw, uwnd, vwnd) on the simulation domain
!  for the ongoing tile.
!------------------------------------------------------------------------------
!  The steps in this interpolation is:
!  - Determination of forcings (tair, rhum, prate, radlw, radsw, uwnd,
!  vwnd) at the first time step of the simulation from the known (tairg, rhumg,
!  prateg, radlwg, radswg, uwndg, vwndg) fields.
!  - Time interpolation of the forcings (tair, rhum, prate, radlw, radsw,
!  uwnd, vwnd) at cff from two known spatially interpolated fields (tairg,rhumg,
!  prateg, radlwg, radswg, uwndg, vwndg) at cff1 and cff2.
!
!  In this subroutine the aformentioned steps are completed for every bulk
!  forcing iterating on the various blkvar_id.
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "grid.h"
# include "online.h"
# ifdef CFB_WIND_TRA
#  include "ocean3d.h"
#  include "params_bulk.h"
# endif
# ifdef ABL1D
#  include "abl1d.h"
# endif
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2, blkvar_id
# ifdef ABL1D
      integer k
# endif
      real cff,cff1,cff2, cff3,cff4
# ifdef SALINITY
      real cff5,cff6
# endif
# ifdef CFB_WIND_TRA
      real cffu1,cffu2,cffv1,cffv2
      real uwnd_r, vwnd_r
# endif
!
# include "compute_extended_bounds.h"
!

      !===== == === === ======= =========
      ! Loop on all the forcing datasets:
      !===== == === === ======= =========
      blkvar_id=0
  10  blkvar_id = blkvar_id+1
!
! Do not process #7 :  upward longwave is obsolete
!
      if (blkvar_id.eq.7) goto 10
# if defined ERA_ECMWF || defined AROME
!
! In the case of ERA the net short wave is used
! Net Short wave - 5 downward short wave not used
!
      if (blkvar_id.eq.4) goto 10
# endif

# ifndef READ_PATM
      if (blkvar_id.eq.10) goto 10
# endif
      if(blkvar_id.gt.nblkvrs) return

      it1=3-itbulkO(blkvar_id)
      it2=itbulkO(blkvar_id)
      ! --------------------------
      ! Times of the interpolation
      ! --------------------------
      cff=time+0.5*dt
      cff1=bulk_timeO(it2,blkvar_id)-cff
      cff2=cff-bulk_timeO(it1,blkvar_id)


      !===== ============= ======= === ==== ======
      ! Time interpolation between two time steps:
      !===== ============= ======= === ==== ======
      ! The roms bulk forcings at cff are interpolated between two time steps
      ! cff1 and cff2
      if (cff1.ge.0. .and. cff2.ge.0.) then
        if (ZEROTH_TILE .and. cff1.lt.dt) synchro_flag=.TRUE.
                                    !note cff order maters
        cff=srf_scale/(cff1+cff2)
        cff3=cff1*cff
        cff4=cff2*cff
# ifdef SALINITY
        cff=stf_scale(isalt)/(cff1+cff2)
        cff5=cff1*cff
        cff6=cff2*cff
# endif
        cff=1./(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
        if (blkvar_id.eq.1) then
          do j=JstrR,JendR
            do i=IstrR,IendR
# ifdef ABL1D
              do k=1,N_abl
                pt_dta(i,j,k)=cff1*pt_dtag(i,j,k,it1)+
     &                        cff2*pt_dtag(i,j,k,it2)
              enddo
# else
              tair(i,j)=cff1*tairg(i,j,it1)+cff2*tairg(i,j,it2)
# endif
            enddo
          enddo
        elseif (blkvar_id.eq.2) then
          do j=JstrR,JendR
            do i=IstrR,IendR
# ifdef ABL1D
              do k=1,N_abl
                pq_dta(i,j,k)=cff1*pq_dtag(i,j,k,it1)+
     &                        cff2*pq_dtag(i,j,k,it2)
              enddo
# else
              rhum(i,j)=cff1*rhumg(i,j,it1)+cff2*rhumg(i,j,it2)
# endif

            enddo
          enddo
# ifdef SALINITY
        elseif (blkvar_id.eq.3) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              prate(i,j)=cff5*prateg(i,j,it1)+cff6*prateg(i,j,it2)
            enddo
          enddo
# endif
        elseif (blkvar_id.eq.5) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              radsw(i,j)=cff3*radswg(i,j,it1)+cff4*radswg(i,j,it2)
              srflx(i,j)=radsw(i,j)
            enddo
          enddo
        elseif (blkvar_id.eq.6) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              radlw(i,j)=cff3*radlwg(i,j,it1)+cff4*radlwg(i,j,it2)
            enddo
          enddo
        elseif (blkvar_id.eq.9) then
          do j=JstrR,JendR
            do i=IstrR,IendR
# ifdef ABL1D
              do k=1,N_abl
                pu_dta(i,j,k)=cff1*pu_dtag(i,j,k,it1)+
     &                        cff2*pu_dtag(i,j,k,it2)
                pv_dta(i,j,k)=cff1*pv_dtag(i,j,k,it1)+
     &                        cff2*pv_dtag(i,j,k,it2)
              enddo
            enddo
          enddo
# else
              uwnd(i,j)=cff1*uwndg(i,j,it1)+cff2*uwndg(i,j,it2)
              vwnd(i,j)=cff1*vwndg(i,j,it1)+cff2*vwndg(i,j,it2)
              wspd(i,j)=cff1*wspdg(i,j,it1)+cff2*wspdg(i,j,it2)
            enddo
          enddo
#  ifdef CFB_WIND_TRA
          cff      = 1.-swparam    ! current-wind coupling parameter: Ua => Ua-(1-sw)Uo

          do j=JstrR,min(JendR,Mm)
            do i=IstrR,min(IendR,Lm)
              cffu1=cff1*uwndg(i  ,j,it1)+cff2*uwndg(i  ,j,it2)
              cffu2=cff1*uwndg(i+1,j,it1)+cff2*uwndg(i+1,j,it2)
              uwnd_r   = 0.5*( cffu1 + cffu2
     &                 - cff*( u(i+1,j,N,nrhs)+u(i,j,N,nrhs))
     &                   )

              cffv1=cff1*vwndg(i,j  ,it1)+cff2*vwndg(i,j  ,it2)
              cffv2=cff1*vwndg(i,j+1,it1)+cff2*vwndg(i,j+1,it2)
              vwnd_r   = 0.5*( cffv1 + cffv2
     &                 - cff*( v(i,j+1,N,nrhs)+v(i,j,N,nrhs))
     &                   )

              wspd_cfb(i,j) = SQRT( uwnd_r*uwnd_r+vwnd_r*vwnd_r )
            enddo
          enddo
#   ifndef EW_PERIODIC
            if (WESTERN_EDGE) then
              do j=Jstr,Jend
                wspd_cfb(Istr-1,j)=wspd_cfb(Istr,j)
              enddo
            endif
            if (EASTERN_EDGE) then
              do j=Jstr,Jend
                wspd_cfb(Iend+1,j)=wspd_cfb(Iend,j)
              enddo
            endif
#   endif
#   ifndef NS_PERIODIC
            if (SOUTHERN_EDGE) then
              do i=Istr,Iend
                wspd_cfb(i,Jstr-1)=wspd_cfb(i,Jstr)
              enddo
            endif
            if (NORTHERN_EDGE) then
              do i=Istr,Iend
                wspd_cfb(i,Jend+1)=wspd_cfb(i,Jend)
              enddo
            endif
#    ifndef EW_PERIODIC
            if (WESTERN_EDGE .and. SOUTHERN_EDGE) then
              wspd_cfb(Istr-1,Jstr-1)=wspd_cfb(Istr,Jstr)
            endif
            if (WESTERN_EDGE .and. NORTHERN_EDGE) then
              wspd_cfb(Istr-1,Jend+1)=wspd_cfb(Istr,Jend)
            endif
            if (EASTERN_EDGE .and. SOUTHERN_EDGE) then
              wspd_cfb(Iend+1,Jstr-1)=wspd_cfb(Iend,Jstr)
            endif
            if (EASTERN_EDGE .and. NORTHERN_EDGE) then
              wspd_cfb(Iend+1,Jend+1)=wspd_cfb(Iend,Jend)
            endif
#    endif
#   endif
#   if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
            call exchange_r2d_tile(Istr,Iend,Jstr,Jend,
     &                             wspd_cfb(START_2D_ARRAY))
#   endif
#  endif
# endif
# ifdef READ_PATM
        elseif (blkvar_id.eq.10) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              patm2d(i,j)=cff1*patmg(i,j,it1)+cff2*patmg(i,j,it2)
            enddo
          enddo
# endif
# ifdef ABL1D
        elseif (blkvar_id.eq.11) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              do k=1,N_abl
                pgu_dta(i,j,k)=cff1*pgu_dtag(i,j,k,it1)+
     &                         cff2*pgu_dtag(i,j,k,it2)
              enddo
            enddo
          enddo
        elseif (blkvar_id.eq.12) then
          do j=JstrR,JendR
            do i=IstrR,IendR
              do k=1,N_abl
                pgv_dta(i,j,k)=cff1*pgv_dtag(i,j,k,it1)+
     &                         cff2*pgv_dtag(i,j,k,it2)
              enddo
            enddo
          enddo
# endif
        endif

      !====== == ==== ==============
      ! Error in time interpolation:
      !====== == ==== ==============
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'bulk_timeO',tdays,bulk_timeO(it2,blkvar_id)*sec2day
  1       format(/,' SET_BULK - current model time exceeds ending',
     &           1x,'value for variable: ',a,/,11x,'TDAYS = ',g12.4,
     &           2x,'TEND = ',g12.4)
          may_day_flag=2
        endif
      endif

      goto 10

      return
      end

#endif /* BULK_FLUX && ONLINE */
