Index: branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 4944)
+++ branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 4946)
@@ -54,6 +54,12 @@
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-   INTEGER(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
+#ifdef key_agrif
+   ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals
+   REAL(wp),    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
                                                                                           !: (first wet level and last level include in the tbl)
+#else
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
+#endif
+
 
    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
@@ -303,8 +309,11 @@
       sbc_isf_alloc = 0       ! set to zero if no array to be allocated
       IF( .NOT. ALLOCATED( qisf ) ) THEN
-         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj), fwfisf(jpi,jpj), &
-               &    fwfisf_b(jpi,jpj), misfkt(jpi,jpj), rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj),     &
-               &    rzisf_tbl(jpi,jpj), misfkb(jpi,jpj), ttbl(jpi,jpj), stbl(jpi,jpj), utbl(jpi,jpj), &
-               &    vtbl(jpi, jpj), risfLeff(jpi,jpj), rhisf_tbl_0(jpi,jpj), ralpha(jpi,jpj), STAT= sbc_isf_alloc )
+         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts)              , &
+               &    qisf(jpi,jpj)     , fwfisf(jpi,jpj)     , fwfisf_b(jpi,jpj)   , &
+               &    rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
+               &    ttbl(jpi,jpj)     , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
+               &    vtbl(jpi, jpj)    , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
+               &    ralpha(jpi,jpj)   , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
+               &    STAT= sbc_isf_alloc )
          !
          IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
@@ -363,5 +372,5 @@
              ! Calculate freezing temperature
                 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                zt_frz = tfreez1D(tsb(ji,jj,ik,jp_sal), zpress) 
+                zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress) 
                 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
              ENDDO
@@ -445,5 +454,5 @@
       zti(:,:)=tinsitu( ttbl, stbl, zpress )
 ! Calculate freezing temperature
-      zfrz(:,:)=tfreez( sss_m(:,:), zpress )
+      zfrz(:,:)=eos_fzp( sss_m(:,:), zpress )
 
       
@@ -526,5 +535,6 @@
                      nit = nit + 1
                      IF (nit .GE. 51) THEN
-                        WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
+                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
+                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
                         CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
                      END IF
@@ -584,10 +594,10 @@
       REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence 
       REAL(wp) :: zcoef                      ! temporary coef
-      REAL(wp) :: zrhos, zalbet, zbeta, zthermal, zhalin
-      REAL(wp) :: zt, zs, zh
+      REAL(wp) :: zdep
       REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
       REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
       REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
+      REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
       !
@@ -656,44 +666,12 @@
 
       !! compute bouyancy 
-               IF( nn_eos < 1) THEN
-                  zt     = ttbl(ji,jj)
-                  zs     = stbl(ji,jj) - 35.0
-                  zh     = fsdepw(ji,jj,ikt)
-                  !  potential volumic mass
-                  zrhos  = rhop(ji,jj,ikt)
-                  zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
-                     &                               - 0.203814e-03 ) * zt   &
-                     &                               + 0.170907e-01 ) * zt   &
-                     &   + 0.665157e-01                                      &
-                     &   +     ( - 0.678662e-05 * zs                         &
-                     &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
-                     &   +   ( ( - 0.302285e-13 * zh                         &
-                     &           - 0.251520e-11 * zs                         &
-                     &           + 0.512857e-12 * zt * zt           ) * zh   &
-                     &           - 0.164759e-06 * zs                         &
-                     &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
-                     &                               + 0.380374e-04 ) * zh
-
-                  zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
-                     &                            - 0.301985e-05 ) * zt      &
-                     &   + 0.785567e-03                                      &
-                     &   + (     0.515032e-08 * zs                           &
-                     &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
-                     &   +(  (   0.121551e-17 * zh                           &
-                     &         - 0.602281e-15 * zs                           &
-                     &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
-                     &                             + 0.408195e-10   * zs     &
-                     &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
-                     &                             - 0.121555e-07 ) * zh
-
-                  zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
-                  zhalin   = zbeta * stbl(ji,jj) * rcs
-               ELSE
-                  zrhos    = rhop(ji,jj,ikt) + rau0 * ( 1. - tmask(ji,jj,ikt) )
-                  zthermal = rn_alpha / ( rcp * zrhos + epsln )
-                  zhalin   = rn_beta * stbl(ji,jj) * rcs
-               ENDIF
+               zts(jp_tem) = ttbl(ji,jj)
+               zts(jp_sal) = stbl(ji,jj)
+               zdep        = fsdepw(ji,jj,ikt)
+               !
+               CALL eos_rab( zts, zdep, zab )
+                  !
       !! compute length scale 
-               zbuofdep = grav * ( zthermal * zqhisf - zhalin * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !! compute Monin Obukov Length
@@ -766,7 +744,10 @@
                ! level partially include in ice shelf boundary layer 
                zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-               IF (cptin == 'T') varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-               IF (cptin == 'U') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
-               IF (cptin == 'V') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
+               IF (cptin == 'T') &
+                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+               IF (cptin == 'U') &
+                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
+               IF (cptin == 'V') &
+                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
             END IF
          END DO
@@ -796,7 +777,7 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
       !!
-      INTEGER(wp)  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER(wp)  ::   ikt, ikb 
-      INTEGER(wp)  ::   nk_isf
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikt, ikb 
+      INTEGER  ::   nk_isf
       REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
       REAL(wp)     ::   zfact     ! local scalar
@@ -834,8 +815,10 @@
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
+                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
                END DO
                ! level partially include in ice shelf boundary layer 
-               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
+               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
+                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
             !==   ice shelf melting mass distributed over several levels   ==!
          END DO
