!======================================================================
! 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 SOLVE3D && !defined ANA_SRFLUX && defined TEMPERATURE

                              ! Read point or grided shortwave
      subroutine get_srflux   ! radiation flux at the appropriate
                              ! time from forcing NetCDF file.
# ifdef OPENACC
      use exchange_device_host
# endif
# define SRFLUX_DATA
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "ncscrum.h"
      real cff
      integer i,ierr, lstr,lvar, lenstr, nf_fread, advance_cycle
# include "netcdf.inc"
!
! Initialization: Inquire about the contents of forcing NetCDF file:
!================  variables and dimensions. Check for consistency.
!
      if (may_day_flag.ne.0) return      !-->  EXIT
      if (itsrf.eq.0 .or. iic.eq.0) then
        lstr=lenstr(frcname)
        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 radiation flux 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, 'srf_time', srf_tid)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) 'srf_time', frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif
!
! srflx
!
        lvar=lenstr(vname(1,indxShflx_rsw))
        ierr=nf_inq_varid (ncidfrc,vname(1,indxShflx_rsw)
     &                                               (1:lvar),srf_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, srf_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              lsrfgrd=1
            else
              lsrfgrd=0
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxShflx_rsw)(1:lvar),frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# ifdef DIURNAL_INPUT_SRFLX
!
! srflxbio : for PISCES daily averaged short-wave (solar) radiation )
!
        lvar=lenstr(vname(1,indxShflx_rswbio))
        ierr=nf_inq_varid (ncidfrc,vname(1,indxShflx_rswbio)(1:lvar)
     &                                              ,srfbio_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, srfbio_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              lsrfbiogrd=1
            else
              lsrfbiogrd=0
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxShflx_rswbio)(1:lvar),
     &                                               frcname(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif /* DIURNAL_INPUT_SRFLX */
!
! Determine whether there is cycling to reuse the input data
! and find cycling period "srf_cycle", set initial cycling
! index "srf_ncycle" and record index "srf_rec".
! Set initial value for time index "itsrf" 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 flux to model units:
! convert from Watts meter-2 to Celsius meter second-1.
!
        call set_cycle (ncidfrc,   srf_tid,    ntsrf,
     &                  srf_cycle, srf_ncycle, srf_rec)
        if (may_day_flag.ne.0) return    !-->  EXIT
        itsrf=2
        srf_time(1)=-1.E+20
        srf_time(2)=-1.E+20
        srf_scale=1./(rho0*Cp)
      endif                   !<-- itsrf.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-itsrf
       cff=time+0.5*dt
        if (srf_time(i).le.cff .and. cff.lt.srf_time(itsrf)) return
        ierr=advance_cycle (srf_cycle, ntsrf, srf_ncycle, srf_rec)
        if (ierr.ne.0) goto 7                           !--> ERROR

        ierr=nf_get_var1_FTYPE(ncidfrc, srf_tid, srf_rec, cff)
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'srf_time', srf_rec
          goto 99                                         !--> ERROR
        endif
# ifdef USE_CALENDAR
        call tool_origindate(ncidfrc,srf_tid,
     &                       srf_origin_date_in_sec)
        cff=cff+srf_origin_date_in_sec*sec2day
# endif
        srf_time(i)=cff*day2sec+srf_cycle*srf_ncycle
        if (srf_time(itsrf).eq.-1.E+20) srf_time(itsrf)=srf_time(i)
!
! srflx
!
        if (lsrfgrd.eq.1) then
          ierr=nf_fread (srflxg(START_2D_ARRAY,i), ncidfrc, srf_id,
     &                                             srf_rec, r2dvar)
# ifdef OPENACC
          call update_device(srflxg(:,:,i))
# endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc,srf_id,srf_rec,srflxp(i))
        endif

        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'srflux', srf_rec
          goto 99                                         !--> ERROR
        endif

# ifdef DIURNAL_INPUT_SRFLX
!
! srflxbio (=radswbio with bulk) for PISCES daily averaged short-wave (solar) radiation)
!
        if (lsrfbiogrd.eq.1) then
          ierr=nf_fread (srflxbiog(START_2D_ARRAY,i), ncidfrc,
     &                        srfbio_id, srf_rec,  r2dvar)
#  ifdef OPENACC
          call update_device(srflxbiog(:,:,i))
#  endif
        else
          ierr=nf_get_var1_FTYPE (ncidfrc,srfbio_id,
     &                            srf_rec ,srflxbiop(i))
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,6) 'srflux_bio', srf_rec
          goto 99                                         !--> ERROR
        endif
# endif /* DIURNAL_INPUT_SRFLX */

       itsrf=i
       MPI_master_only write(stdout,'(6x,A,1x,A,1x,g12.4,1x,I4)')
     &           'GET_SRFLUX --',
     &           'Read solar shortwave radiation for time =', cff
# ifdef USE_CALENDAR
     &  -srf_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                                      , mynode
# endif
      if (ntsrf.gt.1) goto 1
      if (ntsrf.eq.1) return
!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_SRFLUX - unable to find forcing variable: ',a,
     &                          /,15x,'in forcing NetCDF file: ',a)
  4   write(stdout,5) frcname(1:lstr)
  5   format(/,' GET_SRFLUX - unable to open forcing NetCDF file:',
     &                                                        1x,a)
      goto 99

  6   format(/,' GET_SRFLUX - error while reading variable: ',a,2x,
     &                                      ' at TIME index = ',i4)
  7   write(stdout,8) srf_rec, ntsrf, frcname(1:lstr), tdays,
     &                                srf_time(itsrf)*sec2day
# ifdef USE_CALENDAR
     &  -srf_origin_date_in_sec*sec2day
# endif
  8   format(/,' GET_SRFLUX - 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 SRF_TIME = ',g12.4)

  99  may_day_flag=2
      return
      end

      subroutine set_srflux_tile (Istr,Iend,Jstr,Jend)
!
!--------------------------------------------------------------------
!  Set-up shortwave radiation flux for current tile.
!--------------------------------------------------------------------
!
# define SRFLUX_DATA
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "grid.h"
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2
      real cff, cff1, cff2, srf
# ifdef DIURNAL_INPUT_SRFLX
      real srfbio
# endif
!
# include "compute_extended_bounds.h"
!
! Set coefficients for interpolation. Check that for the next time
! step [when time=time+dt] both weights will still be positive, and
! if not, set synchro_flag to signal that new data should be read
! from an appropriate netCDF input file (master thread only).
! After that either load time-invariant data, or interpolate in time
! or complain about error and signal to quit, if interpolation is
! needed, but not possible.
!
!
      it1=3-itsrf
      it2=itsrf
      cff=time+0.5*dt
      cff1=srf_time(it2)-cff
      cff2=cff-srf_time(it1)
      if (ZEROTH_TILE.and. cff1.lt.dt) synchro_flag=.TRUE.
!
! Load time-invariant shortwave radiation flux.
! Time interpolation is not performed in this case.
!
      if (srf_cycle.lt.0.) then
        if (FIRST_RST_TIME_STEP) then
          if (lsrfgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                srflx(i,j)=srf_scale*srflxg(i,j,itsrf)
# ifdef DIURNAL_INPUT_SRFLX
                srflxbio(i,j)=srf_scale*srflxbiog(i,j,itsrf)
# endif
              enddo
            enddo
          else
            srf=srf_scale*srflxp(itsrf)
# ifdef DIURNAL_INPUT_SRFLX
            srflxbio=srf_scale*srflxbiop(itsrf)
# endif
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
                srflx(i,j)=srf
# ifdef DIURNAL_INPUT_SRFLX
                srflxbio(i,j)=srfbio
# endif
              enddo
            enddo
          endif
        endif
!
! Time-interpolate shortwave radiation flux from grided or point
! data. Check that for the next time step [when time=time+dt]
! time+dt is still between srf_time(it1) and srf_time(it2);
! and if not, set synchro_flag to 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
        cff=srf_scale/(cff1+cff2)
        cff1=cff1*cff
        cff2=cff2*cff
        if (lsrfgrd.eq.1) then
!$acc parallel loop if(compute_on_device) default(present)
          do j=JstrR,JendR
            do i=IstrR,IendR
              srflx(i,j)=cff1*srflxg(i,j,it1)+cff2*srflxg(i,j,it2)
# ifdef DIURNAL_INPUT_SRFLX
              srflxbio(i,j)=cff1*srflxbiog(i,j,it1)
     &                                    +cff2*srflxbiog(i,j,it2)
# endif
            enddo
          enddo
        else
          srf=cff1*srflxp(it1)+cff2*srflxp(it2)
# ifdef DIURNAL_INPUT_SRFLX
          srfbio=cff1*srflxbiop(it1)
     &                        +cff2*srflxbiop(it2)
# endif
!$acc parallel loop if(compute_on_device) default(present)
            do j=JstrR,JendR
              do i=IstrR,IendR
              srflx(i,j)=srf
# ifdef DIURNAL_INPUT_SRFLX
              srflxbio(i,j)=srfbio
# endif
            enddo
          enddo
        endif
!
! Unable to set-up shortwave radiation flux:
! Complain about the error and signal to quit (ONE THREAD ONLY).
!
      else
        if (ZEROTH_TILE) then
          write(stdout,1) 'srf_time', tdays, srf_time(it2)*sec2day
# ifdef USE_CALENDAR
     &  -srf_origin_date_in_sec*sec2day
# endif
  1       format(/,' SET_SRFLUX_TILE - current model time exceeds',
     &             ' ending value for variable: ',a,/,14x,'TDAYS = ',
     &                                      g12.4,2x,'TEND = ',g12.4)
          may_day_flag=2
        endif
      endif

      return
      end
#else
      subroutine get_srflux_empty
      return
      end
#endif /* SOLVE3D && !ANA_SRFLUX */
!--------------------------------------------------------------------
#if defined SOLVE3D && defined ANA_DIURNAL_SW && defined TEMPERATURE
      subroutine ana_diurnal_sw_tile (Istr,Iend,Jstr,Jend)
!
!--------------------------------------------------------------------
!
!  DIURNAL CYCLE - USED IN BOTH PHYSICAL AND ECOSYSTEM MODELS
!  Patrick Marchesiello - 99
!
!  Modulate average daily insolation to get diurnal cycle
!  by:  pi*( cos(h)*cos(d)*cos(phi)+   sin(d)sin(phi))/
!          (sin(h0)*cos(d)*cos(phi)+h0*sin(d)sin(phi))
!  h, d, phi are hour, declination, latitude angles
!  h0 is hour angle at sunset and sunrise
!
!--------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "forces.h"
# include "scalars.h"
# include "grid.h"
      integer Istr,Iend,Jstr,Jend, i,j
      real cff,cff1,cff2,
     &     hour,cos_h,phi,h0,dec,cos_dec,sin_dec,tan_dec
!
# include "compute_extended_bounds.h"
!
! Initialize fixed parameters
!
      if (FIRST_RST_TIME_STEP) then
!$acc parallel loop if(compute_on_device) default(present)
        do j=JstrR,JendR
          do i=IstrR,IendR
            phi=latr(i,j)*deg2rad
            cos_phi(i,j)=cos(phi)
            sin_phi(i,j)=sin(phi)
            tan_phi(i,j)=tan(phi)
          enddo
        enddo
      endif
!
! Compute local hour angle (radians)
!
      hour=2.*pi*(tdays+.5-int(tdays+.5))
!
! Compute declination angle
!
      dec=-0.406*cos(deg2rad*(tdays-
     &               int(tdays*day2year)*year2day))
      cos_dec=cos(dec)
      sin_dec=sin(dec)
      tan_dec=sin_dec/cos_dec
!
! Compute daily insolation coefficient
! Modulate srflx and update stflx
!  --> use Greenwich hour angle:
!      cos_h=1 at solar noon in Greenwich
!
!$acc parallel loop if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          cos_h=cos(hour-lonr(i,j)*deg2rad)
          h0=acos(-tan_phi(i,j)*tan_dec)
          cff1=cos_dec*cos_phi(i,j)
          cff2=sin_dec*sin_phi(i,j)
          cff=pi*(cos_h*cff1+cff2)/(sin(h0)*cff1+h0*cff2)

          stflx(i,j,itemp)=stflx(i,j,itemp)-srflx(i,j)
          srflx(i,j)=max(0.,cff*srflx(i,j))
          stflx(i,j,itemp)=stflx(i,j,itemp)+srflx(i,j)

# ifdef BULK_FLUX
!
!  Correct fluxes for diagnostics
!
          shflx_rsw(i,j)=srflx(i,j)
# endif

        enddo
      enddo

      return
      end
#else
      subroutine diurnal_srflux_empty
      return
      end
#endif /* SOLVE3D && ANA_DIURNAL_SW && TEMPERATURE */
