New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3600 for branches – NEMO

Changeset 3600 for branches


Ignore:
Timestamp:
2012-11-19T14:59:22+01:00 (11 years ago)
Author:
rfurner
Message:

Changes from branch dev_r3435_UKMO7_SCOORDS revision 3435 to 3507 merged in

Location:
branches/2012/dev_UKMO_2012
Files:
14 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/2012/dev_UKMO_2012/DOC/TexFiles/Biblio/Biblio.bib

    r3294 r3600  
    25242524} 
    25252525 
     2526@ARTICLE{Siddorn_Furner_OM12, 
     2527  author = {J. Siddorn and R. Furner}, 
     2528  title = {An analytical stretching function that combines the best attributes of geopotential and terrain-following vertical coordinates}, 
     2529  journal = OM, 
     2530  year = {2012}, 
     2531  pages = {submitted}, 
     2532} 
    25262533@ARTICLE{Simmons_al_OM04, 
    25272534  author = {H. L. Simmons and S. R. Jayne and L. C. {St. Laurent} and A. J. Weaver}, 

  • branches/2012/dev_UKMO_2012/DOC/TexFiles/Chapters/Chap_DOM.tex

    r3294 r3600  
    502502time-variation of the free surface so that the transformation is time dependent:  
    503503$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step  
    504 bathymetry or $s$-coordinate (hybride and partial step coordinates have not  
     504bathymetry or $s$-coordinate (hybrid and partial step coordinates have not  
    505505yet been tested in NEMO v2.3).  
    506506 
    … …  
    743743levels are defined from the product of a depth field and either a stretching  
    744744function or its derivative, respectively: 
     745 
    745746\begin{equation} \label{DOM_sco_ana} 
    746747\begin{split} 
    … …  
    749750\end{split} 
    750751\end{equation} 
     752 
    751753where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point  
    752754location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea  
    753755surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean  
    754756depth, since a mixed step-like and bottom-following representation of the  
    755 topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided  
    756 (\rou{zgr\_sco} routine, see \mdl{domzgr}) $h$ is a smooth envelope bathymetry  
    757 and steps are used to represent sharp bathymetric gradients. 
    758  
    759 A new flexible stretching function, modified from \citet{Song_Haidvogel_JCP94} is provided as an example: 
     757topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f). 
     758The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 
     759 
     760Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application. 
     761A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is provided as an example: 
     762 
    760763\begin{equation} \label{DOM_sco_function} 
    761764\begin{split} 
    762 z  &= h_c +( h-h_c)\;c s)  \\ 
     765z  &= h_c +( h-h_c)\;c s   \\ 
    763766c(s)  &=  \frac{ \left[   \tanh{ \left( \theta \, (s+b) \right)}  
    764767               - \tanh{ \left(  \theta \, b      \right)}  \right]} 
    … …  
    766769\end{split} 
    767770\end{equation} 
    768 where $h_c$ is the thermocline depth and $\theta$ and $b$ are the surface and  
     771 
     772 
     773%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     774\begin{figure}[!ht]    \begin{center} 
     775\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf} 
     776\caption{  \label{Fig_sco_function}    
     777Examples of the stretching function applied to a seamount; from left to right:  
     778surface, surface and bottom, and bottom intensified resolutions} 
     779\end{center}   \end{figure} 
     780%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     781 
     782where $h_c$ is the critical depth (\np{rn\_hc}) total depth at which the coordinate transitions from pure $\sigma$ to the stretched coordinate,  and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and  
    769783bottom control parameters such that $0\leqslant \theta \leqslant 20$, and  
    770784$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom  
    771785increase of the vertical resolution (Fig.~\ref{Fig_sco_function}). 
    772786 
    773 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    774 \begin{figure}[!tb]    \begin{center} 
    775 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf} 
    776 \caption{  \label{Fig_sco_function}    
    777 Examples of the stretching function applied to a sea mont; from left to right:  
    778 surface, surface and bottom, and bottom intensified resolutions} 
    779 \end{center}   \end{figure} 
    780 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     787Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that: 
     788 
     789\begin{equation} 
     790z = \gamma\left(h+\zeta\right) \quad \text{ with } \quad 0 \leq \gamma \leq 1 
     791\label{eq:z} 
     792\end{equation} 
     793 
     794The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 
     795 
     796\begin{equation} \label{DOM_gamma_deriv} 
     797\gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right) 
     798\end{equation} 
     799 
     800Where: 
     801\begin{equation} \label{DOM_gamma} 
     802f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} 
     803\end{equation} 
     804 
     805This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth: 
     806 
     807\begin{equation} \label{DOM_zb} 
     808Z_b= h a + b 
     809\end{equation} 
     810 
     811where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. 
     812 
     813%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     814\begin{figure}[!ht] 
     815   \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/FIG_DOM_compare_coordinates_surface.pdf} 
     816        \caption{A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.} 
     817    \label{fig_compare_coordinates_surface} 
     818\end{figure} 
     819%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     820 
     821This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell depths in real space. This is not to be confused with the hybrid schemes that superimpose geopotential coordinates on terrain following coordinates thus creating a non-analytical vertical coordinate that therefore may suffer from large gradients in the vertical resolutions. This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. 
     822 
     823As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than the critical depth $h_c$. In this example two options are available in depths shallower than $h_c$, with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$. 
     824 
     825Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as large slopes lead to hydrostatic consistency. A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, and is output as part of the model mesh file at the start of the run. 
    781826 
    782827% ------------------------------------------------------------------------------------------------------------- 

  • branches/2012/dev_UKMO_2012/DOC/TexFiles/Namelist/namzgr_sco

    r3294 r3600  
    22&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    33!----------------------------------------------------------------------- 
    4    rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m) 
    5    rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    6    rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20) 
    7    rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1) 
    8    rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1) 
    9    ln_s_sigma  = .false.   !  hybrid s-sigma coordinates 
    10    rn_bb       =    0.8    !  stretching with s-sigma 
    11    rn_hc       =  150.0    !  critical depth with s-sigma  
     4   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     5   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     6   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     7                           !  stretching coefficients for all functions 
     8   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     9   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     10   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     11                        !!!!!!!  Envelop bathymetry 
     12   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     13                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
     14   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
     15   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     16                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     17   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     18   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     19   rn_zs       =    1.0    !  depth of surface grid box 
     20                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     21   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     22   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     23                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     24   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    1225/ 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/AMM12/EXP00/namelist

    r3309 r3600  
    2727   cn_exp      =  "AMM12"  !  experience name  
    2828   nn_it000    =       1   !  first time step 
    29    nn_itend    =     576   !  last  time step (std 1 day = 576) 
     29   nn_itend    =    2880   !  last  time step (std 1 day = 288) 
    3030   nn_date0    =  20070101 !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    3131   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    32    ln_rstart   =  .false.  !  start from rest (F) or from a restart file (T) 
     32   ln_rstart   =  .true.  !  start from rest (F) or from a restart file (T) 
    3333   nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T 
    3434                           !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    … …  
    3737   cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input) 
    3838   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
    39    nn_istate   =       0   !  output the initial state (1) or not (0) 
    40    nn_stock    =     576   !  frequency of creation of a restart file (modulo referenced to 1) 
    41    nn_write    =      12   !  frequency of write in the output file   (modulo referenced to nit000) 
     39   nn_istate   =       1   !  output the initial state (1) or not (0) 
     40   nn_stock    =     2880  !  frequency of creation of a restart file (modulo referenced to 1) 
     41   nn_write    =     144   !  frequency of write in the output file   (modulo referenced to nit000)  
    4242   ln_dimgnnn  = .false.   !  DIMG file format: 1 file for all processors (F) or by processor (T) 
    4343   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%) 
    … …  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =   10.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 7000.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     67   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     70                           !  stretching coefficients for all functions 
     71   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     72   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     73   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     74                        !!!!!!!  Envelop bathymetry 
     75   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     76                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
    6977   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
    70    rn_thetb    =    1.00   !  bottom control parameter  (0<=thetb<= 1) 
    71    rn_rmax     =    0.30   !  maximum cut-off r-value allowed (0<r_max<1) 
    72    ln_s_sigma  = .true.    !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
     78   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     79                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     80   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     81   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     82   rn_zs       =    1.0    !  depth of surface grid box 
     83                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     84   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     85   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     86                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     87   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    7588/ 
    7689!----------------------------------------------------------------------- 
    … …  
    7992   nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file 
    8093   nn_closea    =   0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA) 
    81    nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
     94   nn_msh      =    1      !  create (=1) a mesh file or not (=0) 
    8295   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
    8396   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
    8497   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 
    8598                           ! 
    86    rn_rdt      =  150.     !  time step for the dynamics (and tracer if nn_acc=0) 
     99   rn_rdt      =  300.     !  time step for the dynamics (and tracer if nn_acc=0) 
    87100   nn_baro     =   30      !  number of barotropic time step            ("key_dynspg_ts") 
    88101   rn_atfp     =    0.1    !  asselin time filter parameter 
    … …  
    394407    cn_mask_file = ''                     !  name of mask file (if ln_mask_file=.TRUE.) 
    395408    nn_dyn2d      =  2                    !  boundary conditions for barotropic fields 
    396     nn_dyn2d_dta  =  2                    !  = 0, bdy data are equal to the initial state 
     409    nn_dyn2d_dta  =  3                    !  = 0, bdy data are equal to the initial state 
    397410                                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    398411                                          !  = 2, use tidal harmonic forcing data from files 
    … …  
    402415                           !  = 1, bdy data are read in 'bdydata   .nc' files 
    403416    nn_tra        =  1                    !  boundary conditions for T and S 
    404     nn_tra_dta    =  0                    !  = 0, bdy data are equal to the initial state 
     417    nn_tra_dta    =  1                    !  = 0, bdy data are equal to the initial state 
    405418                           !  = 1, bdy data are read in 'bdydata   .nc' files 
    406419    nn_rimwidth  = 10                      !  width of the relaxation zone 
    … …  
    596609   ln_hpg_zco  = .false.   !  z-coordinate - full steps                    
    597610   ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation) 
    598    ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation) 
     611   ln_hpg_sco  = .false.    !  s-coordinate (standard jacobian formulation) 
    599612   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    600    ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
     613   ln_hpg_prj  = .true.   !  s-coordinate (Pressure Jacobian scheme) 
    601614   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
    602615                                 !           centered      time scheme  (F) 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist

    r3309 r3600  
    2727   cn_exp      =  "AMM12"  !  experience name  
    2828   nn_it000    =       1   !  first time step 
    29    nn_itend    =     576   !  last  time step (std 1 day = 576) 
     29   nn_itend    =    2880   !  last  time step (std 1 day = 288) 
    3030   nn_date0    =  20070101 !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    3131   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    32    ln_rstart   =  .false.  !  start from rest (F) or from a restart file (T) 
     32   ln_rstart   =  .true.  !  start from rest (F) or from a restart file (T) 
    3333   nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T 
    3434                           !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    … …  
    3737   cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input) 
    3838   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
    39    nn_istate   =       0   !  output the initial state (1) or not (0) 
    40    nn_stock    =     576   !  frequency of creation of a restart file (modulo referenced to 1) 
    41    nn_write    =      12   !  frequency of write in the output file   (modulo referenced to nit000) 
     39   nn_istate   =       1   !  output the initial state (1) or not (0) 
     40   nn_stock    =     2880  !  frequency of creation of a restart file (modulo referenced to 1) 
     41   nn_write    =     144   !  frequency of write in the output file   (modulo referenced to nit000)  
    4242   ln_dimgnnn  = .false.   !  DIMG file format: 1 file for all processors (F) or by processor (T) 
    4343   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%) 
    … …  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =   10.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 7000.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     67   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     70                           !  stretching coefficients for all functions 
     71   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     72   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     73   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     74                        !!!!!!!  Envelop bathymetry 
     75   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     76                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
    6977   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
    70    rn_thetb    =    1.00   !  bottom control parameter  (0<=thetb<= 1) 
    71    rn_rmax     =    0.30   !  maximum cut-off r-value allowed (0<r_max<1) 
    72    ln_s_sigma  = .true.    !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
    75 / 
     78   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     79                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     80   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     81   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     82   rn_zs       =    1.0    !  depth of surface grid box 
     83                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     84   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     85   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     86                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     87   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
     88/ 
     89 
    7690!----------------------------------------------------------------------- 
    7791&namdom        !   space and time domain (bathymetry, mesh, timestep) 
    … …  
    7993   nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file 
    8094   nn_closea    =   0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA) 
    81    nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
     95   nn_msh      =    1      !  create (=1) a mesh file or not (=0) 
    8296   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
    8397   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
    8498   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 
    8599                           ! 
    86    rn_rdt      =  150.     !  time step for the dynamics (and tracer if nn_acc=0) 
     100   rn_rdt      =  300.     !  time step for the dynamics (and tracer if nn_acc=0) 
    87101   nn_baro     =   30      !  number of barotropic time step            ("key_dynspg_ts") 
    88102   rn_atfp     =    0.1    !  asselin time filter parameter 
    … …  
    102116   ! 
    103117   cn_dir        = './'     !  root directory for the location of the runoff files 
    104    ln_tsd_init   = .false.  !  Initialisation of ocean T & S with T &S input data (T) or not (F) 
    105    ln_tsd_tradmp = .false.  !  damping of ocean T & S toward T &S input data (T) or not (F) 
     118   ln_tsd_init   = .false.   !  Initialisation of ocean T & S with T &S input data (T) or not (F) 
     119   ln_tsd_tradmp = .false.   !  damping of ocean T & S toward T &S input data (T) or not (F) 
    106120/ 
    107121!!====================================================================== 
    … …  
    302316                           !  or to SSS only (=1) or no damping term (=0) 
    303317   rn_dqdt     =   -40.    !  magnitude of the retroaction on temperature   [W/m2/K] 
    304    rn_deds     =  -166.67  !  magnitude of the damping on salinity   [mm/day] 
     318   rn_deds     =  -27.7    !  magnitude of the damping on salinity   [mm/day] 
    305319   ln_sssr_bnd =   .true.  !  flag to bound erp term (associated with nn_sssr=2) 
    306320   rn_sssr_bnd =   4.e0    !  ABS(Max/Min) value of the damping erp term [mm/day] 
    … …  
    596610   ln_hpg_zco  = .false.   !  z-coordinate - full steps                    
    597611   ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation) 
    598    ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation) 
     612   ln_hpg_sco  = .false.    !  s-coordinate (standard jacobian formulation) 
    599613   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    600    ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
     614   ln_hpg_prj  = .true.   !  s-coordinate (Pressure Jacobian scheme) 
    601615   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
    602616                                 !           centered      time scheme  (F) 
    … …  
    875889                           !  0 < n : debug section number n 
    876890/ 
    877  
    878891!!====================================================================== 
    879892!!            ***  Observation & Assimilation namelists *** 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/GYRE/EXP00/namelist

    r3306 r3600  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    69    rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20) 
    70    rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1) 
    71    rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1) 
    72    ln_s_sigma  = .false.   !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
     67   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     70                           !  stretching coefficients for all functions 
     71   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     72   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     73   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     74                        !!!!!!!  Envelop bathymetry 
     75   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     76                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
     77   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
     78   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     79                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     80   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     81   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     82   rn_zs       =    1.0    !  depth of surface grid box 
     83                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     84   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     85   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     86                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     87   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    7588/ 
    7689!----------------------------------------------------------------------- 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist

    r3306 r3600  
    6666&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6767!----------------------------------------------------------------------- 
    68    rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m) 
    69    rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    70    rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20) 
    71    rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1) 
    72    rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1) 
    73    ln_s_sigma  = .false.   !  hybrid s-sigma coordinates 
    74    rn_bb       =    0.8    !  stretching with s-sigma 
    75    rn_hc       =  150.0    !  critical depth with s-sigma  
     68   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     69   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     70   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     71                           !  stretching coefficients for all functions 
     72   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     73   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     74   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     75                        !!!!!!!  Envelop bathymetry 
     76   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     77                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
     78   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
     79   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     80                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     81   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     82   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     83   rn_zs       =    1.0    !  depth of surface grid box 
     84                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     85   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     86   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     87                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     88   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    7689/ 
    7790!----------------------------------------------------------------------- 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3306 r3600  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    69    rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20) 
    70    rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1) 
    71    rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1) 
    72    ln_s_sigma  = .false.   !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
     67   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     70                           !  stretching coefficients for all functions 
     71   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     72   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     73   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     74                        !!!!!!!  Envelop bathymetry 
     75   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     76                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
     77   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
     78   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     79                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     80   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     81   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     82   rn_zs       =    1.0    !  depth of surface grid box 
     83                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     84   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     85   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     86                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     87   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    7588/ 
    7689!----------------------------------------------------------------------- 

  • branches/2012/dev_UKMO_2012/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist

    r3306 r3600  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
    69    rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20) 
    70    rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1) 
    71    rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1) 
    72    ln_s_sigma  = .false.   !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
     67   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
     70                           !  stretching coefficients for all functions 
     71   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     72   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     73   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates 
     74                        !!!!!!!  Envelop bathymetry 
     75   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     76                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.) 
     77   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
     78   rn_bb       =    0.8    !  stretching with SH94 s-sigma    
     79                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.) 
     80   rn_alpha    =    4.4    !  stretching with SF12 s-sigma 
     81   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     82   rn_zs       =    1.0    !  depth of surface grid box 
     83                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     84   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     85   rn_zb_b     =   -0.2    !  offset for calculating Zb 
     86                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
     87   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)  
    7588/ 
    7689!----------------------------------------------------------------------- 

  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r3421 r3600  
    174174   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff     !: interface depth between stretching at  V--F 
    175175   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m) 
     176   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1              !: Maximum grid stiffness ratio 
    176177 
    177178   !!---------------------------------------------------------------------- 
    … …  
    294295         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    295296         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    296          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) ) 
     297         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
    297298 
    298299      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     & 

  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r3421 r3600  
    3636   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
    3737   USE timing          ! Timing 
     38   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    3839 
    3940   IMPLICIT NONE 
    … …  
    8485                             CALL dom_zgr      ! Vertical mesh and bathymetry 
    8586                             CALL dom_msk      ! Masks 
     87      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
    8688      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh 
    8789      ! 
    … …  
    322324   END SUBROUTINE dom_ctl 
    323325 
     326   SUBROUTINE dom_stiff 
     327      !!---------------------------------------------------------------------- 
     328      !!                  ***  ROUTINE dom_stiff  *** 
     329      !!                      
     330      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency 
     331      !! 
     332      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio 
     333      !!                Save the maximum in the vertical direction 
     334      !!                (this number is only relevant in s-coordinates) 
     335      !! 
     336      !!                Haney, R. L., 1991: On the pressure gradient force 
     337      !!                over steep topography in sigma coordinate ocean models.  
     338      !!                J. Phys. Oceanogr., 21, 610???619. 
     339      !!---------------------------------------------------------------------- 
     340      INTEGER  ::   ji, jj, jk  
     341      REAL(wp) ::   zrxmax 
     342      REAL(wp), DIMENSION(4) :: zr1 
     343      !!---------------------------------------------------------------------- 
     344      rx1(:,:) = 0.e0 
     345      zrxmax   = 0.e0 
     346      zr1(:)   = 0.e0 
     347       
     348      DO ji = 2, jpim1 
     349         DO jj = 2, jpjm1 
     350            DO jk = 1, jpkm1 
     351               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  &  
     352                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) & 
     353                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  & 
     354                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) ) 
     355               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  & 
     356                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) & 
     357                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  & 
     358                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) ) 
     359               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  & 
     360                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) & 
     361                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  & 
     362                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) ) 
     363               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  & 
     364                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) & 
     365                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  & 
     366                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) ) 
     367               zrxmax = MAXVAL(zr1(1:4)) 
     368               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) 
     369            END DO 
     370         END DO 
     371      END DO 
     372 
     373      CALL lbc_lnk( rx1, 'T', 1. ) 
     374 
     375      zrxmax = MAXVAL(rx1) 
     376 
     377      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain 
     378 
     379      IF(lwp) THEN 
     380         WRITE(numout,*) 
     381         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 
     382         WRITE(numout,*) '~~~~~~~~~' 
     383      ENDIF 
     384 
     385   END SUBROUTINE dom_stiff 
     386 
     387 
     388 
    324389   !!====================================================================== 
    325390END MODULE domain 

  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r3294 r3600  
    172172             
    173173      IF( ln_sco ) THEN                                         ! s-coordinate 
    174          CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth 
    175          CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )  
     174         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 
     175         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
    176176         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 
    177177         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 
    … …  
    187187         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    188188         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    189          ! 
    190          CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system 
    191          CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 
     189         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio 
     190         ! 
     191         CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept )    !    ! stretched system 
     192         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 
    192193      ENDIF 
    193194       

  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r3421 r3600  
    1515   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
    1616   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
     17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn adn Furner stretching functio 
    1718   !!---------------------------------------------------------------------- 
    1819 
    … …  
    2829   !!       zgr_sco      : s-coordinate 
    2930   !!       fssig        : sigma coordinate non-dimensional function 
     31   !!       fgamma       : Siddorn and Furner stretching function 
    3032   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
    3133   !!--------------------------------------------------------------------- 
    … …  
    4749 
    4850   !                                       !!* Namelist namzgr_sco * 
     51   LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
     52   LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 
     53   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
     54   ! 
    4955   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
    5056   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     57   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
     58   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates 
     59   ! Song and Haidvogel 1994 stretching parameters 
    5160   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20) 
    5261   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1) 
    53    REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    54    LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 
    55    REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching 
     62   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter  
    5663   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    57    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates 
     64   ! Siddorn and Furner stretching parameters 
     65   REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 
     66   REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord 
     67   REAL(wp) ::   rn_zs       =    1.0_wp    !  depth of surface grid box 
     68                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     69   REAL(wp) ::   rn_zb_a     =    0.024_wp  !  bathymetry scaling factor for calculating Zb 
     70   REAL(wp) ::   rn_zb_b     =   -0.2_wp    !  offset for calculating Zb 
    5871 
    5972  !! * Substitutions 
    … …  
    10341047   END SUBROUTINE zgr_zps 
    10351048 
    1036  
    1037    FUNCTION fssig( pk ) RESULT( pf ) 
    1038       !!---------------------------------------------------------------------- 
    1039       !!                 ***  ROUTINE eos_init  *** 
    1040       !!        
    1041       !! ** Purpose :   provide the analytical function in s-coordinate 
    1042       !!           
    1043       !! ** Method  :   the function provide the non-dimensional position of 
    1044       !!                T and W (i.e. between 0 and 1) 
    1045       !!                T-points at integer values (between 1 and jpk) 
    1046       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1047       !!---------------------------------------------------------------------- 
    1048       REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
    1049       REAL(wp)             ::   pf   ! sigma value 
    1050       !!---------------------------------------------------------------------- 
    1051       ! 
    1052       pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
    1053          &     - TANH( rn_thetb * rn_theta                                )  )   & 
    1054          & * (   COSH( rn_theta                           )                      & 
    1055          &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
    1056          & / ( 2._wp * SINH( rn_theta ) ) 
    1057       ! 
    1058    END FUNCTION fssig 
    1059  
    1060  
    1061    FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
    1062       !!---------------------------------------------------------------------- 
    1063       !!                 ***  ROUTINE eos_init  *** 
    1064       !! 
    1065       !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
    1066       !! 
    1067       !! ** Method  :   the function provides the non-dimensional position of 
    1068       !!                T and W (i.e. between 0 and 1) 
    1069       !!                T-points at integer values (between 1 and jpk) 
    1070       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1071       !!---------------------------------------------------------------------- 
    1072       REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
    1073       REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
    1074       REAL(wp)             ::   pf1   ! sigma value 
    1075       !!---------------------------------------------------------------------- 
    1076       ! 
    1077       IF ( rn_theta == 0 ) then      ! uniform sigma 
    1078          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
    1079       ELSE                        ! stretched sigma 
    1080          pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
    1081             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
    1082             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
    1083       ENDIF 
    1084       ! 
    1085    END FUNCTION fssig1 
    1086  
    1087  
    10881049   SUBROUTINE zgr_sco 
    10891050      !!---------------------------------------------------------------------- 
    … …  
    11101071      !!            esigt(k) = fsdsig(k    ) 
    11111072      !!            esigw(k) = fsdsig(k-0.5) 
    1112       !!      This routine is given as an example, it must be modified 
    1113       !!      following the user s desiderata. nevertheless, the output as 
     1073      !!      Three options for stretching are give, and they can be modified 
     1074      !!      following the users requirements. Nevertheless, the output as 
    11141075      !!      well as the way to compute the model levels and scale factors 
    1115       !!      must be respected in order to insure second order a!!uracy 
     1076      !!      must be respected in order to insure second order accuracy 
    11161077      !!      schemes. 
    11171078      !! 
    1118       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1079      !!      The three methods for stretching available are: 
     1080      !!  
     1081      !!           s_sh94 (Song and Haidvogel 1994) 
     1082      !!                a sinh/tanh function that allows sigma and stretched sigma 
     1083      !! 
     1084      !!           s_sf12 (Siddorn and Furner 2012?) 
     1085      !!                allows the maintenance of fixed surface and or 
     1086      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1087      !!                within an analytically derived stretched S-coordinate framework. 
     1088      !!  
     1089      !!          s_tanh  (Madec et al 1996) 
     1090      !!                a cosh/tanh function that gives stretched coordinates         
     1091      !! 
    11191092      !!---------------------------------------------------------------------- 
    11201093      ! 
    11211094      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument 
    11221095      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    1123       REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars 
     1096      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    11241097      ! 
    11251098      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
    1126       REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
    1127       REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
    1128  
    1129       NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 
    1130       !!---------------------------------------------------------------------- 
     1099 
     1100      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     1101                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
     1102     !!---------------------------------------------------------------------- 
    11311103      ! 
    11321104      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11331105      ! 
    11341106      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
    1135       CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1136       CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    11371107      ! 
    11381108      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters 
    … …  
    11441114         WRITE(numout,*) '~~~~~~~~~~~' 
    11451115         WRITE(numout,*) '   Namelist namzgr_sco' 
    1146          WRITE(numout,*) '      sigma-stretching coeffs ' 
    1147          WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max 
    1148          WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min 
    1149          WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta 
    1150          WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb 
    1151          WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax 
    1152          WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma 
    1153          WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb 
    1154          WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc 
    1155       ENDIF 
    1156  
    1157       gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
    1158       esigt3  = 0._wp   ;   esigw3  = 0._wp  
    1159       esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
    1160       esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1116         WRITE(numout,*) '     stretching coeffs ' 
     1117         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max 
     1118         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min 
     1119         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc 
     1120         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax 
     1121         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94 
     1122         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients' 
     1123         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta 
     1124         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb 
     1125         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
     1126         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12 
     1127         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients' 
     1128         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     1129         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
     1130         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
     1131         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
     1132         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
     1133         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
     1134      ENDIF 
    11611135 
    11621136      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate 
    … …  
    13521326      ! non-dimensional "sigma" for model level depth at w- and t-levels 
    13531327 
    1354       IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths 
    1355          !                         ! below rn_hc, with uniform sigma in shallower waters 
    1356          DO ji = 1, jpi 
    1357             DO jj = 1, jpj 
    1358  
    1359                IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    1360                   DO jk = 1, jpk 
    1361                      gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
    1362                      gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    1363                   END DO 
    1364                ELSE ! shallow water, uniform sigma 
    1365                   DO jk = 1, jpk 
    1366                      gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    1367                      gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
    1368                   END DO 
    1369                ENDIF 
    1370                IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 
    1371                ! 
    1372                DO jk = 1, jpkm1 
    1373                   esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1374                   esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
    1375                END DO 
    1376                esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
    1377                esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
    1378                ! 
    1379                ! Coefficients for vertical depth as the sum of e3w scale factors 
    1380                gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
    1381                DO jk = 2, jpk 
    1382                   gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
    1383                END DO 
    1384                ! 
    1385                DO jk = 1, jpk 
    1386                   zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1387                   zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1388                   gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    1389                   gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    1390                   gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    1391                END DO 
    1392                ! 
    1393             END DO   ! for all jj's 
    1394          END DO    ! for all ji's 
    1395  
    1396          DO ji = 1, jpim1 
    1397             DO jj = 1, jpjm1 
    1398                DO jk = 1, jpk 
    1399                   esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
    1400                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1401                   esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
    1402                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1403                   esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
    1404                      &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
    1405                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1406                   esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
    1407                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1408                   esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
    1409                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1410                   ! 
    1411                   e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1412                   e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1413                   e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1414                   e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1415                   ! 
    1416                   e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1417                   e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1418                   e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1419                END DO 
    1420             END DO 
    1421          END DO 
    1422  
    1423          CALL lbc_lnk( e3t , 'T', 1._wp ) 
    1424          CALL lbc_lnk( e3u , 'U', 1._wp ) 
    1425          CALL lbc_lnk( e3v , 'V', 1._wp ) 
    1426          CALL lbc_lnk( e3f , 'F', 1._wp ) 
    1427          CALL lbc_lnk( e3w , 'W', 1._wp ) 
    1428          CALL lbc_lnk( e3uw, 'U', 1._wp ) 
    1429          CALL lbc_lnk( e3vw, 'V', 1._wp ) 
    1430  
    1431          ! 
    1432       ELSE   ! not ln_s_sigma 
    1433          ! 
    1434          DO jk = 1, jpk 
    1435            gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
    1436            gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    1437          END DO 
    1438          IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    1439          ! 
    1440          ! Coefficients for vertical scale factors at w-, t- levels 
    1441 !!gm bug :  define it from analytical function, not like juste bellow.... 
    1442 !!gm        or betteroffer the 2 possibilities.... 
    1443          DO jk = 1, jpkm1 
    1444             esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
    1445             esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
    1446          END DO 
    1447          esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
    1448          esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
    1449  
    1450 !!gm  original form 
    1451 !!org DO jk = 1, jpk 
    1452 !!org    esigt(jk)=fsdsig( FLOAT(jk)     ) 
    1453 !!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 
    1454 !!org END DO 
    1455 !!gm 
    1456          ! 
    1457          ! Coefficients for vertical depth as the sum of e3w scale factors 
    1458          gsi3w(1) = 0.5_wp * esigw(1) 
    1459          DO jk = 2, jpk 
    1460             gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
    1461          END DO 
    1462 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    1463          DO jk = 1, jpk 
    1464             zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1465             zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1466             gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
    1467             gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
    1468             gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
    1469          END DO 
    1470 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    1471          DO jj = 1, jpj 
    1472             DO ji = 1, jpi 
    1473                DO jk = 1, jpk 
    1474                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1475                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1476                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1477                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
    1478                  ! 
    1479                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1480                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1481                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1482                END DO 
    1483             END DO 
    1484          END DO 
    1485          ! 
    1486       ENDIF ! ln_s_sigma 
    1487  
    1488  
     1328 
     1329!======================================================================== 
     1330! Song and Haidvogel  1994 (ln_s_sh94=T) 
     1331! Siddorn and Furner 2012 (ln_sf12=T) 
     1332! or  tanh function       (both false)                     
     1333!======================================================================== 
     1334      IF      ( ln_s_sh94 ) THEN  
     1335                           CALL s_sh94() 
     1336      ELSE IF ( ln_s_sf12 ) THEN 
     1337                           CALL s_sf12() 
     1338      ELSE                  
     1339                           CALL s_tanh() 
     1340      ENDIF  
     1341 
     1342      CALL lbc_lnk( e3t , 'T', 1._wp ) 
     1343      CALL lbc_lnk( e3u , 'U', 1._wp ) 
     1344      CALL lbc_lnk( e3v , 'V', 1._wp ) 
     1345      CALL lbc_lnk( e3f , 'F', 1._wp ) 
     1346      CALL lbc_lnk( e3w , 'W', 1._wp ) 
     1347      CALL lbc_lnk( e3uw, 'U', 1._wp ) 
     1348      CALL lbc_lnk( e3vw, 'V', 1._wp ) 
     1349 
     1350      fsdepw(:,:,:) = gdepw (:,:,:) 
     1351      fsde3w(:,:,:) = gdep3w(:,:,:) 
    14891352      ! 
    14901353      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0 
    … …  
    15441407            &                          ' w ', MAXVAL( fse3w (:,:,:) ) 
    15451408      ENDIF 
    1546       ! 
     1409      !  END DO 
    15471410      IF(lwp) THEN                                  ! selected vertical profiles 
    15481411         WRITE(numout,*) 
    … …  
    15741437      ENDIF 
    15751438 
    1576 !!gm bug?  no more necessary?  if ! defined key_helsinki 
    1577       DO jk = 1, jpk 
     1439!================================================================================ 
     1440! check the coordinate makes sense 
     1441!================================================================================ 
     1442      DO ji = 1, jpi 
    15781443         DO jj = 1, jpj 
    1579             DO ji = 1, jpi 
    1580                IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    1581                   WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1582                   CALL ctl_stop( ctmp1 ) 
    1583                ENDIF 
    1584                IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    1585                   WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    1586                   CALL ctl_stop( ctmp1 ) 
    1587                ENDIF 
    1588             END DO 
    1589          END DO 
    1590       END DO 
    1591 !!gm bug    #endif 
     1444 
     1445            IF( hbatt(ji,jj) > 0._wp) THEN 
     1446               DO jk = 1, mbathy(ji,jj) 
     1447                 ! check coordinate is monotonically increasing 
     1448                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1449                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1450                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1451                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
     1452                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     1453                    CALL ctl_stop( ctmp1 ) 
     1454                 ENDIF 
     1455                 ! and check it has never gone negative 
     1456                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
     1457                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1458                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1459                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1460                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1461                    CALL ctl_stop( ctmp1 ) 
     1462                 ENDIF 
     1463                 ! and check it never exceeds the total depth 
     1464                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1465                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1466                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1467                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1468                    CALL ctl_stop( ctmp1 ) 
     1469                 ENDIF 
     1470               END DO 
     1471 
     1472               DO jk = 1, mbathy(ji,jj)-1 
     1473                 ! and check it never exceeds the total depth 
     1474                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1475                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1476                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1477                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1478                    CALL ctl_stop( ctmp1 ) 
     1479                 ENDIF 
     1480               END DO 
     1481 
     1482            ENDIF 
     1483 
     1484         END DO 
     1485      END DO 
    15921486      ! 
    15931487      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1488      ! 
     1489      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
     1490      ! 
     1491   END SUBROUTINE zgr_sco 
     1492 
     1493!!====================================================================== 
     1494   SUBROUTINE s_sh94() 
     1495 
     1496      !!---------------------------------------------------------------------- 
     1497      !!                  ***  ROUTINE s_sh94  *** 
     1498      !!                      
     1499      !! ** Purpose :   stretch the s-coordinate system 
     1500      !! 
     1501      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994 
     1502      !!                mixed S/sigma coordinate 
     1503      !! 
     1504      !! Reference : Song and Haidvogel 1994.  
     1505      !!---------------------------------------------------------------------- 
     1506      ! 
     1507      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1508      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1509      ! 
     1510      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
     1511      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
     1512 
     1513      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
     1514      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1515 
     1516      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
     1517      esigt3  = 0._wp   ;   esigw3  = 0._wp  
     1518      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
     1519      esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1520 
     1521      DO ji = 1, jpi 
     1522         DO jj = 1, jpj 
     1523 
     1524            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     1525               DO jk = 1, jpk 
     1526                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
     1527                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
     1528               END DO 
     1529            ELSE ! shallow water, uniform sigma 
     1530               DO jk = 1, jpk 
     1531                  gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
     1532                  gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
     1533                  END DO 
     1534            ENDIF 
     1535            ! 
     1536            DO jk = 1, jpkm1 
     1537               esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
     1538               esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
     1539            END DO 
     1540            esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
     1541            esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
     1542            ! 
     1543            ! Coefficients for vertical depth as the sum of e3w scale factors 
     1544            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
     1545            DO jk = 2, jpk 
     1546               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1547            END DO 
     1548            ! 
     1549            DO jk = 1, jpk 
     1550               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1551               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1552               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
     1553               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
     1554               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
     1555            END DO 
     1556           ! 
     1557         END DO   ! for all jj's 
     1558      END DO    ! for all ji's 
     1559 
     1560      DO ji = 1, jpim1 
     1561         DO jj = 1, jpjm1 
     1562            DO jk = 1, jpk 
     1563               esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
     1564                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1565               esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
     1566                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1567               esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
     1568                  &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
     1569                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1570               esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
     1571                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1572               esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
     1573                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1574               ! 
     1575               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1576               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1577               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1578               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1579               ! 
     1580               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1581               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1582               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1583            END DO 
     1584        END DO 
     1585      END DO 
     1586 
    15941587      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    15951588      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    1596       ! 
    1597       IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
    1598       ! 
    1599    END SUBROUTINE zgr_sco 
     1589 
     1590   END SUBROUTINE s_sh94 
     1591 
     1592   SUBROUTINE s_sf12 
     1593 
     1594      !!---------------------------------------------------------------------- 
     1595      !!                  ***  ROUTINE s_sf12 ***  
     1596      !!                      
     1597      !! ** Purpose :   stretch the s-coordinate system 
     1598      !! 
     1599      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     1600      !!                mixed S/sigma/Z coordinate 
     1601      !! 
     1602      !!                This method allows the maintenance of fixed surface and or 
     1603      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1604      !!                within an analytically derived stretched S-coordinate framework. 
     1605      !! 
     1606      !! 
     1607      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
     1608      !!---------------------------------------------------------------------- 
     1609      ! 
     1610      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1611      REAL(wp) ::   fsmth          ! smoothing around critical depth 
     1612      REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space 
     1613      ! 
     1614      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
     1615      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
     1616 
     1617      ! 
     1618      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
     1619      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1620 
     1621      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
     1622      esigt3  = 0._wp   ;   esigw3  = 0._wp  
     1623      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
     1624      esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1625 
     1626      DO ji = 1, jpi 
     1627         DO jj = 1, jpj 
     1628 
     1629          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
     1630               
     1631              zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
     1632                                                     ! could be changed by users but care must be taken to do so carefully 
     1633              zbb = 1.0_wp-(zbb/hbatt(ji,jj)) 
     1634             
     1635              zss = rn_zs / hbatt(ji,jj)  
     1636               
     1637              IF (rn_efold /= 0.0_wp) THEN 
     1638                fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
     1639              ELSE 
     1640                fsmth = 1.0_wp  
     1641              ENDIF 
     1642                
     1643              DO jk = 1, jpk 
     1644                gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
     1645                gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 
     1646              ENDDO 
     1647              gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  ) 
     1648              gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  ) 
     1649  
     1650          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 
     1651 
     1652            DO jk = 1, jpk 
     1653              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
     1654              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
     1655            END DO 
     1656 
     1657          ELSE  ! shallow water, z coordinates 
     1658 
     1659            DO jk = 1, jpk 
     1660              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     1661              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     1662            END DO 
     1663 
     1664          ENDIF 
     1665 
     1666          DO jk = 1, jpkm1 
     1667             esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
     1668             esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
     1669          END DO 
     1670          esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  )) 
     1671          esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 
     1672 
     1673          ! Coefficients for vertical depth as the sum of e3w scale factors 
     1674          gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 
     1675          DO jk = 2, jpk 
     1676             gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1677          END DO 
     1678 
     1679          DO jk = 1, jpk 
     1680             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk) 
     1681             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk) 
     1682             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk) 
     1683          END DO 
     1684 
     1685        ENDDO   ! for all jj's 
     1686      ENDDO    ! for all ji's 
     1687 
     1688      DO ji=1,jpi 
     1689        DO jj=1,jpj 
     1690 
     1691          DO jk = 1, jpk 
     1692                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 
     1693                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1694                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 
     1695                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1696                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  & 
     1697                                      hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 
     1698                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1699                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 
     1700                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1701                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 
     1702                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1703 
     1704             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk) 
     1705             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk) 
     1706             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk) 
     1707             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk) 
     1708             ! 
     1709             e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk) 
     1710             e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk) 
     1711             e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk) 
     1712          END DO 
     1713 
     1714        ENDDO 
     1715      ENDDO 
     1716 
     1717      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
     1718      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1719 
     1720   END SUBROUTINE s_sf12 
     1721 
     1722   SUBROUTINE s_tanh() 
     1723 
     1724      !!---------------------------------------------------------------------- 
     1725      !!                  ***  ROUTINE s_tanh***  
     1726      !!                      
     1727      !! ** Purpose :   stretch the s-coordinate system 
     1728      !! 
     1729      !! ** Method  :   s-coordinate stretch  
     1730      !! 
     1731      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1732      !!---------------------------------------------------------------------- 
     1733 
     1734      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1735      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1736 
     1737      DO jk = 1, jpk 
     1738        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
     1739        gsigt(jk) = -fssig( REAL(jk,wp)        ) 
     1740      END DO 
     1741      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
     1742      ! 
     1743      ! Coefficients for vertical scale factors at w-, t- levels 
     1744!!gm bug :  define it from analytical function, not like juste bellow.... 
     1745!!gm        or betteroffer the 2 possibilities.... 
     1746      DO jk = 1, jpkm1 
     1747         esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
     1748         esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
     1749      END DO 
     1750      esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
     1751      esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
     1752      ! 
     1753      ! Coefficients for vertical depth as the sum of e3w scale factors 
     1754      gsi3w(1) = 0.5_wp * esigw(1) 
     1755      DO jk = 2, jpk 
     1756         gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
     1757      END DO 
     1758!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
     1759      DO jk = 1, jpk 
     1760         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1761         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1762         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
     1763         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
     1764         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
     1765      END DO 
     1766!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     1767      DO jj = 1, jpj 
     1768         DO ji = 1, jpi 
     1769            DO jk = 1, jpk 
     1770              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1771              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1772              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1773              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
     1774              ! 
     1775              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1776              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1777              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1778            END DO 
     1779         END DO 
     1780      END DO 
     1781   END SUBROUTINE s_tanh 
     1782 
     1783   FUNCTION fssig( pk ) RESULT( pf ) 
     1784      !!---------------------------------------------------------------------- 
     1785      !!                 ***  ROUTINE fssig *** 
     1786      !!        
     1787      !! ** Purpose :   provide the analytical function in s-coordinate 
     1788      !!           
     1789      !! ** Method  :   the function provide the non-dimensional position of 
     1790      !!                T and W (i.e. between 0 and 1) 
     1791      !!                T-points at integer values (between 1 and jpk) 
     1792      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1793      !!---------------------------------------------------------------------- 
     1794      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
     1795      REAL(wp)             ::   pf   ! sigma value 
     1796      !!---------------------------------------------------------------------- 
     1797      ! 
     1798      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
     1799         &     - TANH( rn_thetb * rn_theta                                )  )   & 
     1800         & * (   COSH( rn_theta                           )                      & 
     1801         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
     1802         & / ( 2._wp * SINH( rn_theta ) ) 
     1803      ! 
     1804   END FUNCTION fssig 
     1805 
     1806 
     1807   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
     1808      !!---------------------------------------------------------------------- 
     1809      !!                 ***  ROUTINE fssig1 *** 
     1810      !! 
     1811      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
     1812      !! 
     1813      !! ** Method  :   the function provides the non-dimensional position of 
     1814      !!                T and W (i.e. between 0 and 1) 
     1815      !!                T-points at integer values (between 1 and jpk) 
     1816      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1817      !!---------------------------------------------------------------------- 
     1818      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
     1819      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
     1820      REAL(wp)             ::   pf1   ! sigma value 
     1821      !!---------------------------------------------------------------------- 
     1822      ! 
     1823      IF ( rn_theta == 0 ) then      ! uniform sigma 
     1824         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
     1825      ELSE                        ! stretched sigma 
     1826         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
     1827            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
     1828            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
     1829      ENDIF 
     1830      ! 
     1831   END FUNCTION fssig1 
     1832 
     1833 
     1834   FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam ) 
     1835      !!---------------------------------------------------------------------- 
     1836      !!                 ***  ROUTINE fgamma  *** 
     1837      !! 
     1838      !! ** Purpose :   provide analytical function for the s-coordinate 
     1839      !! 
     1840      !! ** Method  :   the function provides the non-dimensional position of 
     1841      !!                T and W (i.e. between 0 and 1) 
     1842      !!                T-points at integer values (between 1 and jpk) 
     1843      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1844      !! 
     1845      !!                This method allows the maintenance of fixed surface and or 
     1846      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1847      !!                within an analytically derived stretched S-coordinate framework. 
     1848      !! 
     1849      !! Reference  :   Siddorn and Furner, in prep 
     1850      !!---------------------------------------------------------------------- 
     1851      REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate 
     1852      REAL(wp)                ::   gam(jpk)   ! stretched coordinate 
     1853      REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth 
     1854      REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth 
     1855      REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter 
     1856      REAL(wp)                ::   a1,a2,a3 ! local variables 
     1857      REAL(wp)                ::   n1,n2    ! local variables 
     1858      REAL(wp)                ::   A,B,X    ! local variables 
     1859      integer                 ::   jk 
     1860      !!---------------------------------------------------------------------- 
     1861      ! 
     1862 
     1863      n1  =  1./(jpk-1.) 
     1864      n2  =  1. -  n1 
     1865 
     1866      a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp)  
     1867      a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp) 
     1868      a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1) 
     1869      
     1870      A = Zbb - a3*(Zss-a1)-a2 
     1871      A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) ) 
     1872      B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1) 
     1873      X = 1.0_wp-A/2.0_wp-B 
     1874  
     1875      DO jk = 1, jpk 
     1876        gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 
     1877        gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F) 
     1878      ENDDO  
     1879 
     1880      ! 
     1881   END FUNCTION fgamma 
    16001882 
    16011883   !!====================================================================== 

  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90

    r3294 r3600  
    1919      jpidta  = 198,        &  !: first horizontal dimension > or = to jpi 
    2020      jpjdta  = 224,        &  !: second                     > or = to jpj 
    21       jpkdta  = 33,         &  !: number of levels           > or = to jpk 
     21      jpkdta  = 51,         &  !: number of levels           > or = to jpk 
    2222      ! total domain matrix size 
    2323      jpiglo  = jpidta,      &  !: first  dimension of global domain --> i 
Note: See TracChangeset for help on using the changeset viewer.