!======================================================================
! 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"
#ifdef SOLVE3D

      subroutine set_scoord     ! Input:   hmin,   Tcline,
!                               !          theta_s, theta_b
! Define S-coordinate system.   !
!                               ! Output:  hc, sc_w(0:N), Cs_w(0:N)
                                !              sc_r(1:N), Cs_r(1:N)
# ifdef OPENACC
      use exchange_device_host
# endif
      implicit none
# include "param.h"
# include "scalars.h"
      integer k
      real cff,cff1,cff2,cff3, ds, sc, csf
# ifdef BSTRESS_FAST
#  include "nbq.h"
      real Hz0,r_D,umag
# endif
!
! Set S-Curves in domain [-1 < sc < 0] at vertical W- and RHO-points.
!
# ifndef NEW_S_COORD
      hc=min(hmin,Tcline)
      cff1=1./sinh(theta_s)
      cff2=0.5/tanh(0.5*theta_s)

      sc_w(0)=-1.0
      Cs_w(0)=-1.0

      cff=1./float(N)
      do k=1,N,+1
        sc_w(k)=cff*float(k-N)
        Cs_w(k)=(1.-theta_b)*cff1*sinh(theta_s*sc_w(k))
     &             +theta_b*(cff2*tanh(theta_s*(sc_w(k)+0.5))-0.5)

        sc_r(k)=cff*(float(k-N)-0.5)
        Cs_r(k)=(1.-theta_b)*cff1*sinh(theta_s*sc_r(k))
     &             +theta_b*(cff2*tanh(theta_s*(sc_r(k)+0.5))-0.5)
      enddo
#  ifdef M3FAST_SEDLAYERS
      if (N_sl.gt.0) then
        sc_r(0)=cff*(float(0-N)-0.5)
        Cs_r(0)=(1.-theta_b)*cff1*sinh(theta_s*sc_r(0))
     &             +theta_b*(cff2*tanh(theta_s*(sc_r(0)+0.5))-0.5)
     
       do k=-N_sl+1,-1,+1
        sc_w(k)=cff*float(k-N)
        Cs_w(k)=(1.-theta_b)*cff1*sinh(theta_s*sc_w(k))
     &             +theta_b*(cff2*tanh(theta_s*(sc_w(k)+0.5))-0.5)

        sc_r(k)=cff*(float(k-N)-0.5)
        Cs_r(k)=(1.-theta_b)*cff1*sinh(theta_s*sc_r(k))
     &             +theta_b*(cff2*tanh(theta_s*(sc_r(k)+0.5))-0.5)
       enddo
       k=-N_sl
       sc_w(k)=cff*float(k-N)
       Cs_w(k)=(1.-theta_b)*cff1*sinh(theta_s*sc_w(k))
     &             +theta_b*(cff2*tanh(theta_s*(sc_w(k)+0.5))-0.5)
     
      endif
#  endif
# else
      hc= Tcline
      ds=1./float(N)
      do k=1,N,+1
        sc     = ds*(float(k-N)-0.5)
        Cs_r(k)=CSF(sc, theta_s,theta_b)
        sc_r(k)=sc
      enddo

      sc_w(0) = -1.0
      sc_w(N) =  0.
      Cs_w(0) = -1.0
      Cs_w(N) =  0.
#  ifdef M3FAST_SEDLAYERS
      do k=-N_sl+1,0
        sc     = ds*(float(k-N)-0.5)
        Cs_r(k)=CSF(sc, theta_s,theta_b)
        sc_r(k)=sc      
      enddo
#  endif
      do k=1,N-1,+1
        sc     = ds*float(k-N)
        Cs_w(k)=CSF(sc, theta_s,theta_b)
        sc_w(k)=sc
      enddo
#  ifdef M3FAST_SEDLAYERS
      do k=-N_sl,-1,+1
        sc     = ds*float(k-N)
        Cs_w(k)=CSF(sc, theta_s,theta_b)
        sc_w(k)=sc      
      enddo    
#  endif     
# endif
!
! Report information about vertical S-levels.
!
# ifdef WET_DRY
      if (hmin .lt. D_wetdry) hmin=D_wetdry
# endif
      MPI_master_only write(stdout,'(/1x,A/,/1x,A,10x,A/)')
     &                       'Vertical S-coordinate System:',
     &                       'level   S-coord     Cs-curve',
     &                       'at_hmin  over_slope     at_hmax'

      do k=N,-N_sl,-1
# ifdef NEW_S_COORD
        sc=ds*(k-N)
        cff1 = hmin*(hc*sc + hmin*Cs_w(k))/(hc+hmin)
        cff2 = 0.5*hmax*(hc*sc + 0.5*hmax*Cs_w(k))/(hc+0.5*hmax)
        cff3 = hmax*(hc*sc + hmax*Cs_w(k))/(hc+hmax)
# else
        cff1=sc_w(k)*hc+(hmin-hc)*Cs_w(k)
        cff2=sc_w(k)*hc+(0.5*(hmin+hmax)-hc)*Cs_w(k)
        cff3=sc_w(k)*hc+(hmax-hc)*Cs_w(k)
# endif
        MPI_master_only write(stdout,'(I6,2F12.7,4x,3F12.3)')
     &                     k, sc_w(k),Cs_w(k), cff1,cff2,cff3
      enddo

# ifdef BSTRESS_FAST
#  ifdef BBL
      inc_faststep=1
#  else
      Hz0=(hc*(Sc_r(1)-Sc_w(0))+hmin*(Cs_r(1)-Cs_w(0)))
#   ifdef NEW_S_COORD
     &                                  *hmin/(hc+hmin)
#   endif
      umag=0.5
      if (Zobt.ne.0.) then
        cff=vonKar/log(max(Hz0,Zobt+1.e-4)/Zobt)
        r_D=umag*MIN(Cdb_max,MAX(Cdb_min,cff*cff))
      else if (rdrg2.gt.0.) then
        r_D=rdrg+rdrg2*umag
      else if (rdrg.gt.0.0) then
        r_D=rdrg
      endif
      inc_faststep=MAX(1,MIN(inc_faststep_max,
     &                       INT(ndtfast*MIN(0.5, Hz0/(r_D*dt) ))))
#  endif
      MPI_master_only write(stdout,'(/1x,A,I3)')
     &         'Friction 3D Fast Timestep Increment: inc_faststep = ',
     &                                               inc_faststep
# endif /* BSTRESS_FAST */

# ifdef OPENACC
      call update_device(sc_w)
      call update_device(sc_r)
      call update_device(Cs_w)
      call update_device(Cs_r)
# endif
      return
      end
!======================================================================
# ifdef NEW_S_COORD
                                            ! NOTE: Mathematical
      function CSF (sc, theta_s,theta_b)    ! limits of CSF,csrf for
      implicit none                         ! theta_s, theta_b --> 0
      real CSF, sc, theta_s,theta_b,csrf    ! match that under "else"
                                            ! logical branches.
      if (theta_s.gt.0.D0) then
        csrf=(1.D0-cosh(theta_s*sc))/(cosh(theta_s)-1.D0)
      else
        csrf=-sc**2
      endif
      if (theta_b.gt.0.D0) then
        CSF=(exp(theta_b*csrf)-1.D0)/(1.D0-exp(-theta_b))
      else
        CSF=csrf
      endif
      return
      end
# endif /* NEW_S_COORD */
!======================================================================
#else
      subroutine set_scoord_empty
      return
      end
#endif /* SOLVE3D */
