Index: /branches/2012/dev_UKMO_2012/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- /branches/2012/dev_UKMO_2012/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3617)
+++ /branches/2012/dev_UKMO_2012/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3618)
@@ -759,10 +759,22 @@
 
 Options 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.
-A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is provided as an example:
+
+The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true (\np{ln\_sco\_SH94}~=~false and \np{ln\_sco\_SF12}~=~false.) This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
+
+\begin{equation}
+  z = s_{min}+C\left(s\right)\left(H-s_{min}\right)
+  \label{eq:SH94_1}
+\end{equation}
+
+where $s_{min}$ is the depth at which the s-coordinate stretching starts and allows a z-coordinate to placed on top of the stretched coordinate, and z is the depth (negative down from the asea surface).
+
+\begin{equation}
+  s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1
+  \label{eq:s}
+\end{equation}
 
 \begin{equation} \label{DOM_sco_function}
 \begin{split}
-z	&= h_c +( h-h_c)\;c s	\\
-c(s)	&=  \frac{ \left[	  \tanh{ \left( \theta \, (s+b) \right)} 
+C(s)	&=  \frac{ \left[	  \tanh{ \left( \theta \, (s+b) \right)} 
 	  	 			- \tanh{ \left(  \theta \, b      \right)}  \right]}
 		      {2\;\sinh \left( \theta \right)}
@@ -770,4 +782,11 @@
 \end{equation}
 
+A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling:
+
+\begin{equation}
+  C\left(s\right) =   \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} +      \\
+  b\frac{ \tanh \left[ \theta \left(s + \frac{1}{2} \right)\right] - \tanh\left(\frac{\theta}{2}\right)}{ 2\tanh\left (\frac{\theta}{2}\right)}
+  \label{eq:SH94_2}
+\end{equation}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -780,5 +799,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-where $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 
+where $H_c$ is the critical depth (\np{rn\_hc}) 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 
 bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 
 $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom 
@@ -788,5 +807,5 @@
 
 \begin{equation}
-z = \gamma\left(h+\zeta\right) \quad \text{ with } \quad 0 \leq \gamma \leq 1
+z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1
 \label{eq:z}
 \end{equation}
@@ -800,5 +819,5 @@
 Where:
 \begin{equation} \label{DOM_gamma}
-f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2}
+f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1} 
 \end{equation}
 
@@ -819,5 +838,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-This 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.
+This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell thicknesses 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.
 
 As 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$.
Index: /branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3617)
+++ /branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3618)
@@ -15,5 +15,5 @@
    !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
    !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
-   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn adn Furner stretching functio
+   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
    !!----------------------------------------------------------------------
 
@@ -28,7 +28,7 @@
    !!       zgr_zps      : z-coordinate with partial steps
    !!       zgr_sco      : s-coordinate
-   !!       fssig        : sigma coordinate non-dimensional function
-   !!       fgamma       : Siddorn and Furner stretching function
-   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
+   !!       fssig        : tanh stretch function
+   !!       fssig1       : Song and Haidvogel 1994 stretch function
+   !!       fgamma       : Siddorn and Furner 2012 stretching function
    !!---------------------------------------------------------------------
    USE oce               ! ocean variables
@@ -51,5 +51,4 @@
    LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
    LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
-   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
    !
    REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m)
@@ -63,4 +62,5 @@
    !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
    ! Siddorn and Furner stretching parameters
+   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
    REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
    REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord
@@ -1065,10 +1065,10 @@
       !!            hbatv = mj( hbatt )
       !!            hbatf = mi( mj( hbatt ) )
-      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
+      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
       !!         function and its derivative given as function.
-      !!            gsigt(k) = fssig (k    )
-      !!            gsigw(k) = fssig (k-0.5)
-      !!            esigt(k) = fsdsig(k    )
-      !!            esigw(k) = fsdsig(k-0.5)
+      !!            z_gsigt(k) = fssig (k    )
+      !!            z_gsigw(k) = fssig (k-0.5)
+      !!            z_esigt(k) = fsdsig(k    )
+      !!            z_esigw(k) = fsdsig(k-0.5)
       !!      Three options for stretching are give, and they can be modified
       !!      following the users requirements. Nevertheless, the output as
@@ -1125,4 +1125,5 @@
          WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
          WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
+         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
          WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
          WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
@@ -1383,12 +1384,4 @@
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
 
-      !                                               ! =============
-      IF(lwp) THEN                                    ! Control print
-         !                                            ! =============
-         WRITE(numout,*) 
-         WRITE(numout,*) ' domzgr: vertical coefficients for model level'
-         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
-         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
-      ENDIF
       IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
          WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
@@ -1508,14 +1501,14 @@
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
-
-      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
-      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
-
-      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
-      esigt3  = 0._wp   ;   esigw3  = 0._wp 
-      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
-      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
+      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
+      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
+      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
 
       DO ji = 1, jpi
@@ -1524,25 +1517,25 @@
             IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
                DO jk = 1, jpk
-                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
-                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
+                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
+                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
                END DO
             ELSE ! shallow water, uniform sigma
                DO jk = 1, jpk
-                  gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
-                  gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
+                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
+                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
                   END DO
             ENDIF
             !
             DO jk = 1, jpkm1
-               esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
-               esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
-            END DO
-            esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
-            esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
+               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
+               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
+            END DO
+            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
+            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
             !
             ! Coefficients for vertical depth as the sum of e3w scale factors
-            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
+            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
             DO jk = 2, jpk
-               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
+               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
             END DO
             !
@@ -1550,7 +1543,7 @@
                zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
                zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
-               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
+               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
             END DO
            !
@@ -1561,30 +1554,30 @@
          DO jj = 1, jpjm1
             DO jk = 1, jpk
-               esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
+               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
+               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-               esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
-                  &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
+               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
+                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-               esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
+               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
+               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
                !
-               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
                !
-               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
-               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
             END DO
         END DO
       END DO
 
-      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
+      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
 
    END SUBROUTINE s_sh94
@@ -1609,18 +1602,18 @@
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
-      REAL(wp) ::   fsmth          ! smoothing around critical depth
-      REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space
-      !
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
-
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
-      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
-
-      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
-      esigt3  = 0._wp   ;   esigw3  = 0._wp 
-      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
-      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
+      REAL(wp) ::   zsmth               ! smoothing around critical depth
+      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
+      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
+      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
+      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
 
       DO ji = 1, jpi
@@ -1629,28 +1622,28 @@
           IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
               
-              zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
+              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
                                                      ! could be changed by users but care must be taken to do so carefully
-              zbb = 1.0_wp-(zbb/hbatt(ji,jj))
+              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
             
-              zss = rn_zs / hbatt(ji,jj) 
+              zzs = rn_zs / hbatt(ji,jj) 
               
               IF (rn_efold /= 0.0_wp) THEN
-                fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
+                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
               ELSE
-                fsmth = 1.0_wp 
+                zsmth = 1.0_wp 
               ENDIF
                
               DO jk = 1, jpk
-                gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
-                gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
+                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
               ENDDO
-              gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  )
-              gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  )
+              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
+              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
  
           ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
 
             DO jk = 1, jpk
-              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
-              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
             END DO
 
@@ -1658,6 +1651,6 @@
 
             DO jk = 1, jpk
-              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
-              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
             END DO
 
@@ -1665,20 +1658,20 @@
 
           DO jk = 1, jpkm1
-             esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
-             esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
+             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
+             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
           END DO
-          esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ))
-          esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
+          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
+          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
 
           ! Coefficients for vertical depth as the sum of e3w scale factors
-          gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
+          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
           DO jk = 2, jpk
-             gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
+             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
           END DO
 
           DO jk = 1, jpk
-             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk)
-             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk)
-             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk)
+             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
           END DO
 
@@ -1690,31 +1683,32 @@
 
           DO jk = 1, jpk
-                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &
+                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &
+                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  &
-                                      hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &
+                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
+                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &
+                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &
+                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
 
-             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk)
-             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk)
-             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk)
-             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk)
+             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
+             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
+             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
+             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
              !
-             e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk)
-             e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk)
-             e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk)
+             e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
+             e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
+             e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
           END DO
 
         ENDDO
       ENDDO
-
-      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
-      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
+      !                                               ! =============
+
+      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
+      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
 
    END SUBROUTINE s_sf12
@@ -1735,9 +1729,18 @@
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
 
+      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
+      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
+
+      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
+      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               )
+
+      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
+      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
+
       DO jk = 1, jpk
-        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
-        gsigt(jk) = -fssig( REAL(jk,wp)        )
-      END DO
-      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
+        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
+        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
+      END DO
+      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
       !
       ! Coefficients for vertical scale factors at w-, t- levels
@@ -1745,14 +1748,14 @@
 !!gm        or betteroffer the 2 possibilities....
       DO jk = 1, jpkm1
-         esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
-         esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
-      END DO
-      esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
-      esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
+         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
+         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
+      END DO
+      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
+      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
       !
       ! Coefficients for vertical depth as the sum of e3w scale factors
-      gsi3w(1) = 0.5_wp * esigw(1)
+      z_gsi3w(1) = 0.5_wp * z_esigw(1)
       DO jk = 2, jpk
-         gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
+         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
       END DO
 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
@@ -1760,7 +1763,7 @@
          zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
          zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
-         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
-         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
+         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
       END DO
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
@@ -1768,15 +1771,19 @@
          DO ji = 1, jpi
             DO jk = 1, jpk
-              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
-              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
+              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
               !
-              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
-            END DO
-         END DO
-      END DO
+              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+            END DO
+         END DO
+      END DO
+
+      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      )
+      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               )
+
    END SUBROUTINE s_tanh
 
@@ -1832,5 +1839,5 @@
 
 
-   FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam )
+   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE fgamma  ***
@@ -1849,31 +1856,31 @@
       !! Reference  :   Siddorn and Furner, in prep
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate
-      REAL(wp)                ::   gam(jpk)   ! stretched coordinate
-      REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth
-      REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth
-      REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter
-      REAL(wp)                ::   a1,a2,a3 ! local variables
-      REAL(wp)                ::   n1,n2    ! local variables
-      REAL(wp)                ::   A,B,X    ! local variables
+      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
+      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
+      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
+      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
+      REAL(wp)                ::   za1,za2,za3    ! local variables
+      REAL(wp)                ::   zn1,zn2        ! local variables
+      REAL(wp)                ::   za,zb,zx       ! local variables
       integer                 ::   jk
       !!----------------------------------------------------------------------
       !
 
-      n1  =  1./(jpk-1.)
-      n2  =  1. -  n1
-
-      a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp) 
-      a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp)
-      a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1)
+      zn1  =  1./(jpk-1.)
+      zn2  =  1. -  zn1
+
+      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
+      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
+      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
      
-      A = Zbb - a3*(Zss-a1)-a2
-      A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) )
-      B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1)
-      X = 1.0_wp-A/2.0_wp-B
+      za = pzb - za3*(pzs-za1)-za2
+      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
+      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
+      zx = 1.0_wp-za/2.0_wp-zb
  
       DO jk = 1, jpk
-        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) )
-        gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F)
+        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
+        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
       ENDDO 
 
