Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5904)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5905)
@@ -440,13 +440,12 @@
 !              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  !
 ! nn_isf == 4
-   sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
-   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_fwfisf   = 'rnfisf'  ,         -12       ,'sowflisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 3
-   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_rnfisf   = 'rnfisf'  ,         -12       ,'sofwfisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2 and 3
-   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
-   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_depmax_isf = 'rnfisf' ,        -12       ,'sozisfmax' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_depmin_isf = 'rnfisf' ,        -12       ,'sozisfmin' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''
 ! nn_isf == 2
-   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   ''
+   sn_Leff_isf = 'rnfisf'  ,         -12       ,'Leff'    ,     .false.    , .true.  , 'yearly'  ,  ''      ,   ''
 ! for all case
    nn_isf      = 1         !  ice shelf melting/freezing
@@ -458,12 +457,12 @@
    rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula
 ! only for nn_isf = 1 or 4
-   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008)
+   rn_hisf_tbl =  30.      ! thickness of the top boundary layer    (Losh et al. 2008)
                           ! 0 => thickness of the tbl = thickness of the first wet cell
 ! only for nn_isf = 1
-   nn_isfblk   = 1        ! 1 ISOMIP ; 2 : 3 equation formulation, Jenkins et al. 1991
+   nn_isfblk   = 1        ! 1 ISOMIP  like: 2 equations formulation (Hunter et al., 2006)
+                          ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015)
    nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s)
                           ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
-                          !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
-                          ! 2 = velocity and stability dependent Gamma    Holland et al. 1999
+                          ! 2 = velocity and stability dependent Gamma    (Holland et al. 1999)
 /
 !-----------------------------------------------------------------------
Index: /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5904)
+++ /branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5905)
@@ -68,6 +68,6 @@
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
    CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
-   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
+   TYPE(FLD_N)       , PUBLIC ::   sn_fwfisf            !: information about the isf file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_fwfisf
    TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
@@ -85,17 +85,20 @@
  
   SUBROUTINE sbc_isf(kt)
-    INTEGER, INTENT(in) :: kt         ! ocean time step
-    INTEGER             :: ji, jj, jk, inum, ierror
-    INTEGER             :: ikt, ikb   ! top and bottom level of the isf boundary layer
+    INTEGER, INTENT(in) :: kt                   ! ocean time step
+    INTEGER             :: ji, jj, jk           ! loop index
+    INTEGER             :: ikt, ikb             ! top and bottom level of the isf boundary layer
+    INTEGER             :: inum, ierror
+    INTEGER             :: ios                  ! Local integer output status for namelist read
     REAL(wp)            :: zhk
-    REAL(wp)            ::   zt_frz, zpress
+    REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
     CHARACTER(len=256)  :: cvarzisf, cvarhisf   ! name for isf file
-    CHARACTER(LEN=32 )  :: cvarLeff                    ! variable name for efficient Length scale
-    INTEGER             :: ios           ! Local integer output status for namelist read
+    CHARACTER(LEN=32 )  :: cvarLeff             ! variable name for efficient Length scale
       !
       !!---------------------------------------------------------------------
       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf      , &
-                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
-      !
+                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !
+      ! allocation
+      CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
       !
       !                                         ! ====================== !
@@ -148,5 +151,5 @@
                CALL iom_close(inum)
                !
-               risfLeff = risfLeff*1000           !: convertion in m
+               risfLeff = risfLeff*1000.0_wp           !: convertion in m
             END IF
 
@@ -179,9 +182,7 @@
             
             ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
+            ALLOCATE( sf_fwfisf(1), STAT=ierror )
             ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
-            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
             CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
          END IF
          
@@ -197,6 +198,5 @@
 
                ! determine the deepest level influenced by the boundary layer
-               ! test on tmask useless ?????
-               DO jk = ikt, mbkt(ji,jj)
+               DO jk = ikt+1, mbkt(ji,jj)
                   IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
@@ -247,6 +247,6 @@
             ! specified runoff in depth (Mathiot et al., XXXX in preparation)
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
-            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting) 
+            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat flux
             stbl(:,:)   = soce
 
@@ -254,18 +254,21 @@
             ! specified fwf and heat flux forcing beneath the ice shelf
             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
-            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
-            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)           ! fwf
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
-            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)            ! heat flux
+            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting)
+            qisf(:,:)   = fwfisf(:,:) * lfusisf                ! heat flux
             stbl(:,:)   = soce
-
          END IF
+
          ! compute tsc due to isf
-         ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable).
-!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
-         zt_frz = -1.9_wp !CALL eos_fzp( tsn(ji,jj,jk,jp_sal), zt_frz, zpress )
-         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz * r1_rau0 !
+         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
+         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
+         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj))
+            END DO
+         END DO
+         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
          
-         ! salt effect already take into account in vertical advection
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
          risf_tsc(:,:,jp_sal) = 0.0_wp
 
@@ -276,5 +279,5 @@
          CALL lbc_lnk(qisf(:,:)     ,'T',1.)
 
-         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
             IF( ln_rstart .AND.    &                     ! Restart: read in restart file
                  & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
@@ -293,4 +296,7 @@
          IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf)
       END IF
+
+      ! deallocation
+      CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
   
   END SUBROUTINE sbc_isf
@@ -360,6 +366,5 @@
              ! after verif with UNESCO, wrong sign in BG eq. 2
              ! Calculate freezing temperature
-                zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
-                CALL eos_fzp(tsb(ji,jj,jk,jp_sal), zt_frz, zpress) 
+                CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
                 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
              ENDDO
@@ -398,5 +403,5 @@
       INTEGER, INTENT(in)          ::   kt         ! ocean time step
       !
-      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz,zpress,zti
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
       REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
       REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
@@ -411,35 +416,26 @@
       !!---------------------------------------------------------------------
       !
-      ! coeficient for linearisation of tfreez
-      zlamb1=-0.0575_wp
-      zlamb2= 0.0901_wp
-      zlamb3=-7.61e-04
+      ! coeficient for linearisation of potential tfreez
+      ! Crude approximation for pressure (but commonly used)
+      zlamb1=-0.0573_wp
+      zlamb2= 0.0832_wp
+      zlamb3=-7.53e-08 * grav * rau0
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
-      CALL wrk_alloc( jpi,jpj, zfrz  , zpress, zti     , zgammat, zgammas )
-      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b                   )
+      CALL wrk_alloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
 
       ! initialisation
-      zpress (:,:)=0.0_wp
       zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0;
       zhtflx (:,:)=0.0_wp     ; zhtflx_b(:,:)=0.0_wp    ;
       zfwflx (:,:)=0.0_wp
-      !
-      ! Crude approximation for pressure (but commonly used)
-      ! 1e-04 to convert from Pa to dBar
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
-         END DO
-      END DO
-
-      ! Calculate in-situ temperature (ref to surface)
-      zti(:,:) = tinsitu( ttbl, stbl, zpress )
-
+
+      ! compute ice shelf melting
       nit = 1 ; lit = .TRUE.
       DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
-         IF (nn_isfblk == 1) THEN ! ISOMIP case
+         IF (nn_isfblk == 1) THEN 
+            ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
             ! Calculate freezing temperature
-            CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress(:,:) )
+            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
 
             ! compute gammat every where (2d)
@@ -449,5 +445,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
+                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
                   zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf
                END DO
@@ -455,14 +451,14 @@
 
             ! Compute heat flux and upward fresh water flux
-            ! For genuine ISOMIP protocol this should probably be something like
-            qisf  (:,:) = - zhtflx(:,:)                                 * (1._wp - tmask(:,:,1)) * ssmask(:,:)
-            fwfisf(:,:) =   zfwflx(:,:) * ( soce / MAX(stbl(:,:),zeps)) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
 
          ELSE IF (nn_isfblk == 2 ) THEN
-            ! More complicated 3 equation thermodynamics as in MITgcm
+            ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
             ! compute gammat every where (2d)
             CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
 
             ! compute upward heat flux zhtflx and upward water flux zwflx
+            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -472,5 +468,5 @@
                   zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
                   zeps4=zlamb2+zlamb3*risfdep(ji,jj)
-                  zeps6=zeps4-zti(ji,jj)
+                  zeps6=zeps4-ttbl(ji,jj)
                   zeps7=zeps4-tsurf
                   zaqe=zlamb1 * (zeps1 + zeps3)
@@ -485,12 +481,12 @@
                   IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
 
-                  ! compute t freeze
+                  ! compute t freeze (eq. 22)
                   zfrz(ji,jj)=zeps4+zlamb1*zsfrz
   
                   ! zfwflx is upward water flux
                   ! zhtflx is upward heat flux (out of ocean)
-                  ! compute the upward water and heat flux
+                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
                   zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
-                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
+                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
                END DO
             END DO
@@ -523,6 +519,6 @@
       IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas)
       ! 
-      CALL wrk_dealloc( jpi,jpj, zfrz  , zpress, zti     , zgammat, zgammas )
-      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b                   )
+      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  )
+      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
@@ -563,4 +559,5 @@
       IF      ( nn_gammablk == 0 ) THEN
       !! gamma is constant (specified in namelist)
+      !! ISOMIP formulation (Hunter et al, 2006)
          gt(:,:) = rn_gammat0
          gs(:,:) = rn_gammas0
@@ -569,6 +566,7 @@
       !! gamma is assume to be proportional to u* 
       !! Jenkins et al., 2010, JPO, p2298-2312
-
-      !! compute ustar
+      !! Adopted by Asay-Davis et al. (2015)
+
+      !! compute ustar (eq. 24)
          zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
 
@@ -604,7 +602,7 @@
                   zcoef = 0.5 / fse3w(ji,jj,ikt)
                   !                                            ! shear of horizontal velocity
-                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
+                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
                      &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
-                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
+                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
                      &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
                   !                                            ! richardson number (minimum value set to zero)
@@ -664,13 +662,16 @@
 
       REAL(wp) :: ze3, zhk
-      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl
-
-      INTEGER :: ji,jj,jk
-      INTEGER :: ikt,ikb
-
+      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl
+
+      INTEGER :: ji,jj,jk                   ! loop index
+      INTEGER :: ikt,ikb                    ! top and bottom index of the tbl
+ 
+      ! allocation
       CALL wrk_alloc( jpi,jpj, zhisf_tbl)
-
-      ! get first and last level of tbl
+      
+      ! initialisation
       varout(:,:)=0._wp
+   
+      ! compute U in the top boundary layer at T- point 
       IF (cptin == 'U') THEN
          DO jj = 1,jpj
@@ -681,6 +682,5 @@
 
             ! determine the deepest level influenced by the boundary layer
-            ! test on tmask useless ?????
-               DO jk = ikt, mbku(ji,jj)
+               DO jk = ikt+1, mbku(ji,jj)
                   IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
                END DO
@@ -706,4 +706,5 @@
       END IF
 
+      ! compute V in the top boundary layer at T- point 
       IF (cptin == 'V') THEN
          DO jj = 1,jpj
@@ -714,6 +715,5 @@
 
             ! determine the deepest level influenced by the boundary layer
-            ! test on tmask useless ?????
-               DO jk = ikt, mbkv(ji,jj)
+               DO jk = ikt+1, mbkv(ji,jj)
                   IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
@@ -739,4 +739,5 @@
       END IF
 
+      ! compute T in the top boundary layer at T- point 
       IF (cptin == 'T') THEN
          DO jj = 1,jpj
@@ -757,6 +758,9 @@
          END DO
       END IF
+
+      ! mask mean tbl value
       varout(:,:) = varout(:,:) * ssmask(:,:)
 
+      ! deallocation
       CALL wrk_dealloc( jpi,jpj, zhisf_tbl )      
 
@@ -795,5 +799,5 @@
 
                ! determine the deepest level influenced by the boundary layer
-               DO jk = ikt, mbkt(ji,jj)
+               DO jk = ikt+1, mbkt(ji,jj)
                   IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
@@ -825,77 +829,4 @@
       !
    END SUBROUTINE sbc_isf_div
-                        
-   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_init  ***
-      !!
-      !! ** Purpose :   Compute the in-situ temperature [Celcius]
-      !!
-      !! ** Method  :   
-      !!
-      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
-      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
-      REAL(wp) :: zt, zs, zp, zh, zq, zxk
-      INTEGER  :: ji, jj            ! dummy loop indices
-      !
-      CALL wrk_alloc( jpi,jpj, pti  )
-      ! 
-      DO jj=1,jpj
-         DO ji=1,jpi
-            zh = ppress(ji,jj)
-! Theta1
-            zt = ptem(ji,jj)
-            zs = psal(ji,jj)
-            zp = 0.0
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 0.5 * zxk
-            zq = zxk
-! Theta2
-            zp = zp + 0.5 * zh
-            zxk= zh*fsatg( zs, zt, zp )
-            zt = zt + 0.29289322 * ( zxk - zq )
-            zq = 0.58578644 * zxk + 0.121320344 * zq
-! Theta3
-            zxk= zh * fsatg( zs, zt, zp )
-            zt = zt + 1.707106781 * ( zxk - zq )
-            zq = 3.414213562 * zxk - 4.121320344 * zq
-! Theta4
-            zp = zp + 0.5 * zh
-            zxk= zh * fsatg( zs, zt, zp )
-            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, pti  )
-      !
-   END FUNCTION tinsitu
-   !
-   FUNCTION fsatg( pfps, pfpt, pfphp )
-      !!----------------------------------------------------------------------
-      !!                 ***  FUNCTION fsatg  ***
-      !!
-      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
-      !!
-      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
-      !! 
-      !! ** units      :   pressure        pfphp    decibars
-      !!                   temperature     pfpt     deg celsius (ipts-68)
-      !!                   salinity        pfps     (ipss-78)
-      !!                   adiabatic       fsatg    deg. c/decibar
-      !!----------------------------------------------------------------------
-      REAL(wp) :: pfps, pfpt, pfphp 
-      REAL(wp) :: fsatg
-      !
-      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
-        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
-        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
-        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
-        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
-      !
-    END FUNCTION fsatg
-    !!======================================================================
+   !!======================================================================
 END MODULE sbcisf
