#include "cppdefs.h" #ifdef SOLVE3D SUBROUTINE set_scoord (ng) ! !svn $Id: set_scoord.F 374 2009-07-24 18:54:26Z 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 ! !----------------------------------------------------------------------- ! Set thickness controlling vertical coordinate stretching. !----------------------------------------------------------------------- ! ! Set hc <= hmin, in the original formulation (Vtransform=1) to avoid ! [h(x,y)-hc] to be negative which results in dz/ds to be negative. ! Notice that this restriction is REMOVED in the new transformation ! (Vtransform=2): hc can be any value. It works for both hc < hmin ! and hc > hmin. ! IF (Vtransform(ng).eq.1) THEN # if defined WET_DRY hc(ng)=MIN(MAX(hmin(ng),0.0_r8),Tcline(ng)) # else hc(ng)=MIN(hmin(ng),Tcline(ng)) # endif ELSE IF (Vtransform(ng).eq.2) THEN hc(ng)=Tcline(ng) END IF ! !----------------------------------------------------------------------- ! Original vertical strectching function, Song and Haidvogel (1994). !----------------------------------------------------------------------- ! 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,10) DO k=N(ng),0,-1 IF (Vstretching(ng).eq.2) THEN cff=0.5_r8*(hmax(ng)+hmin(ng)) cff1=hmin(ng)*(SCALARS(ng)%sc_w(k)*hc(ng)+ & & SCALARS(ng)%Cs_w(k)*hmin(ng))/ & & (hc(ng)+hmin(ng)) cff2=cff *(SCALARS(ng)%sc_w(k)*hc(ng)+ & & SCALARS(ng)%Cs_w(k)*cff)/ & & (hc(ng)+cff) cff3=hmax(ng)*(SCALARS(ng)%sc_w(k)*hc(ng)+ & & SCALARS(ng)%Cs_w(k)*hmax(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,20) k, SCALARS(ng)%sc_w(k), & & SCALARS(ng)%Cs_w(k), cff1, cff2, cff3 END DO END IF 10 FORMAT (/,' Vertical S-coordinate System: ',/,/, & & ' level S-coord Cs-curve',10x, & & 'at_hmin over_slope at_hmax',/) 20 FORMAT (i6,2f12.7,4x,3f12.3) RETURN END SUBROUTINE set_scoord #else SUBROUTINE set_scoord RETURN END SUBROUTINE set_scoord #endif