#include "cppdefs.h"
#ifdef SOLVE3D
      SUBROUTINE set_scoord (ng)
!
!svn $Id: set_scoord.F 332 2009-03-13 03:38:51Z arango $
!=======================================================================
!  Copyright (c) 2002-2009 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                           Hernan G. Arango   !
!========================================== Alexander F. Shchepetkin ===
!                                                                      !
!  This routine sets and initializes relevant variables associated     !
!  with the vertical terrain-following coordinates transformation.     !
!                                                                      !
!  Definitions:                                                        !
!                                                                      !
!    N(ng) : Number of vertical levels for each nested grid.           !
!                                                                      !
!     zeta : time-varying free-surface, zeta(x,y,t), (m)               !
!                                                                      !
!        h : bathymetry, h(x,y), (m, positive, maybe time-varying)     !
!                                                                      !
!       hc : critical (thermocline, pycnocline) depth (m, positive)    !
!                                                                      !
!        z : vertical depths, z(x,y,s,t), meters, negative             !
!              z_w(x,y,0:N(ng))      at   W-points  (top/bottom cell)  !
!              z_r(z,y,1:N(ng))      at RHO-points  (cell center)      !
!                                                                      !
!              z_w(x,y,0    ) = -h(x,y)                                !
!              z_w(x,y,N(ng)) = zeta(x,y,t)                            !
!                                                                      !
!        s : nondimensional stretched vertical coordinate,             !
!             -1 <= s <= 0                                             !
!                                                                      !
!              s = 0   at the free-surface, z(x,y, 0,t) = zeta(x,y,t)  !
!              s = -1  at the bottom,       z(x,y,-1,t) = - h(x,y,t)   !
!                                                                      !
!              sc_w(k) = (k-N(ng))/N(ng)       k=0:N,    W-points      !
!              sc_r(k) = (k-N(ng)-0.5)/N(ng)   k=1:N,  RHO-points      !
!                                                                      !
!        C : nondimensional vertical stretching function, C(s),        !
!              -1 <= C(s) <= 0                                         !
!                                                                      !
!              C(s) = 0    for s = 0,  at the free-surface             !
!              C(s) = -1   for s = -1, at the bottom                   !
!                                                                      !
!              Cs_w(k) = F(s,theta_s,theta_b)  k=0:N,    W-points      !
!              Cs_r(k) = C(s,theta_s,theta_b)  k=1:N,  RHO-points      !
!                                                                      !
!       Zo : vertical transformation functional, Zo(x,y,s):            !
!                                                                      !
!              Zo(x,y,s) = H(x,y)C(s)      separable functions         !
!                                                                      !
!                                                                      !
!  Two vertical transformations are supported, z => z(x,y,s,t):        !
!                                                                      !
!  (1) Original transformation (Shchepetkin and McWilliams, 2005): In  !
!      ROMS since 1999 (version 1.8):                                  !
!                                                                      !
!        z(x,y,s,t) = Zo(x,y,s) + zeta(x,y,t) * [1 + Zo(x,y,s)/h(x,y)] !
!                                                                      !
!      where                                                           !
!                                                                      !
!        Zo(x,y,s) = hc * s + [h(x,y) - hc] * C(s)                     !
!                                                                      !
!        Zo(x,y,s) = 0         for s = 0,  C(s) = 0,  at the surface   !
!        Zo(x,y,s) = -h(x,y)   for s = -1, C(s) = -1, at the bottom    !
!                                                                      !
!  (2) New transformation: In UCLA-ROMS since 2005:                    !
!                                                                      !
!        z(x,y,s,t) = zeta(x,y,t) + [zeta(x,y,t) + h(x,y)] * Zo(x,y,s) !
!                                                                      !
!      where                                                           !
!                                                                      !
!        Zo(x,y,s) = [hc * s(k) + h(x,y) * C(k)] / [hc + h(x,y)]       !
!                                                                      !
!        Zo(x,y,s) = 0         for s = 0,  C(s) = 0,  at the surface   !
!        Zo(x,y,s) = -1        for s = -1, C(s) = -1, at the bottom    !
!                                                                      !
!      At the rest state, corresponding to zero free-surface, this     !
!      transformation yields the following unperturbed depths, zhat:   !
!                                                                      !
!        zhat = z(x,y,s,0) = h(x,y) * Zo(x,y,s)                        !
!                                                                      !
!             = h(x,y) * [hc * s(k) + h(x,y) * C(k)] / [hc + h(x,y)]   !
!                                                                      !
!      and                                                             !
!                                                                      !
!        d(zhat) = ds * h(x,y) * hc / [hc + h(x,y)]                    !
!                                                                      !
!      As a consequence, the uppermost grid box retains very little    !
!      dependency from bathymetry in the areas where hc << h(x,y),     !
!      that is deep areas. For example, if hc=250 m, and  h(x,y)       !
!      changes from 2000 to 6000 meters, the uppermost grid box        !
!      changes only by a factor of 1.08 (less than 10%).               !
!                                                                      !
!      Notice that:                                                    !
!                                                                      !
!      * Regardless of the design of C(s), transformation (2) behaves  !
!        like equally-spaced sigma-coordinates in shallow areas, where !
!        h(x,y) << hc.  This is advantageous because high vertical     !
!        resolution and associated CFL limitation is avoided in these  !
!        areas.                                                        !
!                                                                      !
!      * Near-surface refinement is close to geopotential coordinates  !
!        in deep areas (level thickness do not depend or weakly-depend !
!        on the bathymetry).  Contrarily,  near-bottom refinement is   !
!        like sigma-coordinates with thicknesses roughly proportional  !
!        to depth reducing high r-factors in these areas.              !
!                                                                      !
!                                                                      !
!  This generic transformation design facilitates numerous vertical    !
!  stretching functions, C(s).  These functions are set-up in this     !
!  routine in terms of several stretching parameters specified in      !
!  the standard input file.                                            !
!                                                                      !
!  C(s) vertical stretching function properties:                       !
!                                                                      !
!  * a nonlinear, monotonic function                                   !
!  * a continuous differentiable function, or                          !
!  * a piecewise function with smooth transition and differentiable    !
!  * must be constrained by -1 <= C(s) <= 0, with C(0)=0 at the        !
!    free-surface and C(-1)=-1 at the bottom (bathymetry).             !
!                                                                      !
!  References:                                                         !
!                                                                      !
!    Shchepetkin, A.F. and J.C. McWilliams, 2005: The regional oceanic !
!         modeling system (ROMS): a split-explicit, free-surface,      !
!         topography-following-coordinate oceanic model, Ocean         !
!         Modelling, 9, 347-404.                                       !
!                                                                      !
!    Song, Y. and D. Haidvogel, 1994: A semi-implicit ocean            !
!         circulation model using a generalized topography-            !
!         following coordinate system,  J.  Comp.  Physics,            !
!         115, 228-244.                                                !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_grid
      USE mod_iounits
      USE mod_scalars
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      integer :: k

      real(r8) :: Aweight, Bweight, Cweight, Cbot, Csur, Hscale
      real(r8) :: ds, exp_bot, exp_sur, sc_r, sc_w
      real(r8) :: cff, cff1, cff2, cff3
!
!-----------------------------------------------------------------------
!  Original vertical strectching function, Song and Haidvogel (1994).
!-----------------------------------------------------------------------
!
# if defined WET_DRY
        hc(ng)=MIN(MAX(hmin(ng),0.0_r8),Tcline(ng))
# else
        hc(ng)=MIN(hmin(ng),Tcline(ng))
# endif
        IF (Vstretching(ng).eq.1) THEN
!
!  This vertical stretching function is defined as:
!      
!      C(s) = (1 - b) * [SINH(s * a) / SINH(a)] +
!
!             b * [-0.5 + 0.5 * TANH(a * (s + 0.5)) / TANH(0.5 * a)]
!           
!  where the stretching parameters (a, b) are specify at input:
!
!         a = theta_s               0 < theta_s  <= 8
!         b = theta_b               0 <= theta_b <= 1
!
!  If theta_b=0, the refinement is surface intensified as theta_s is
!  increased.
!  If theta_b=1, the refinement is both bottom ans surface intensified
!  as theta_s is increased.
!

        IF (theta_s(ng).ne.0.0_r8) THEN
          cff1=1.0_r8/SINH(theta_s(ng))
          cff2=0.5_r8/TANH(0.5_r8*theta_s(ng))
        END IF
        SCALARS(ng)%sc_w(0)=-1.0_r8
        SCALARS(ng)%Cs_w(0)=-1.0_r8
        ds=1.0_r8/REAL(N(ng),r8)
        DO k=1,N(ng)
          SCALARS(ng)%sc_w(k)=ds*REAL(k-N(ng),r8)
          SCALARS(ng)%sc_r(k)=ds*(REAL(k-N(ng),r8)-0.5_r8)
          IF (theta_s(ng).ne.0.0_r8) THEN
            SCALARS(ng)%Cs_w(k)=(1.0_r8-theta_b(ng))*                   &
     &                          cff1*SINH(theta_s(ng)*                  &
     &                                    SCALARS(ng)%sc_w(k))+         &
     &                          theta_b(ng)*                            &
     &                          (cff2*TANH(theta_s(ng)*                 &
     &                                     (SCALARS(ng)%sc_w(k)+        &
     &                                      0.5_r8))-                   &
     &                           0.5_r8)
            SCALARS(ng)%Cs_r(k)=(1.0_r8-theta_b(ng))*                   &
     &                          cff1*SINH(theta_s(ng)*                  &
     &                                    SCALARS(ng)%sc_r(k))+         &
     &                          theta_b(ng)*                            &
     &                          (cff2*TANH(theta_s(ng)*                 &
     &                                     (SCALARS(ng)%sc_r(k)+        &
     &                                      0.5_r8))-                   &
     &                           0.5_r8)
          ELSE
            SCALARS(ng)%Cs_w(k)=SCALARS(ng)%sc_w(k)
            SCALARS(ng)%Cs_r(k)=SCALARS(ng)%sc_r(k)
          END IF
        END DO
!
!-----------------------------------------------------------------------
!  A. Shchepetkin new vertical stretching function.
!-----------------------------------------------------------------------
!
      ELSE IF (Vstretching(ng).eq.2) THEN
!
!  This vertical stretching function is defined, in the simplest form,
!  as:
!
!      C(s) = [1.0 - COSH(theta_s * s)] / [COSH(theta_s) - 1.0]
!
!  it is similar in meaning to the original vertical stretcing function
!  (Song and Haidvogel, 1994), but note that hyperbolic functions are
!  COSH, and not SINH.
!
!  Note that the above definition results in
!
!         -1 <= C(s) <= 0
!
!  as long as
!
!         -1 <= s <= 0
!
!  and, unlike in any previous definition
!
!         d[C(s)]/ds  -->  0      if  s -->  0
!
!  For the purpose of bottom boundary layer C(s) is further modified
!  to allow near-bottom refinement.  This is done by blending it with
!  another function.
!
    
          Aweight=1.0_r8
          Bweight=1.0_r8
          ds=1.0_r8/REAL(N(ng),r8)
!
          SCALARS(ng)%sc_w(N(ng))=0.0_r8
          SCALARS(ng)%Cs_w(N(ng))=0.0_r8
          DO k=N(ng)-1,1,-1
            sc_w=ds*REAL(k-N(ng),r8)
            SCALARS(ng)%sc_w(k)=sc_w
            IF (theta_s(ng).gt.0.0_r8) THEN
              Csur=(1.0_r8-COSH(theta_s(ng)*sc_w))/                     &
     &             (COSH(theta_s(ng))-1.0_r8)
              IF (theta_b(ng).gt.0.0_r8) THEN
                Cbot=SINH(theta_b(ng)*(sc_w+1.0_r8))/                   &
     &               SINH(theta_b(ng))-1.0_r8
                Cweight=(sc_w+1.0_r8)**Aweight*                         &
     &                  (1.0_r8+(Aweight/Bweight)*                      &
     &                          (1.0_r8-(sc_w+1.0_r8)**Bweight))
                SCALARS(ng)%Cs_w(k)=Cweight*Csur+(1.0_r8-Cweight)*Cbot
              ELSE
                SCALARS(ng)%Cs_w(k)=Csur
              END IF
            ELSE
              SCALARS(ng)%Cs_w(k)=sc_w
            END IF
          END DO
          SCALARS(ng)%sc_w(0)=-1.0_r8
          SCALARS(ng)%Cs_w(0)=-1.0_r8
!
          DO k=1,N(ng)
            sc_r=ds*(REAL(k-N(ng),r8)-0.5_r8)
            SCALARS(ng)%sc_r(k)=sc_r
            IF (theta_s(ng).gt.0.0_r8) THEN
              Csur=(1.0_r8-COSH(theta_s(ng)*sc_r))/                     &
     &             (COSH(theta_s(ng))-1.0_r8)
              IF (theta_b(ng).gt.0.0_r8) THEN
                Cbot=SINH(theta_b(ng)*(sc_r+1.0_r8))/                   &
     &               SINH(theta_b(ng))-1.0_r8
                Cweight=(sc_r+1.0_r8)**Aweight*                         &
     &                  (1.0_r8+(Aweight/Bweight)*                      &
     &                          (1.0_r8-(sc_r+1.0_r8)**Bweight))
                SCALARS(ng)%Cs_r(k)=Cweight*Csur+(1.0_r8-Cweight)*Cbot
              ELSE
                SCALARS(ng)%Cs_r(k)=Csur
              END IF
            ELSE
              SCALARS(ng)%Cs_r(k)=sc_r
            END IF
          END DO

!
!-----------------------------------------------------------------------
!  R. Geyer stretching function for high bottom boundary layer
!  resolution.
!-----------------------------------------------------------------------
!
      ELSE IF (Vstretching(ng).eq.3) THEN
!
!  This stretching function is intended for very shallow coastal
!  applications, like gravity sediment flows.
!
!  At the surface, C(s=0)=0
!
!      C(s) = - LOG(COSH(Hscale * ABS(s) ** alpha)) /
!               LOG(COSH(Hscale))
!
!  At the bottom, C(s=-1)=-1
!
!      C(s) = LOG(COSH(Hscale * (s + 1) ** beta)) /
!             LOG(COSH(Hscale)) - 1
!
!  where
!
!       Hscale : scale value for all hyperbolic functions
!                  Hscale = 3.0    set internally here
!        alpha : surface stretching exponent
!                  alpha = 0.65   minimal increase of surface resolution
!                          1.0    significant amplification
!         beta : bottoom stretching exponent
!                  beta  = 0.58   no amplification
!                          1.0    significant amplification
!                          3.0    super-high bottom resolution
!            s : stretched vertical coordinate, -1 <= s <= 0
!                  s(k) = (k-N)/N       k=0:N,    W-points  (s_w)
!                  s(k) = (k-N-0.5)/N   k=1:N,  RHO-points  (s_rho)
!
!  The stretching exponents (alpha, beta) are specify at input:
!
!         alpha = theta_s
!         beta  = theta_b
!
  
          exp_sur=theta_s(ng)
          exp_bot=theta_b(ng)
          Hscale=3.0_r8
          ds=1.0_r8/REAL(N(ng),r8)
!
          SCALARS(ng)%sc_w(N(ng))=0.0_r8
          SCALARS(ng)%Cs_w(N(ng))=0.0_r8
          DO k=N(ng)-1,1,-1
            sc_w=ds*REAL(k-N(ng),r8)
            SCALARS(ng)%sc_w(k)=sc_w
            Cbot= LOG(COSH(Hscale*(sc_w+1.0_r8)**exp_bot))/             &
     &            LOG(COSH(Hscale))-1.0_r8
            Csur=-LOG(COSH(Hscale*ABS(sc_w)**exp_sur))/                 &
     &            LOG(COSH(Hscale))
            Cweight=0.5_r8*(1.0_r8-TANH(Hscale*(sc_w+0.5_r8)))
            SCALARS(ng)%Cs_w(k)=Cweight*Cbot+(1.0_r8-Cweight)*Csur
          END DO
          SCALARS(ng)%sc_w(0)=-1.0_r8
          SCALARS(ng)%Cs_w(0)=-1.0_r8
!
          DO k=1,N(ng)
            sc_r=ds*(REAL(k-N(ng),r8)-0.5_r8)
            SCALARS(ng)%sc_r(k)=sc_r
            Cbot= LOG(COSH(Hscale*(sc_r+1.0_r8)**exp_bot))/             &
     &            LOG(COSH(Hscale))-1.0_r8
            Csur=-LOG(COSH(Hscale*ABS(sc_r)**exp_sur))/                 &
     &            LOG(COSH(Hscale))
            Cweight=0.5_r8*(1.0_r8-TANH(Hscale*(sc_r+0.5_r8)))
            SCALARS(ng)%Cs_r(k)=Cweight*Cbot+(1.0_r8-Cweight)*Csur
          END DO

      END IF
!
!-----------------------------------------------------------------------
!  Report information about vertical transformation.
!-----------------------------------------------------------------------
!
      IF (Master.and.LwrtInfo(ng)) THEN
        WRITE (stdout,20)
        DO k=N(ng),0,-1
          IF (Vstretching(ng).eq.2) THEN
            cff=0.5_r8*(hmax(ng)+hmin(ng))
            cff1=hc(ng)  *(SCALARS(ng)%sc_w(k)*hc(ng)+                  &
     &                     SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+hc(ng))
            cff2=cff     *(SCALARS(ng)%sc_w(k)*cff+                     &
     &                     SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+cff)
            cff3=hmax(ng)*(SCALARS(ng)%sc_w(k)*hmax(ng)+                &
     &                     SCALARS(ng)%Cs_w(k)*hc(ng))/(hc(ng)+hmax(ng))
          ELSE
            cff1=SCALARS(ng)%sc_w(k)*hc(ng)+                            &
     &           (hmin(ng)-hc(ng))*SCALARS(ng)%Cs_w(k)
            cff2=SCALARS(ng)%sc_w(k)*hc(ng)+                            &
     &           (0.5*(hmin(ng)+hmax(ng))-hc(ng))*SCALARS(ng)%Cs_w(k)
            cff3=SCALARS(ng)%sc_w(k)*hc(ng)+                            &
     &           (hmax(ng)-hc(ng))*SCALARS(ng)%Cs_w(k)
          END IF
          WRITE (stdout,30) k, SCALARS(ng)%sc_w(k),                     &
     &                         SCALARS(ng)%Cs_w(k), cff1, cff2, cff3
        END DO
      END IF

  10  FORMAT (/,' SET_SCOORD - Specified critical depth parameter,  ',  &
     &        '  hc =',f7.3,/,14x,                                      &
     &        'exceeds minimum unmasked bathymetry, hmin = ',f7.3)
  20  FORMAT (/,' Vertical S-coordinate System: ',/,/,                  &
     &          ' level   S-coord     Cs-curve',10x,                    &
     &          'at_hmin  over_slope     at_hmax',/)
  30  FORMAT (i6,2f12.7,4x,3f12.3)

      RETURN
      END SUBROUTINE set_scoord
#else
      SUBROUTINE set_scoord
      RETURN
      END SUBROUTINE set_scoord
#endif

