!======================================================================
! 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 LMD_SKPP || defined LMD_BKPP || defined GLS_MIXING

      subroutine lmd_swfrac_tile (Istr,Iend,Jstr,Jend, Zscale,Z,swdk)
!
!--------------------------------------------------------------------
!  This subroutine computes the  fraction of solar shortwave flux
!  penetrating to specified depth (times Zscale) due to exponential
!  decay in Jerlov water type.
!
!  Input:
!          Zscale   scale factor to apply to depth array
!          Z        vertical height [meters, negative] for
!                        desired solar short-wave fraction.
!  Output:
!          swdk     shortwave (radiation) fractional decay.
!
!  Reference:
!
!  Paulson, C.A., and J.J. Simpson, 1977: Irradiance meassurements
!  in the upper ocean, J. Phys. Oceanogr., 7, 952-956.
!
!  This routine was adapted from Bill Large 1995 code.
!--------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "scalars.h"
# include "lmd_kpp.h"
      integer Istr,Iend,Jstr,Jend, i,j,indx
      real Zscale, Z(PRIVATE_2D_SCRATCH_ARRAY),
     &             swdk(PRIVATE_2D_SCRATCH_ARRAY)
      integer lmd_Jwt
      !real lmd_mu1(5), lmd_mu2(5), lmd_r1(5), cff1, cff2, cff3 ,cff4
      real  cff1, cff2, cff3 ,cff4
      real,parameter,dimension(5)::
     &  lmd_mu1=(/0.35,0.6,1.0,1.5,1.4/)
      real,parameter,dimension(5)::
     &  lmd_mu2=(/23.,20.,17.,14.,7.9/)
      real,parameter,dimension(5)::
     &  lmd_r1=(/0.58,0.62,0.67,0.77,0.78/)

                        !  Set Jerlov water type to assign
      lmd_Jwt=1         !  everywhere; an integer from 1 to 5.

!$acc kernels if(compute_on_device) default(present) async(1)
!
!  On first call: initialize (for now) initialize index for Jerlov
!  water type to 1 (type I). This variable is used in the computation
!  of the fraction solar shortwave flux.
!
      if (FIRST_RST_TIME_STEP) then
        do j=Jstr,Jend
          do i=Istr,Iend
            Jwtype(i,j)=lmd_Jwt
          enddo
        enddo
      endif
!
!  Use Paulson and Simpson (1977) two wavelength bands solar
!  absorption model.
!
      do j=Jstr,Jend
        do i=Istr,Iend
          indx=Jwtype(i,j)
          cff1=Z(i,j)*Zscale/lmd_mu1(indx)
          cff2=Z(i,j)*Zscale/lmd_mu2(indx)
          if (cff1.ge.-20.) then            ! Quick fix to avoid computing
            cff3=lmd_r1(indx) *exp(cff1)    ! exp(-15000)
          else
            cff3=0.
          endif
          if (cff2.ge.-20.) then
            cff4=(1.-lmd_r1(indx)) *exp(cff2)
          else
            cff4=0.
          endif
          swdk(i,j)=cff3+cff4
        enddo
      enddo
!$acc end kernels
      return
      end

      subroutine lmd_swfrac_ER_tile (Istr,Iend,Jstr,Jend, Zscale,Z,swdk)
!
!--------------------------------------------------------------------
!  This subroutine computes the  fraction of solar shortwave flux
!  penetrating to specified depth (times Zscale) due to exponential
!  decay in Jerlov water type.
!
!  Input:
!          Zscale   scale factor to apply to depth array
!          Z        vertical height [meters, negative] for
!                        desired solar short-wave fraction.
!  Output:
!          swdk     shortwave (radiation) fractional decay.
!
!  Reference:
!
!  Paulson, C.A., and J.J. Simpson, 1977: Irradiance meassurements
!  in the upper ocean, J. Phys. Oceanogr., 7, 952-956.
!
!  This routine was adapted from Bill Large 1995 code.
!--------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "scalars.h"
# include "lmd_kpp.h"
      integer Istr,Iend,Jstr,Jend, i,j,indx
      integer imin,imax,jmin,jmax
      real Zscale, Z(PRIVATE_2D_SCRATCH_ARRAY),
     &             swdk(PRIVATE_2D_SCRATCH_ARRAY)
      integer lmd_Jwt
      !real lmd_mu1(5), lmd_mu2(5), lmd_r1(5), cff1, cff2, cff3 ,cff4
      real  cff1, cff2, cff3 ,cff4
!    real, parameter,dimension(5) :: lmd_mu1, lmd_mu2, lmd_r1
      real,parameter,dimension(5)::
     &  lmd_mu1=(/0.35,0.6,1.0,1.5,1.4/)
      real,parameter,dimension(5)::
     &  lmd_mu2=(/23.,20.,17.,14.,7.9/)
      real,parameter,dimension(5)::
     &  lmd_r1=(/0.58,0.62,0.67,0.77,0.78/)

                        !  Set Jerlov water type to assign
      lmd_Jwt=1         !  everywhere; an integer from 1 to 5.

# ifdef EW_PERIODIC
#  define I_EXT_RANGE Istr-1,Iend+1
# else
      if (WESTERN_EDGE) then
        imin=Istr
      else
        imin=Istr-1
      endif
      if (EASTERN_EDGE) then
        imax=Iend
      else
        imax=Iend+1
      endif
#  define I_EXT_RANGE imin,imax
# endif
# ifdef NS_PERIODIC
#  define J_EXT_RANGE Jstr-1,Jend+1
# else
      if (SOUTHERN_EDGE) then
        jmin=Jstr
      else
        jmin=Jstr-1
      endif
      if (NORTHERN_EDGE) then
        jmax=Jend
      else
        jmax=Jend+1
      endif
#  define J_EXT_RANGE jmin,jmax
# endif

!$acc kernels if(compute_on_device) default(present) async(1)

!
!  On first call: initialize (for now) initialize index for Jerlov
!  water type to 1 (type I). This variable is used in the computation
!  of the fraction solar shortwave flux.
!
      if (FIRST_RST_TIME_STEP) then
        do j=J_EXT_RANGE
          do i=I_EXT_RANGE
            Jwtype(i,j)=lmd_Jwt
          enddo
        enddo
      endif
!
!  Use Paulson and Simpson (1977) two wavelength bands solar
!  absorption model.
!
      do j=J_EXT_RANGE
        do i=I_EXT_RANGE
          indx=Jwtype(i,j)
          cff1=Z(i,j)*Zscale/lmd_mu1(indx)
          cff2=Z(i,j)*Zscale/lmd_mu2(indx)
          if (cff1.ge.-20.) then            ! Quick fix to avoid computing
            cff3=lmd_r1(indx) *exp(cff1)    ! exp(-15000)
          else
            cff3=0.
          endif
          if (cff2.ge.-20.) then
            cff4=(1.-lmd_r1(indx)) *exp(cff2)
          else
            cff4=0.
          endif
          swdk(i,j)=cff3+cff4
        enddo
      enddo
!$acc end kernels
      return
      end
#else
      subroutine lmd_swfrac
      return
      end
#endif /* LMD_SKPP || LMD_BKPP || GLS_MIXING */


