Index: trunk/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1577)
+++ trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1601)
@@ -42,10 +42,14 @@
 
 !!gm   DOCTOR name for the namelist parameter...
-   !                                 !!! ** Namelist nam_zgr_sco **
-   REAL(wp) ::   sbot_min =  300.     ! minimum depth of s-bottom surface (>0) (m)
-   REAL(wp) ::   sbot_max = 5250.     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   REAL(wp) ::   theta    =    6.0    ! surface control parameter (0<=theta<=20)
-   REAL(wp) ::   thetb    =    0.75   ! bottom control parameter  (0<=thetb<= 1)
-   REAL(wp) ::   r_max    =    0.15   ! maximum cut-off r-value allowed (0<r_max<1)
+   !                                    !!! ** Namelist namzgr_sco **
+   REAL(wp) ::   rn_sbot_min =  300.     ! minimum depth of s-bottom surface (>0) (m)
+   REAL(wp) ::   rn_sbot_max = 5250.     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   REAL(wp) ::   rn_theta    =    6.0    ! surface control parameter (0<=rn_theta<=20)
+   REAL(wp) ::   rn_thetb    =    0.75   ! bottom control parameter  (0<=rn_thetb<= 1)
+   REAL(wp) ::   rn_rmax     =    0.15   ! maximum cut-off r-value allowed (0<rn_rmax<1)
+   LOGICAL  ::   ln_s_sigma  = .false.   ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
+   REAL(wp) ::   rn_bb       =    0.8    ! stretching parameter for song and haidvogel stretching
+   !                                     ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
+   REAL(wp) ::   rn_hc       = 150.      ! Critical depth for s-sigma coordinates
  
    !! * Substitutions
@@ -79,9 +83,9 @@
       INTEGER ::   ioptio = 0   ! temporary integer
       !!
-      NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
-      !!----------------------------------------------------------------------
-
-      REWIND ( numnam )                ! Read Namelist nam_zgr : vertical coordinate'
-      READ   ( numnam, nam_zgr )
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+      !!----------------------------------------------------------------------
+
+      REWIND ( numnam )                ! Read Namelist namzgr : vertical coordinate'
+      READ   ( numnam, namzgr )
 
       IF(lwp) THEN                     ! Control print
@@ -89,5 +93,5 @@
          WRITE(numout,*) 'dom_zgr : vertical coordinate'
          WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '          Namelist nam_zgr : set vertical coordinate'
+         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
          WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
          WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
@@ -232,8 +236,10 @@
       ENDIF
 
+!!gm BUG in s-coordinate this does not work!
       ! deepest/shallowest W level Above/Bellow ~10m
       zrefdep = 10. - ( 0.1*MINVAL(e3w_0) )                          ! ref. depth with tolerance (10% of minimum layer thickness)
       nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Bellow ~10m
       nla10 = nlb10 - 1                                              ! deepest    W level Above  ~10m
+!!gm end bug
 
       IF(lwp) THEN                        ! control print
@@ -1001,14 +1007,14 @@
       !!----------------------------------------------------------------------
       !
-      pf =   (   TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb )  )      &
-         &     - TANH( thetb * theta                                )  )   &
-         & * (   COSH( theta                           )                   &
-         &     + COSH( theta * ( 2.e0 * thetb - 1.e0 ) )  )                &
-         & / ( 2.e0 * SINH( theta ) )
+      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb )  )      &
+         &     - TANH( rn_thetb * rn_theta                                )  )   &
+         & * (   COSH( rn_theta                           )                   &
+         &     + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) )  )                &
+         & / ( 2.e0 * SINH( rn_theta ) )
       !
    END FUNCTION fssig
 
 
-   FUNCTION fssig1( pk1, bb ) RESULT( pf1 )
+   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE eos_init  ***
@@ -1024,14 +1030,14 @@
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ) ::   pk1   ! continuous "k" coordinate
-      REAL(wp), INTENT(in   ) ::   bb    ! Stretching coefficient
+      REAL(wp), INTENT(in   ) ::   pbb   ! Stretching coefficient
       REAL(wp)                ::   pf1   ! sigma value
       !!----------------------------------------------------------------------
       !
-      IF ( theta == 0 ) then      ! uniform sigma
+      IF ( rn_theta == 0 ) then      ! uniform sigma
          pf1 = -(pk1-0.5) / REAL( jpkm1 )
       ELSE                        ! stretched sigma
-         pf1 =   (1.0-bb) * (sinh( theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(theta) + &
-            &    bb * ( (tanh( theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*theta) ) / &
-            &    (2*tanh(0.5*theta) ) )
+         pf1 =   (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + &
+            &    pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) / &
+            &    (2*tanh(0.5*rn_theta) ) )
       ENDIF
       !
@@ -1077,7 +1083,4 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
       !!
-      LOGICAL  :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true.
-      REAL(wp) :: bb = 0.8   ! stretching parameter for song and haidvogel stretching, bb=0; top only, bb =1; top and bottom
-      REAL(wp) :: hc = 150   ! Critical depth for s-sigma coordinates
 !!gm never do that !!!!   ==> Pb at compilation phase on several computer
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigw3 = 0.0d0
@@ -1093,9 +1096,9 @@
 !!gm end
       !!
-      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max, ln_s_sigma, bb, hc
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam )                        ! Read Namelist nam_zgr_sco : sigma-stretching parameters
-      READ  ( numnam, nam_zgr_sco )
+      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam )                        ! Read Namelist namzgr_sco : sigma-stretching parameters
+      READ  ( numnam, namzgr_sco )
 
       IF(lwp) THEN                            ! control print
@@ -1103,27 +1106,28 @@
          WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '         Namelist nam_zgr_sco'
-         WRITE(numout,*) '            sigma-stretching coeffs '
-         WRITE(numout,*) '            maximum depth of s-bottom surface (>0)    sbot_max   = ' ,sbot_max
-         WRITE(numout,*) '            minimum depth of s-bottom surface (>0)    sbot_min   = ' ,sbot_min
-         WRITE(numout,*) '            surface control parameter (0<=theta<=20)  theta      = ', theta
-         WRITE(numout,*) '            bottom control parameter  (0<=thetb<= 1)  thetb      = ', thetb
-         WRITE(numout,*) '            maximum cut-off r-value allowed           r_max      = ' , r_max
-         WRITE(numout,*) '            Critical depth                            hc         = ', hc
-         WRITE(numout,*) '            Hybrid s-sigma-coordinate                 ln_s_sigma = ', ln_s_sigma
-      ENDIF
-
-      hift(:,:) = sbot_min                     ! set the minimum depth for the s-coordinate
-      hifu(:,:) = sbot_min
-      hifv(:,:) = sbot_min
-      hiff(:,:) = sbot_min
+         WRITE(numout,*) '   Namelist namzgr_sco'
+         WRITE(numout,*) '      sigma-stretching coeffs '
+         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
+         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
+         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
+         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
+         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
+         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
+         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
+         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
+      ENDIF
+
+      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
+      hifu(:,:) = rn_sbot_min
+      hifv(:,:) = rn_sbot_min
+      hiff(:,:) = rn_sbot_min
 
       !                                        ! set maximum ocean depth
-      bathy(:,:) = MIN( sbot_max, bathy(:,:) )
+      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
       DO jj = 1, jpj
          DO ji = 1, jpi
            IF (bathy(ji,jj) .gt. 0.0) THEN
-              bathy(ji,jj) = MAX( sbot_min, bathy(ji,jj) )
+              bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
            ENDIF
          END DO
@@ -1139,5 +1143,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
+            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
          END DO
       END DO
@@ -1145,5 +1149,5 @@
       zrmax = 1.e0
       !                                                     ! ================ !
-      DO WHILE ( jl <= 10000 .AND. zrmax > r_max )          !  Iterative loop  !
+      DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax )          !  Iterative loop  !
          !                                                  ! ================ !
          jl = jl + 1
@@ -1157,8 +1161,8 @@
                zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
                zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
-               IF( zri(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
-               IF( zri(ji,jj) > r_max )   zmsk(iip1,jj  ) = 1.0
-               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
-               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,ijp1) = 1.0
+               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0
+               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1.0
+               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1.0
+               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1.0
             END DO
          END DO
@@ -1218,5 +1222,5 @@
             DO ji = 1, jpi
                ztaper = EXP( -(gphit(ji,jj)/8)**2 )
-               hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
+               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
             END DO
          END DO
@@ -1239,9 +1243,9 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
-      ENDIF
-      hbatu(:,:) = sbot_min
-      hbatv(:,:) = sbot_min
-      hbatf(:,:) = sbot_min
+         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+      ENDIF
+      hbatu(:,:) = rn_sbot_min
+      hbatv(:,:) = rn_sbot_min
+      hbatf(:,:) = rn_sbot_min
       DO jj = 1, jpjm1
         DO ji = 1, jpim1
@@ -1259,5 +1263,5 @@
          DO ji = 1, jpi
             IF( hbatu(ji,jj) == 0.e0 ) THEN
-               IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = sbot_min
+               IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0.e0 )   hbatu(ji,jj) = zhbat(ji,jj)
             ENDIF
@@ -1268,5 +1272,5 @@
          DO ji = 1, jpi
             IF( hbatv(ji,jj) == 0.e0 ) THEN
-               IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = sbot_min
+               IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0.e0 )   hbatv(ji,jj) = zhbat(ji,jj)
             ENDIF
@@ -1277,5 +1281,5 @@
          DO ji = 1, jpi
             IF( hbatf(ji,jj) == 0.e0 ) THEN
-               IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = sbot_min
+               IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0.e0 )   hbatf(ji,jj) = zhbat(ji,jj)
             ENDIF
@@ -1307,13 +1311,13 @@
       ! non-dimensional "sigma" for model level depth at w- and t-levels
 
-      IF ( ln_s_sigma ) THEN  !Song and Haidvogel style stretched sigma for depths below hc, with uniform sigma in shallower waters
-
-         DO ji=1,jpi
-           DO jj=1,jpj
-
-             IF (hbatt(ji,jj).GT.hc) THEN !deep water, stretched sigma
+      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
+         !                         ! below rn_hc, with uniform sigma in shallower waters
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+
+             IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
                DO jk = 1, jpk
-                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, bb )
-                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , bb )
+                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
+                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
                END DO
              ELSE ! shallow water, uniform sigma
@@ -1342,7 +1346,7 @@
                 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
                 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
-                gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigt3(ji,jj,jk)+hc*zcoeft)
-                gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigw3(ji,jj,jk)+hc*zcoefw)
-                gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsi3w3(ji,jj,jk)+hc*zcoefw)
+                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*zcoefw)
              END DO
 
@@ -1367,12 +1371,12 @@
                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
 
-                e3t(ji,jj,jk)=((hbatt(ji,jj)-hc)*esigt3(ji,jj,jk) + hc/FLOAT(jpkm1))
-                e3u(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigtu3(ji,jj,jk) + hc/FLOAT(jpkm1))
-                e3v(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigtv3(ji,jj,jk) + hc/FLOAT(jpkm1))
-                e3f(ji,jj,jk)=((hbatf(ji,jj)-hc)*esigtf3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
+                e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
+                e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
+                e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
                 !
-                e3w (ji,jj,jk)=((hbatt(ji,jj)-hc)*esigw3(ji,jj,jk) + hc/FLOAT(jpkm1))
-                e3uw(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigwu3(ji,jj,jk) + hc/FLOAT(jpkm1))
-                e3vw(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigwv3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
+                e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
+                e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
              END DO
 
