Index: branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2783)
+++ branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2840)
@@ -3,5 +3,5 @@
    !!                       ***  MODULE  eosbn2  ***
    !! Ocean diagnostic variable : equation of state - in situ and potential density
-   !!                                               - Brunt-Vaisala frequency 
+   !!                                               - Brunt-Vaisala frequency
    !!==============================================================================
    !! History :  OPA  ! 1989-03  (O. Marti)  Original code
@@ -27,5 +27,5 @@
    !!   eos_insitu_2d  : Compute the in situ density for 2d fields
    !!   eos_bn2        : Compute the Brunt-Vaisala frequency
-   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff.
+   !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio
    !!   tfreez         : Compute the surface freezing temperature
    !!   eos_init       : set eos parameters (namelist)
@@ -41,11 +41,11 @@
    PRIVATE
 
-   !                   !! * Interface 
+   !                   !! * Interface
    INTERFACE eos
       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
-   END INTERFACE 
+   END INTERFACE
    INTERFACE bn2
       MODULE PROCEDURE eos_bn2
-   END INTERFACE 
+   END INTERFACE
 
    PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules
@@ -61,5 +61,5 @@
 
    REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio
-   
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -75,6 +75,6 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE eos_insitu  ***
-      !! 
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
       !!       potential temperature and salinity using an equation of state
       !!       defined through the namelist parameter nn_eos.
@@ -134,5 +134,5 @@
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
-         !  
+         !
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -199,8 +199,8 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
-      !!           
+      !!
       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
       !!      potential volumic mass (Kg/m3) from potential temperature and
-      !!      salinity fields using an equation of state defined through the 
+      !!      salinity fields using an equation of state defined through the
       !!     namelist parameter nn_eos.
       !!
@@ -230,5 +230,5 @@
       !!      nn_eos = 2 : linear equation of state function of temperature and
       !!               salinity
-      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
+      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0
       !!                       = rn_beta * s - rn_alpha * tn - 1.
       !!              rhop(t,s)  = rho(t,s)
@@ -265,5 +265,5 @@
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
-         !  
+         !
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -336,5 +336,5 @@
       !!                  ***  ROUTINE eos_insitu_2d  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
+      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
       !!      potential temperature and salinity using an equation of state
       !!      defined through the namelist parameter nn_eos. * 2D field case
@@ -374,5 +374,5 @@
       !                                                           ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m]
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
       !!
       INTEGER  ::   ji, jj                    ! dummy loop indices
@@ -449,5 +449,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
+               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -468,5 +468,5 @@
       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time-
       !!      step of the input arguments
-      !!      
+      !!
       !! ** Method :
       !!       * nn_eos = 0  : UNESCO sea water properties
@@ -482,5 +482,5 @@
       !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
       !!      The use of potential density to compute N^2 introduces e r r o r
-      !!      in the sign of N^2 at great depths. We recommand the use of 
+      !!      in the sign of N^2 at great depths. We recommand the use of
       !!      nn_eos = 0, except for academical studies.
       !!        Macro-tasked on horizontal slab (jk-loop)
@@ -497,5 +497,5 @@
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars 
+      REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars
 #if defined key_zdfddm
       REAL(wp) ::   zds   ! local scalars
@@ -504,5 +504,5 @@
 
       ! pn2 : interior points only (2=< jk =< jpkm1 )
-      ! -------------------------- 
+      ! --------------------------
       !
       SELECT CASE( nn_eos )
@@ -542,5 +542,5 @@
                      &                                - 0.121555e-07_wp ) * zh
                      !
-                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
+                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2
                      &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
                      &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
@@ -565,10 +565,10 @@
                &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   &
                &               / fse3w(:,:,jk) * tmask(:,:,jk)
-         END DO 
+         END DO
 #if defined key_zdfddm
          DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s])
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  
+                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
                   IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
                   rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
@@ -587,31 +587,31 @@
 
 
-   SUBROUTINE eos_alpbet( pts, palph, pbeta )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE ldf_slp_grif  ***
-      !!
-      !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points. 
-      !!
-      !! ** Method  :   calculates alpha and beta at T-points 
+   SUBROUTINE eos_alpbet( pts, palpbet, beta0 )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_alpbet  ***
+      !!
+      !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points
+      !!
+      !! ** Method  :   calculates alpha / beta ratio at T-points
       !!       * nn_eos = 0  : UNESCO sea water properties
-      !!         The brunt-vaisala frequency is computed using the polynomial
-      !!      polynomial expression of McDougall (1987):
-      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
-      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
-      !!      computed and used in zdfddm module :
-      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
+      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
+      !!                       polynomial expression of McDougall (1987).
+      !!                       Scalar beta0 is returned = 1.
       !!       * nn_eos = 1  : linear equation of state (temperature only)
-      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
+      !!                       The ratio is undefined, so we return alpha as palpbet
+      !!                       Scalar beta0 is returned = 0.
       !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
-      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
-      !!       * nn_eos = 3  : Jackett JAOT 2003 ???
-      !!
-      !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff.
-      !
+      !!                       The alpha/beta ratio is returned as ralpbet
+      !!                       Scalar beta0 is returned = 1.
+      !!
+      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points
+      !!            :   beta0   : 1. or 0.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio
+      REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T)
+      !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh   ! local scalars 
+      REAL(wp) ::   zt, zs, zh   ! local scalars
       !!----------------------------------------------------------------------
       !
@@ -624,42 +624,31 @@
                   zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
                   zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
-                  zh = fsdept(ji,jj,jk)              ! depth in meters 
-                  !
-                  pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   &
-                     &                                        - 0.301985e-05_wp ) * zt   &
-                     &           + 0.785567e-03_wp                                       &
-                     &           + (     0.515032e-08_wp * zs                            &
-                     &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   &
-                     &           + ( (   0.121551e-17_wp * zh                            &
-                     &                 - 0.602281e-15_wp * zs                            &
-                     &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   &
-                     &                                        + 0.408195e-10_wp   * zs   &
-                     &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   &
-                     &                                        - 0.121555e-07_wp ) * zh
-                     !
-                  palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             &
-                      &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &
-                      &                                  - 0.203814e-03_wp ) * zt   &
-                      &                                  + 0.170907e-01_wp ) * zt   &
-                      &   + 0.665157e-01_wp                                         &
-                      &   +     ( - 0.678662e-05_wp * zs                            &
-                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
-                      &   +   ( ( - 0.302285e-13_wp * zh                            &
-                      &           - 0.251520e-11_wp * zs                            &
-                      &           + 0.512857e-12_wp * zt * zt              ) * zh   &
-                      &           - 0.164759e-06_wp * zs                            &
-                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
-                      &                                  + 0.380374e-04_wp ) * zh)
+                  zh = fsdept(ji,jj,jk)               ! depth in meters
+                  !
+                  palpbet(ji,jj,jk) =                                              &
+                     &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &
+                     &                                  - 0.203814e-03_wp ) * zt   &
+                     &                                  + 0.170907e-01_wp ) * zt   &
+                     &   + 0.665157e-01_wp                                         &
+                     &   +     ( - 0.678662e-05_wp * zs                            &
+                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
+                     &   +   ( ( - 0.302285e-13_wp * zh                            &
+                     &           - 0.251520e-11_wp * zs                            &
+                     &           + 0.512857e-12_wp * zt * zt              ) * zh   &
+                     &           - 0.164759e-06_wp * zs                            &
+                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
+                     &                                  + 0.380374e-04_wp ) * zh
                END DO
             END DO
          END DO
-         !
-      CASE ( 1 )
-         palph(:,:,:) = - rn_alpha
-         pbeta(:,:,:) =   0._wp
-         !
-      CASE ( 2 )
-         palph(:,:,:) = - rn_alpha
-         pbeta(:,:,:) =   rn_beta
+         beta0 = 1._wp
+         !
+      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
+         palpbet(:,:,:) = rn_alpha
+         beta0 = 0._wp
+         !
+      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
+         palpbet(:,:,:) = ralpbet
+         beta0 = 1._wp
          !
       CASE DEFAULT
@@ -747,3 +736,3 @@
 
    !!======================================================================
-END MODULE eosbn2  
+END MODULE eosbn2
