Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 5944)
@@ -137,5 +137,5 @@
       DO jk=1,nb_ana
        DO ji=1,jpmax_harmo
-          IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
+          IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN
              name(jk) = ji
              EXIT
@@ -362,6 +362,6 @@
                X1=ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj)
-               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
+               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj)
+               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj)
             END DO
          END DO
@@ -492,5 +492,5 @@
             DO jj_sd = ji_sd, ninco
                zval2 = ABS(ztmp3(ji_sd,jj_sd))
-               IF( zval2.GE.zval1 )THEN
+               IF( zval2 >= zval1 )THEN
                   ipivot(ji_sd) = jj_sd
                   zval1         = zval2
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5944)
@@ -211,6 +211,6 @@
         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
-        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
+        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9   )              ! Salt content variation (psu*km3)
+        CALL iom_put( 'bgvolssh' , zdiff_v1  * 1.e-9   )              ! volume ssh variation (km3)  
         CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
         CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5944)
@@ -319,5 +319,5 @@
          WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
       ENDIF
-      IF ( nn_cla .EQ. 1 ) THEN
+      IF ( nn_cla == 1 ) THEN
          IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2 
             CONTINUE
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 5944)
@@ -332,5 +332,5 @@
          ! --------------------------------------------------
          IF( ln_vvl_ztilde ) THEN
-            IF( kt .GT. nit000 ) THEN
+            IF( kt > nit000 ) THEN
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
@@ -426,5 +426,5 @@
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          ! - ML - test: for the moment, stop simulation for too large e3_t variations
-         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
+         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
             IF( lk_mpp ) THEN
                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5944)
@@ -824,7 +824,7 @@
    END SUBROUTINE zgr_bot_level
 
-      SUBROUTINE zgr_top_level
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bot_level  ***
+   SUBROUTINE zgr_top_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_top_level  ***
       !!
       !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
@@ -1149,5 +1149,5 @@
                   gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
                END DO
-               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
+               IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
                DO jk = misfdep(ji,jj) + 1, jpk
                   gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
@@ -1212,8 +1212,11 @@
       INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices
       INTEGER  ::   ik, it               ! temporary integers
-      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
+      INTEGER  ::   icompt, ibtest       ! (ISF)
+      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF)
+      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF)
       REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t
       REAL(wp) ::   zmax             ! Maximum and minimum depth
-      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
+      REAL(wp) ::   zbathydiff       ! isf temporary scalar
+      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
@@ -1230,6 +1233,6 @@
 
       ! (ISF) compute misfdep
-      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
-      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
+      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
+      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
       END WHERE  
 
@@ -1242,19 +1245,17 @@
          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
       END DO 
-      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
+      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )
          risfdep(:,:) = 0. ; misfdep(:,:) = 1
       END WHERE
 
       ! remove very shallow ice shelf (less than ~ 10m if 75L)
-      IF ( cp_cfg .NE. "isomip" ) THEN
-         WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1)
-            misfdep = 0; risfdep = 0.0_wp;
-            mbathy  = 0; bathy   = 0.0_wp;
-         END WHERE
-         WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60)
-            misfdep = 0; risfdep = 0.0_wp;
-            mbathy  = 0; bathy   = 0.0_wp;
-         END WHERE
-      END IF
+      WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
+      WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp)
+         misfdep = 0; risfdep = 0.0_wp;
+         mbathy  = 0; bathy   = 0.0_wp;
+      END WHERE
  
 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
@@ -1262,5 +1263,5 @@
 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
       DO jl = 1, 10     
-         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
+         WHERE (bathy(:,:) == risfdep(:,:) )
             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
@@ -1271,21 +1272,20 @@
          ENDWHERE
          IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
-            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west 
+            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west 
             misfdep(jpi,:) = misfdep(  2  ,:) 
-         ENDIF
-
-         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
-            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
-            mbathy(jpi,:) = mbathy(  2  ,:)
+            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west
+            mbathy(jpi,:)  = mbathy(  2  ,:)
          ENDIF
 
@@ -1314,5 +1314,5 @@
                ! find the minimum change option:
                ! test bathy
-               IF (risfdep(ji,jj) .GT. 1) THEN
+               IF (risfdep(ji,jj) > 1) THEN
                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
                  &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
@@ -1320,6 +1320,6 @@
                  &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
  
-                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
-                     IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
+                     IF (zbathydiff <= zrisfdepdiff) THEN
                         bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
                         mbathy(ji,jj)= mbathy(ji,jj) + 1
@@ -1338,10 +1338,10 @@
                ! find the minimum change option:
                ! test bathy
-               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
+               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
                   zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
                     & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
                   zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
                     & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (zbathydiff <= zrisfdepdiff) THEN
                      mbathy(ji,jj) = mbathy(ji,jj) + 1
                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
@@ -1354,13 +1354,13 @@
          END DO
  
- ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+ ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
+               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
                   zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
                   zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (zbathydiff <= zrisfdepdiff) THEN
                      mbathy(ji,jj) = mbathy(ji,jj) + 1
                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
@@ -1376,22 +1376,24 @@
  
          IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
- ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 
+ ! point V misdep(ji,jj) == mbathy(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
+               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN
                   zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
                   zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (zbathydiff <= zrisfdepdiff) THEN
                      mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
                      bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
@@ -1406,23 +1408,25 @@
  
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF 
  
- ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+ ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
+               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
                   zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (zbathydiff <= zrisfdepdiff) THEN
                      mbathy(ji,jj) = mbathy(ji,jj) + 1
                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
@@ -1436,23 +1440,25 @@
  
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF 
  
- ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 
+ ! point U misfdep(ji,jj) == bathy(ji,jj+1) 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
+               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN
                   zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
                   zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
                       &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                  IF (zbathydiff <= zrisfdepdiff) THEN
                      mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
                      bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
@@ -1468,12 +1474,14 @@
  
          IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
       END DO
@@ -1485,13 +1493,13 @@
          DO jk = 2, jpk
             WHERE (misfdep==0) misfdep=jpk
-            zmask=0
-            WHERE (misfdep .LE. jk) zmask=1
+            zmask=0._wp
+            WHERE (misfdep <= jk) zmask=1
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  IF (misfdep(ji,jj) .EQ. jk) THEN
+                  IF (misfdep(ji,jj) == jk) THEN
                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
+                     IF (ibtest <= 1) THEN
                         risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk
                      END IF
                   END IF
@@ -1500,28 +1508,30 @@
          END DO
          WHERE (misfdep==jpk)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
          END WHERE
          IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
  
  ! remove single point "bay" on bathy coast line beneath an ice shelf'
          DO jk = jpk,1,-1
-            zmask=0
-            WHERE (mbathy .GE. jk ) zmask=1
+            zmask=0._wp
+            WHERE (mbathy >= jk ) zmask=1
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
+                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN
                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
+                     IF (ibtest <= 1) THEN
                         bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0
                      END IF
                   END IF
@@ -1530,15 +1540,17 @@
          END DO
          WHERE (mbathy==0)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
          END WHERE
          IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
+            mbathy(:,:)  = INT( zbathy(:,:) )
          ENDIF
  
@@ -1546,18 +1558,18 @@
          zmisfdep = misfdep
          zrisfdep = risfdep
-         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
+         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
+               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN
                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
                END IF
-               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
+               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN
                   misfdep(ji,jj) = ibtest
                   risfdep(ji,jj) = gdepw_1d(ibtest)
@@ -1567,11 +1579,13 @@
  
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
+
             CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+            CALL lbc_lnk( bathy,   'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
+            CALL lbc_lnk( zbathy,  'T', 1. )
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
@@ -1583,13 +1597,13 @@
                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
+               IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0
                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
+               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN
                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
                END IF
-               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
+               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN
                   mbathy(ji,jj) = ibtest
                   bathy(ji,jj)  = gdepw_1d(ibtest+1) 
@@ -1598,11 +1612,13 @@
          END DO
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
+
             CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+            CALL lbc_lnk( bathy,   'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
+            CALL lbc_lnk( zbathy,  'T', 1. )
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
@@ -1610,5 +1626,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
+               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
                END IF
@@ -1616,11 +1632,13 @@
          END DO
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
+
             CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+            CALL lbc_lnk( bathy,   'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
+            CALL lbc_lnk( zbathy,  'T', 1. )
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF 
@@ -1628,5 +1646,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
+               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
                   mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
                END IF
@@ -1634,9 +1652,11 @@
          END DO
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
             CALL lbc_lnk( zbathy, 'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
@@ -1646,5 +1666,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
+               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
                END IF
@@ -1652,9 +1672,11 @@
          END DO
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
             CALL lbc_lnk( zbathy, 'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
@@ -1664,15 +1686,17 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
-                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
+               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
+                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
                END IF
             END DO
          END DO
          IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
             CALL lbc_lnk( zbathy, 'T', 1. ) 
             misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
             zbathy(:,:) = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
@@ -1694,5 +1718,5 @@
 ! end check compatibility ice shelf/bathy
       ! remove very shallow ice shelf (less than ~ 10m if 75L)
-      WHERE (misfdep(:,:) <= 5)
+      WHERE (risfdep(:,:) <= 10._wp)
          misfdep = 1; risfdep = 0.0_wp;
       END WHERE
@@ -1814,5 +1838,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
       
-   END SUBROUTINE
+   END SUBROUTINE zgr_isf
 
    SUBROUTINE zgr_sco
@@ -2143,11 +2167,11 @@
       fsde3w(:,:,:) = gdep3w_0(:,:,:)
       !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      where (e3t_0   (:,:,:) == 0.0)  e3t_0(:,:,:)  = 1.0_wp
+      where (e3u_0   (:,:,:) == 0.0)  e3u_0(:,:,:)  = 1.0_wp
+      where (e3v_0   (:,:,:) == 0.0)  e3v_0(:,:,:)  = 1.0_wp
+      where (e3f_0   (:,:,:) == 0.0)  e3f_0(:,:,:)  = 1.0_wp
+      where (e3w_0   (:,:,:) == 0.0)  e3w_0(:,:,:)  = 1.0_wp
+      where (e3uw_0  (:,:,:) == 0.0)  e3uw_0(:,:,:) = 1.0_wp
+      where (e3vw_0  (:,:,:) == 0.0)  e3vw_0(:,:,:) = 1.0_wp
 
 #if defined key_agrif
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5944)
@@ -235,5 +235,5 @@
                      &                              * (1._wp - tmask(ji,jj,jk))
                END DO
-               IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
                                   &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
             END DO
@@ -521,10 +521,10 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce, zpshpi, zpshpj
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,  2, ztstop) 
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj)
-      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce, zpshpi, zpshpj) 
+      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce )
       !
       ! Local constant initialization
@@ -605,5 +605,5 @@
       CALL wrk_dealloc( jpi,jpj,2  , ztstop)
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)
-      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce, zpshpi, zpshpj)
+      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce )
       !
    END SUBROUTINE hpg_isf
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5944)
@@ -100,5 +100,5 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   iku, ikv     ! local integers
+      INTEGER  ::   ikt          ! local integers
 #if ! defined key_dynspg_flt
       REAL(wp) ::   z2dt         ! temporary scalar
@@ -268,9 +268,9 @@
                DO jj = 1,jpj
                   DO ji = 1,jpi
-                     jk = mikt(ji,jj)
-                     fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                  & 
+                     ikt = mikt(ji,jj)
+                     fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0                  & 
                                        &                          * ( (emp_b(ji,jj) - emp(ji,jj) ) &
                                        &                            - (rnf_b(ji,jj) - rnf(ji,jj) ) &    
-                                       &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj,jk)
+                                       &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj,ikt)
                   END DO
                END DO
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5944)
@@ -934,5 +934,5 @@
       !
       IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
-         IF ( Agrif_NbStepint().EQ.0 ) THEN
+         IF ( Agrif_NbStepint() == 0 ) THEN
             ub2_i_b(:,:) = 0.e0
             vb2_i_b(:,:) = 0.e0
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5944)
@@ -86,5 +86,5 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi
-               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
@@ -193,5 +193,5 @@
          DO jj = 2, jpj                   
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj,jk) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5944)
@@ -586,5 +586,5 @@
                   !
                   jk = nmln(ji+ip,jj) + 1
-                  IF( jk .GT. mbkt(ji+ip,jj) ) THEN  !ML reaches bottom
+                  IF( jk >  mbkt(ji+ip,jj) ) THEN  !ML reaches bottom
                     zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp
                   ELSE
@@ -596,5 +596,5 @@
                   !
                   jk = nmln(ji,jj+jp) + 1
-                  IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom
+                  IF( jk >  mbkt(ji,jj+jp) ) THEN  !ML reaches bottom
                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp
                   ELSE
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5944)
@@ -91,5 +91,5 @@
          IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
          !
-         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
+         area = glob_sum( e12t(:,:) * tmask(:,:,1))           ! interior global domain surface
          ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
          ! and in case of no melt, it can generate HSSW.
@@ -108,5 +108,5 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
@@ -133,5 +133,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( e12t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
@@ -159,8 +159,8 @@
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
             !
-            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
-            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
+            zsurf_neg = glob_sum( e12t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp 
+            zsurf_pos = glob_sum( e12t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
-            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
+            z_fwf     = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
@@ -172,9 +172,9 @@
             ENDIF
             !
-            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
+            zsum_fwf   = glob_sum( e12t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
 !!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
             z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
             !                                                  ! weight to respect erp field 2D structure 
-            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
+            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e12t(:,:) )
             z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
             !                                                  ! final correction term to apply
@@ -191,10 +191,10 @@
                IF( z_fwf < 0._wp ) THEN
                   WRITE(numout,*)'   z_fwf < 0'
-                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv'
                ELSE
                   WRITE(numout,*)'   z_fwf >= 0'
-                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv'
                ENDIF
-               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
+               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e12t(:,:) )*1.e-9,' Sv'
                WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
                WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
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 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5944)
@@ -56,21 +56,23 @@
                                                                                           !: (first wet level and last level include in the tbl)
 #else
-   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base
 #endif
 
 
-   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
-   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! phycst ?
 
 !: 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_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           
-   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
+   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
+   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           
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the runoff file to be read ??
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the runoff file to be read ??
+   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the runoff file to be read ??
    
    !! * Substitutions
@@ -85,145 +87,41 @@
  
   SUBROUTINE sbc_isf(kt)
-    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), 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
-      !
       !!---------------------------------------------------------------------
-      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf      , &
-                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
-      !
-      ! allocation
-      CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
+      !!                  ***  ROUTINE sbc_isf  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt                   ! ocean time step
+      !
+      INTEGER               :: ji, jj               ! loop index
+      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
+      !!---------------------------------------------------------------------
       !
       !                                         ! ====================== !
       IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
          !                                      ! ====================== !
-         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
-         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
-
-         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
-         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
-         IF(lwm) WRITE ( numond, namsbc_isf )
-
-
-         IF ( lwp ) WRITE(numout,*)
-         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
-         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
-         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
-         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
-         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
-         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
-         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
-         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
-         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+         CALL sbc_isf_init
+      !                                         ! ---------------------------------------- !
+      ELSE                                      !          Swap of forcing fields          !
+         !                                      ! ---------------------------------------- !
+         fwfisf_b  (:,:  ) = fwfisf  (:,:  )    ! Swap the ocean forcing fields except at nit000
+         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)    ! where before fields are set at the end of the routine
          !
-         ! Allocate public variable
-         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
-         !
-         ! initialisation
-         qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp
-         risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp
-         !
-         ! define isf tbl tickness, top and bottom indice
-         IF      (nn_isf == 1) THEN
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
-            ALLOCATE( sf_rnfisf(1), STAT=ierror )
-            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-
-            !: read effective lenght (BG03)
-            IF (nn_isf == 2) THEN
-               cvarLeff = 'soLeff' 
-               CALL iom_open( sn_Leff_isf%clname, inum )
-               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
-               CALL iom_close(inum)
-               !
-               risfLeff = risfLeff*1000.0_wp           !: convertion in m
-            END IF
-
-           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
-            CALL iom_open( sn_depmax_isf%clname, inum )
-            cvarhisf = TRIM(sn_depmax_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
-            CALL iom_close(inum)
-            !
-            CALL iom_open( sn_depmin_isf%clname, inum )
-            cvarzisf = TRIM(sn_depmin_isf%clvar)
-            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
-            CALL iom_close(inum)
-            !
-            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
-
-           !! compute first level of the top boundary layer
-           DO ji = 1, jpi
-              DO jj = 1, jpj
-                  jk = 2
-                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
-                  misfkt(ji,jj) = jk-1
-               END DO
-            END DO
-
-         ELSE IF ( nn_isf == 4 ) THEN
-            ! as in nn_isf == 1
-            rhisf_tbl(:,:) = rn_hisf_tbl
-            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
-            
-            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), STAT=ierror )
-            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(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' )
-         END IF
-         
-         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
-
-         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               ikt = misfkt(ji,jj)
-               ikb = misfkt(ji,jj)
-               ! thickness of boundary layer at least the top level thickness
-               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
-
-               ! determine the deepest level influenced by the boundary layer
-               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
-               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
-               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
-               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
-
-               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
-               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
-            END DO
-         END DO
-         
       END IF
 
-      !                                            ! ---------------------------------------- !
-      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
-         !                                         ! ---------------------------------------- !
-         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
-         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
-         !
-      ENDIF
-
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
-
-
-         ! compute salf and heat flux
-         IF (nn_isf == 1) THEN
-            ! realistic ice shelf formulation
+         ! allocation
+         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
+
+         ! compute salt and heat flux
+         SELECT CASE ( nn_isf )
+         CASE ( 1 )    ! realistic ice shelf formulation
             ! compute T/S/U/V for the top boundary layer
             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
@@ -239,23 +137,21 @@
             CALL sbc_isf_cav (kt)
 
-         ELSE IF (nn_isf == 2) THEN
-            ! Beckmann and Goosse parametrisation 
+         CASE ( 2 )    ! Beckmann and Goosse parametrisation 
             stbl(:,:)   = soce
             CALL sbc_isf_bg03(kt)
 
-         ELSE IF (nn_isf == 3) THEN
-            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
+         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting) 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat flux
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
             stbl(:,:)   = soce
 
-         ELSE IF (nn_isf == 4) THEN
-            ! specified fwf and heat flux forcing beneath the ice shelf
+         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
             fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting)
-            qisf(:,:)   = fwfisf(:,:) * lfusisf                ! heat flux
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
             stbl(:,:)   = soce
-         END IF
+
+         END SELECT
 
          ! compute tsc due to isf
@@ -276,6 +172,6 @@
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
          CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
-         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
-         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
+         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)           ,'T',1.)
 
          IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
@@ -290,13 +186,17 @@
                risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
             END IF
-         ENDIF
+         END IF
          ! 
          ! output
-         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf)
+         ! JMM : iom_use not necessary here, qisf and fwfisf are always computed. 
+         !   If not required in iodef.xml, iom_put does not do anything
+!        IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
+!        IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf)
+         CALL iom_put('qisf'  , qisf)
+         CALL iom_put('fwfisf', fwfisf)
+
+         ! deallocation
+         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
       END IF
-
-      ! deallocation
-      CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
   
   END SUBROUTINE sbc_isf
@@ -315,78 +215,205 @@
                &    STAT= sbc_isf_alloc )
          !
-         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
          IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
          !
-      ENDIF
+      END IF
   END FUNCTION
 
+  SUBROUTINE sbc_isf_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_init  ***
+      !!
+      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER               :: ji, jj, jk           ! loop index
+      INTEGER               :: ik                   ! current level 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
+      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
+      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_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_isf )
+
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
+      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
+      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
+      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
+      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
+      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+      !
+      ! Allocate public variable
+      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
+      !
+      ! initialisation
+      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp
+      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp
+      !
+      ! define isf tbl tickness, top and bottom indice
+      SELECT CASE ( nn_isf )
+      CASE ( 1 ) 
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+
+      CASE ( 2 , 3 )
+         ALLOCATE( sf_rnfisf(1), STAT=ierror )
+         ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+
+         !  read effective lenght (BG03)
+         IF (nn_isf == 2) THEN
+            CALL iom_open( sn_Leff_isf%clname, inum )
+            cvarLeff = TRIM(sn_Leff_isf%clvar)
+            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
+            CALL iom_close(inum)
+            !
+            risfLeff = risfLeff*1000.0_wp           !: convertion in m
+         END IF
+
+         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
+         CALL iom_open( sn_depmax_isf%clname, inum )
+         cvarhisf = TRIM(sn_depmax_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
+         CALL iom_close(inum)
+         !
+         CALL iom_open( sn_depmin_isf%clname, inum )
+         cvarzisf = TRIM(sn_depmin_isf%clvar)
+         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
+         CALL iom_close(inum)
+         !
+         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
+
+         !! compute first level of the top boundary layer
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+                ik = 2
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                misfkt(ji,jj) = ik-1
+            END DO
+         END DO
+
+      CASE ( 4 ) 
+         ! as in nn_isf == 1
+         rhisf_tbl(:,:) = rn_hisf_tbl
+         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         
+         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
+         ALLOCATE( sf_fwfisf(1), STAT=ierror )
+         ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(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' )
+
+      END SELECT
+         
+      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
+
+      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = misfkt(ji,jj)
+            ikb = misfkt(ji,jj)
+            ! thickness of boundary layer at least the top level thickness
+            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
+
+            ! determine the deepest level influenced by the boundary layer
+            DO jk = ikt+1, mbkt(ji,jj)
+               IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+            END DO
+            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
+            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
+            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
+
+            zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
+            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+         END DO
+      END DO
+
+  END SUBROUTINE sbc_isf_init
+
   SUBROUTINE sbc_isf_bg03(kt)
-   !!==========================================================================
-   !!                 *** SUBROUTINE sbcisf_bg03  ***
-   !! add net heat and fresh water flux from ice shelf melting
-   !! into the adjacent ocean using the parameterisation by
-   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
-   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
-   !!  (hereafter BG)
-   !!==========================================================================
-   !!----------------------------------------------------------------------
-   !!   sbc_isf_bg03      : routine called from sbcmod
-   !!----------------------------------------------------------------------
-   !!
-   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
-   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
-   !!
-   !! History :
-   !!      !  06-02  (C. Wang) Original code
-   !!----------------------------------------------------------------------
-
-    INTEGER, INTENT ( in ) :: kt
-
-    INTEGER :: ji, jj, jk !temporary integer
-
-    REAL(wp) :: zt_sum    ! sum of the temperature between 200m and 600m
-    REAL(wp) :: zt_ave    ! averaged temperature between 200m and 600m
-    REAL(wp) :: zt_frz    ! freezing point temperature at depth z
-    REAL(wp) :: zpress    ! pressure to compute the freezing point in depth
-    
-    !!----------------------------------------------------------------------
-    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
-     !
-
-    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
-    DO ji = 1, jpi
-       DO jj = 1, jpj
-          jk = misfkt(ji,jj)
-          !! Initialize arrays to 0 (each step)
-          zt_sum = 0.e0_wp
-          IF ( jk .GT. 1 ) THEN
-    ! 1. -----------the average temperature between 200m and 600m ---------------------
-             DO jk = misfkt(ji,jj),misfkb(ji,jj)
-             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
-             ! after verif with UNESCO, wrong sign in BG eq. 2
-             ! Calculate freezing temperature
-                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
-             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
-    
-    ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
-          ! For those corresponding to zonal boundary    
-             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
-                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_bg03  ***
+      !!
+      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
+      !!          into the adjacent ocean
+      !!
+      !! ** Method  :   See reference
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!         (hereafter BG)
+      !! History :
+      !!         06-02  (C. Wang) Original code
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT ( in ) :: kt
+      !
+      INTEGER  :: ji, jj, jk ! dummy loop index
+      INTEGER  :: ik         ! current level
+      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
+      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
+      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
+      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
+      !!----------------------------------------------------------------------
+
+      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
+      !
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+            ik = misfkt(ji,jj)
+            !! Initialize arrays to 0 (each step)
+            zt_sum = 0.e0_wp
+            IF ( ik > 1 ) THEN
+               ! 1. -----------the average temperature between 200m and 600m ---------------------
+               DO jk = misfkt(ji,jj),misfkb(ji,jj)
+                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
+                  ! after verif with UNESCO, wrong sign in BG eq. 2
+                  ! Calculate freezing temperature
+                  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
+               END DO
+               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
+               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
+               ! For those corresponding to zonal boundary    
+               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
+                           & * r1_e12t(ji,jj) * tmask(ji,jj,jk)
              
-             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                  
-             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
-             !add to salinity trend
-          ELSE
-             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
-          END IF
-       ENDDO
-    ENDDO
-    !
-    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                  
+               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
+               !add to salinity trend
+            ELSE
+               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
+            END IF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+
   END SUBROUTINE sbc_isf_bg03
 
-   SUBROUTINE sbc_isf_cav( kt )
+  SUBROUTINE sbc_isf_cav( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_isf_cav  ***
@@ -403,7 +430,6 @@
       INTEGER, INTENT(in)          ::   kt         ! ocean time step
       !
-      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
-      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
-      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   nit
       REAL(wp) ::   zlamb1, zlamb2, zlamb3
       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
@@ -411,14 +437,14 @@
       REAL(wp) ::   zeps = 1.e-20_wp        
       REAL(wp) ::   zerr
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      INTEGER  ::   nit
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
+      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
       LOGICAL  ::   lit
       !!---------------------------------------------------------------------
-      !
       ! 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
+      zlamb1 =-0.0573_wp
+      zlamb2 = 0.0832_wp
+      zlamb3 =-7.53e-08_wp * grav * rau0
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
@@ -427,13 +453,13 @@
 
       ! initialisation
-      zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0;
-      zhtflx (:,:)=0.0_wp     ; zhtflx_b(:,:)=0.0_wp    ;
-      zfwflx (:,:)=0.0_wp
+      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
+      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp    
+      zfwflx (:,:) = 0.0_wp
 
       ! 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 formulation (2 equations) for volume flux (Hunter et al., 2006)
+         SELECT CASE ( nn_isfblk )
+         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
             ! Calculate freezing temperature
             CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
@@ -446,5 +472,5 @@
                DO ji = 1, jpi
                   zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
-                  zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf
+                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
                END DO
             END DO
@@ -454,6 +480,5 @@
             fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
 
-         ELSE IF (nn_isfblk == 2 ) THEN
-            ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
+         CASE ( 2 )  ! 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)
@@ -464,20 +489,20 @@
                DO ji = 1, jpi
                   ! compute coeficient to solve the 2nd order equation
-                  zeps1=rcp*rau0*zgammat(ji,jj)
-                  zeps2=lfusisf*rau0*zgammas(ji,jj)
-                  zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
-                  zeps4=zlamb2+zlamb3*risfdep(ji,jj)
-                  zeps6=zeps4-ttbl(ji,jj)
-                  zeps7=zeps4-tsurf
-                  zaqe=zlamb1 * (zeps1 + zeps3)
-                  zaqer=0.5/MIN(zaqe,-zeps)
-                  zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
-                  zcqe=zeps2*stbl(ji,jj)
-                  zdis=zbqe*zbqe-4.0*zaqe*zcqe               
+                  zeps1 = rcp*rau0*zgammat(ji,jj)
+                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
+                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
+                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
+                  zeps6 = zeps4-ttbl(ji,jj)
+                  zeps7 = zeps4-tsurf
+                  zaqe  = zlamb1 * (zeps1 + zeps3)
+                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
+                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
+                  zcqe  = zeps2*stbl(ji,jj)
+                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
 
                   ! Presumably zdis can never be negative because gammas is very small compared to gammat
                   ! compute s freeze
                   zsfrz=(-zbqe-SQRT(zdis))*zaqer
-                  IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
+                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
 
                   ! compute t freeze (eq. 22)
@@ -496,26 +521,28 @@
             fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
 
-         ENDIF
+         END SELECT
 
          ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
-         IF ( nn_gammablk .LT. 2 ) THEN ; lit = .FALSE.
+         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
          ELSE                           
             ! check total number of iteration
-            IF (nit .GE. 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-            ELSE                   ; nit = nit + 1
-            ENDIF
+            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            ELSE                 ; nit = nit + 1
+            END IF
 
             ! compute error between 2 iterations
             ! if needed save gammat and compute zhtflx_b for next iteration
             zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
-            IF ( zerr .LE. 0.01 ) THEN ; lit = .FALSE.
-            ELSE                       ; zhtflx_b(:,:) = zhtflx(:,:)
-            ENDIF
+            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
+            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
+            END IF
          END IF
       END DO
       !
-      ! output
-      IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat)
-      IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas)
+      ! output  see JMM comment above
+!     IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat)
+!     IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas)
+      CALL iom_put('isfgammat', zgammat)
+      CALL iom_put('isfgammas', zgammas)
       ! 
       CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  )
@@ -526,5 +553,5 @@
    END SUBROUTINE sbc_isf_cav
 
-   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf )
+   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -535,9 +562,9 @@
       !!                Jenkins et al., 2010, JPO, p2298-2312
       !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(out) :: gt, gs
-      REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf
-
+      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs
+      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf
+      !
       INTEGER  :: ikt                        
-      INTEGER  :: ji,jj                      ! loop index
+      INTEGER  :: ji, jj                     ! loop index
       REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
@@ -551,44 +578,42 @@
       REAL(wp) :: zdep
       REAL(wp) :: zeps = 1.0e-20_wp    
-      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
-      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
       CALL wrk_alloc( jpi,jpj, zustar )
       !
-      IF      ( nn_gammablk == 0 ) THEN
-      !! gamma is constant (specified in namelist)
-      !! ISOMIP formulation (Hunter et al, 2006)
-         gt(:,:) = rn_gammat0
-         gs(:,:) = rn_gammas0
-
-      ELSE IF ( nn_gammablk == 1 ) THEN
-      !! gamma is assume to be proportional to u* 
-      !! Jenkins et al., 2010, JPO, p2298-2312
-      !! Adopted by Asay-Davis et al. (2015)
-
-      !! compute ustar (eq. 24)
+      SELECT CASE ( nn_gammablk )
+      CASE ( 0 ) ! gamma is constant (specified in namelist)
+         !! ISOMIP formulation (Hunter et al, 2006)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+
+      CASE ( 1 ) ! gamma is assume to be proportional to u*
+         !! Jenkins et al., 2010, JPO, p2298-2312
+         !! Adopted by Asay-Davis et al. (2015)
+
+         !! compute ustar (eq. 24)
          zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
 
-      !! Compute gammats
-         gt(:,:) = zustar(:,:) * rn_gammat0
-         gs(:,:) = zustar(:,:) * rn_gammas0
+         !! Compute gammats
+         pgt(:,:) = zustar(:,:) * rn_gammat0
+         pgs(:,:) = zustar(:,:) * rn_gammas0
       
-      ELSE IF ( nn_gammablk == 2 ) THEN
-      !! gamma depends of stability of boundary layer
-      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
-      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
-      !! compute ustar
+      CASE ( 2 ) ! gamma depends of stability of boundary layer
+         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
+         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+         !! compute ustar
          zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
 
-      !! compute Pr and Sc number (can be improved)
-         zPr =   13.8
-         zSc = 2432.0
-
-      !! compute gamma mole
-         zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
-         zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
-
-      !! compute gamma
+         !! compute Pr and Sc number (can be improved)
+         zPr =   13.8_wp
+         zSc = 2432.0_wp
+
+         !! compute gamma mole
+         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+
+         !! compute gamma
          DO ji=2,jpi
             DO jj=2,jpj
@@ -596,9 +621,9 @@
 
                IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
-                  gt = rn_gammat0
-                  gs = rn_gammas0
+                  pgt = rn_gammat0
+                  pgs = rn_gammas0
                ELSE
-      !! compute Rc number (as done in zdfric.F90)
-                  zcoef = 0.5 / fse3w(ji,jj,ikt)
+                  !! compute Rc number (as done in zdfric.F90)
+                  zcoef = 0.5_wp / fse3w(ji,jj,ikt)
                   !                                            ! shear of horizontal velocity
                   zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
@@ -609,5 +634,5 @@
                   zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
 
-      !! compute bouyancy 
+                  !! compute bouyancy 
                   zts(jp_tem) = ttbl(ji,jj)
                   zts(jp_sal) = stbl(ji,jj)
@@ -616,8 +641,8 @@
                   CALL eos_rab( zts, zdep, zab )
                   !
-      !! compute length scale 
-                  zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      !! compute Monin Obukov Length
+                  !! compute length scale 
+                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+                  !! compute Monin Obukov Length
                   ! Maximum boundary layer depth
                   zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp
@@ -626,28 +651,28 @@
                   zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
 
-      !! compute eta* (stability parameter)
+                  !! compute eta* (stability parameter)
                   zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
 
-      !! compute the sublayer thickness
+                  !! compute the sublayer thickness
                   zhnu = 5 * znu / zustar(ji,jj)
 
-      !! compute gamma turb
+                  !! compute gamma turb
                   zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
                   &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
 
-      !! compute gammats
-                  gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
-                  gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
+                  !! compute gammats
+                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
                END IF
             END DO
          END DO
-         CALL lbc_lnk(gt(:,:),'T',1.)
-         CALL lbc_lnk(gs(:,:),'T',1.)
-      END IF
+         CALL lbc_lnk(pgt(:,:),'T',1.)
+         CALL lbc_lnk(pgs(:,:),'T',1.)
+      END SELECT
       CALL wrk_dealloc( jpi,jpj, zustar )
 
-   END SUBROUTINE
-
-   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
+   END SUBROUTINE sbc_isf_gammats
+
+   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_tbl  ***
@@ -656,89 +681,84 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
-      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
-      
-      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
-
+      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin
+      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout
+      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out
+      !
       REAL(wp) :: ze3, zhk
       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
- 
+      INTEGER :: ji, jj, jk                  ! loop index
+      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl
+      !!----------------------------------------------------------------------
       ! allocation
       CALL wrk_alloc( jpi,jpj, zhisf_tbl)
       
       ! initialisation
-      varout(:,:)=0._wp
+      pvarout(:,:)=0._wp
    
-      ! compute U in the top boundary layer at T- point 
-      IF (cptin == 'U') THEN
+      SELECT CASE ( cd_ptin )
+      CASE ( 'U' ) ! compute U in the top boundary layer at T- point 
          DO jj = 1,jpj
             DO ji = 1,jpi
                ikt = miku(ji,jj) ; ikb = miku(ji,jj)
-            ! thickness of boundary layer at least the top level thickness
+               ! thickness of boundary layer at least the top level thickness
                zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
 
-            ! determine the deepest level influenced by the boundary layer
+               ! determine the deepest level influenced by the boundary layer
                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
+                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
-            ! level fully include in the ice shelf boundary layer
+               ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
                   ze3 = fse3u_n(ji,jj,jk)
-                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
-               END DO
-
-            ! level partially include in ice shelf boundary layer 
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
                zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
-               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
          END DO
          DO jj = 2,jpj
             DO ji = 2,jpi
-               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj))
-            END DO
-         END DO
-         CALL lbc_lnk(varout,'T',-1.)
-      END IF
-
-      ! compute V in the top boundary layer at T- point 
-      IF (cptin == 'V') THEN
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+      
+      CASE ( 'V' ) ! compute V in the top boundary layer at T- point 
          DO jj = 1,jpj
             DO ji = 1,jpi
                ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
-           ! thickness of boundary layer at least the top level thickness
+               ! thickness of boundary layer at least the top level thickness
                zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
 
-            ! determine the deepest level influenced by the boundary layer
+               ! determine the deepest level influenced by the boundary layer
                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
+                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
-            ! level fully include in the ice shelf boundary layer
+               ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
                   ze3 = fse3v_n(ji,jj,jk)
-                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
-               END DO
-
-            ! level partially include in ice shelf boundary layer 
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
                zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
-               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
          END DO
          DO jj = 2,jpj
             DO ji = 2,jpi
-               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1))
-            END DO
-         END DO
-         CALL lbc_lnk(varout,'T',-1.)
-      END IF
-
-      ! compute T in the top boundary layer at T- point 
-      IF (cptin == 'T') THEN
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+
+      CASE ( 'T' ) ! compute T in the top boundary layer at T- point 
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -746,19 +766,19 @@
                ikb = misfkb(ji,jj)
 
-            ! level fully include in the ice shelf boundary layer
+               ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
                   ze3 = fse3t_n(ji,jj,jk)
-                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
-               END DO
-
-            ! level partially include in ice shelf boundary layer 
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
                zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
-               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
-            END DO
-         END DO
-      END IF
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+      END SELECT
 
       ! mask mean tbl value
-      varout(:,:) = varout(:,:) * ssmask(:,:)
+      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
 
       ! deallocation
@@ -780,6 +800,6 @@
       !! ** Action  :   phdivn   decreased by the runoff inflow
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
+      ! 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt, ikb 
@@ -790,5 +810,5 @@
       zfact   = 0.5_wp
       !
-      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
+      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -800,5 +820,5 @@
                ! determine the deepest level influenced by the boundary layer
                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
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -820,9 +840,9 @@
                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
+                    &              * 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) 
+                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
          END DO
       END DO
Index: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5921)
+++ branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5944)
@@ -141,10 +141,10 @@
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
             rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
                rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
             END WHERE
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp
             END WHERE
          ELSE                                                        ! use SST as runoffs temperature
