Index: /NEMO/trunk/src/OCE/ISF/isfcav.F90
===================================================================
--- /NEMO/trunk/src/OCE/ISF/isfcav.F90	(revision 15083)
+++ /NEMO/trunk/src/OCE/ISF/isfcav.F90	(revision 15084)
@@ -22,5 +22,5 @@
    USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
    !
-   USE oce      , ONLY: ts                              ! ocean tracers
+   USE oce      , ONLY: ts, uu, vv, rn2                 ! ocean dynamics and tracers
    USE par_oce  , ONLY: jpi,jpj                         ! ocean space and time domain
    USE phycst   , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp  ! physical constants
@@ -31,4 +31,5 @@
    USE fldread        ! read input field at current time step
    USE lbclnk         ! lbclnk
+   USE lib_mpp        ! MPP library
 
    IMPLICIT NONE
@@ -38,4 +39,6 @@
    PUBLIC   isf_cav, isf_cav_init ! routine called in isfmlt
 
+   !! * Substitutions   
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -71,8 +74,9 @@
       !!---------------------------------------------------------------------
       LOGICAL :: lit
-      INTEGER :: nit
+      INTEGER :: nit, ji, jj, ikt
       REAL(wp) :: zerr
+      REAL(wp) :: zcoef, zdku, zdkv
       REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh  ! heat fluxes
-      REAL(wp), DIMENSION(jpi,jpj) :: zqh_b                    !
+      REAL(wp), DIMENSION(jpi,jpj) :: zqh_b, zRc               !
       REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas         ! exchange coeficient
       REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl             ! temp. and sal. in top boundary layer
@@ -91,4 +95,19 @@
          zqoce(:,:) = -pqfwf(:,:) * rLfusisf ! 
          zqh_b(:,:) =  ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
+
+         DO_2D( 0, 0, 0, 0 )
+            ikt = mikt(ji,jj)
+            ! compute Rc number (as done in zdfric.F90)
+!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
+            zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
+            !                                            ! shear of horizontal velocity
+            zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
+               &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
+            zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
+               &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
+            !                                            ! richardson number (minimum value set to zero)
+            zRc(ji,jj) = MAX(rn2(ji,jj,ikt+1), 1.e-20_wp) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
+         END_2D
+         CALL lbc_lnk( 'isfmlt', zRc, 'T', 1._wp )
       ENDIF
       !
@@ -100,5 +119,5 @@
          ! useless if melt specified
          IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
-            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf,  &
+            CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce  , pqfwf, zRc,  &
                &                                    zgammat, zgammas )
          END IF
@@ -115,5 +134,9 @@
          CASE ( 'vel_stab' )
             ! compute error between 2 iterations
-            zerr = MAXVAL(ABS(zqhc(:,:)+zqoce(:,:) - zqh_b(:,:)))
+            zerr = 0._wp
+            DO_2D( 0, 0, 0, 0 )
+               zerr = MAX( zerr, ABS(zqhc(ji,jj)+zqoce(ji,jj) - zqh_b(ji,jj)) )
+            END_2D
+            CALL mpp_max( 'isfcav', zerr )   ! max over the global domain
             !
             ! define if iteration needed
@@ -130,23 +153,22 @@
       END DO
       !
-      ! compute heat and water flux ( > 0 from isf to oce)
-      pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:)
-      zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:)
-      zqhc (:,:) = zqhc(:,:)  * mskisf_cav(:,:)
-      !
-      ! compute heat content flux ( > 0 from isf to oce)
-      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
-      !
-      ! total heat flux ( > 0 from isf to oce)
-      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
-      !
-      ! lbclnk on melt
-      CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ! compute heat and water flux ( > 0 from isf to oce)
+         pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_cav(ji,jj)
+         zqoce(ji,jj) = zqoce(ji,jj) * mskisf_cav(ji,jj)
+         zqhc (ji,jj) = zqhc(ji,jj)  * mskisf_cav(ji,jj)
+         !
+         ! compute heat content flux ( > 0 from isf to oce)
+         zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf    ! 2d latent heat flux (W/m2)
+         !
+         ! total heat flux ( > 0 from isf to oce)
+         zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
+         !
+         ! set temperature content
+         ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
+      END_2D
       !
       ! output fluxes
       CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
-      !
-      ! set temperature content
-      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
       !
       ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
Index: /NEMO/trunk/src/OCE/ISF/isfcavgam.F90
===================================================================
--- /NEMO/trunk/src/OCE/ISF/isfcavgam.F90	(revision 15083)
+++ /NEMO/trunk/src/OCE/ISF/isfcavgam.F90	(revision 15084)
@@ -14,5 +14,5 @@
    USE isftbl  , ONLY: isf_tbl
 
-   USE oce     , ONLY: uu, vv, rn2         ! ocean dynamics and tracers
+   USE oce     , ONLY: uu, vv              ! ocean dynamics
    USE phycst  , ONLY: grav, vkarmn        ! physical constant
    USE eosbn2  , ONLY: eos_rab             ! equation of state
@@ -44,5 +44,5 @@
    !!-----------------------------------------------------------------------------------------------------
    !
-   SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
+   SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat and fwf flux
@@ -57,4 +57,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pRc             ! Richardson number
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)                :: zutbl, zvtbl    ! top boundary layer velocity
@@ -94,7 +95,7 @@
          pgs(:,:) = rn_gammas0
       CASE ( 'vel' ) ! gamma is proportional to u*
-         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
+         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,                    pgt, pgs )
       CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs )
       CASE DEFAULT
          CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
@@ -153,6 +154,6 @@
    END SUBROUTINE gammats_vel
 
-   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
-      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s]
+   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, pRc, &  ! <<== in
+      &                                                                     pgt  , pgs         )  ! ==>> out gammats [m/s]
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange coefficient 
@@ -171,9 +172,10 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl   ! velocity in the losch top boundary layer
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl   ! tracer   in the losch top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pRc            ! Richardson number
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj                     ! loop index
       INTEGER  :: ikt                        ! local integer
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
-      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
+      REAL(wp) :: zPr, zSc                   ! Prandtl and Scmidth number 
       REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
       REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
@@ -190,12 +192,4 @@
       !!---------------------------------------------------------------------
       !
-      ! compute ustar
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )
-      END_2D
-      !
-      ! output ustar
-      CALL iom_put('isfustar',zustar(:,:))
-      !
       ! compute Pr and Sc number (eq ??)
       zPr =   13.8_wp
@@ -207,6 +201,10 @@
       !
       ! compute gamma
-      DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+
          ikt = mikt(ji,jj)
+
+         ! compute ustar
+         zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )
 
          IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
@@ -214,15 +212,4 @@
             pgs(ji,jj) = rn_gammas0
          ELSE
-            ! compute Rc number (as done in zdfric.F90)
-!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
-            zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
-            !                                            ! shear of horizontal velocity
-            zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
-               &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
-            zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
-               &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
-            !                                            ! richardson number (minimum value set to zero)
-            zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
-            
             ! compute bouyancy 
             zts(jp_tem) = pttbl(ji,jj)
@@ -244,11 +231,12 @@
             !
             ! compute eta* (stability parameter) (Eq ??)
-            zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
+            zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * zustar(ji,jj) &
+               &                                        / MAX( 1.e-20, ABS(ff_t(ji,jj)) * zmols * pRc(ji,jj) ) )))
             !
             ! compute the sublayer thickness (Eq ??)
-            zhnu = 5 * znu / zustar(ji,jj)
+            zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) )
             !
             ! compute gamma turb (Eq ??)
-            zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
+            zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) &
                &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
             !
@@ -258,4 +246,6 @@
          END IF
       END_2D
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
 
    END SUBROUTINE gammats_vel_stab
Index: /NEMO/trunk/src/OCE/ISF/isfpar.F90
===================================================================
--- /NEMO/trunk/src/OCE/ISF/isfpar.F90	(revision 15083)
+++ /NEMO/trunk/src/OCE/ISF/isfpar.F90	(revision 15084)
@@ -30,5 +30,4 @@
    USE iom            ! I/O library
    USE fldread        ! read input field at current time step
-   USE lbclnk         ! lbc_lnk
 
    IMPLICIT NONE
@@ -37,4 +36,6 @@
    PUBLIC   isf_par, isf_par_init
 
+   !! * Substitutions   
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -67,4 +68,5 @@
       INTEGER, INTENT(in) ::   Kmm                                    ! ocean time level index
       !!---------------------------------------------------------------------
+      INTEGER ::   ji, jj
       REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
       !!---------------------------------------------------------------------
@@ -73,23 +75,22 @@
       CALL isfpar_mlt( kt, Kmm, zqhc, zqoce, pqfwf  )
       !
-      ! compute heat and water flux (from isf to oce)
-      pqfwf(:,:) = pqfwf(:,:) * mskisf_par(:,:)
-      zqoce(:,:) = zqoce(:,:) * mskisf_par(:,:)
-      zqhc (:,:) = zqhc(:,:)  * mskisf_par(:,:)
-      !
-      ! compute latent heat flux (from isf to oce)
-      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
-      !
-      ! total heat flux (from isf to oce)
-      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
-      !
-      ! lbclnk on melt and heat fluxes
-      CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ! compute heat and water flux (from isf to oce)
+         pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_par(ji,jj)
+         zqoce(ji,jj) = zqoce(ji,jj) * mskisf_par(ji,jj)
+         zqhc (ji,jj) = zqhc(ji,jj)  * mskisf_par(ji,jj)
+         !
+         ! compute latent heat flux (from isf to oce)
+         zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf    ! 2d latent heat flux (W/m2)
+         !
+         ! total heat flux (from isf to oce)
+         zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) )
+         !
+         ! set temperature content
+         ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
+      END_2D
       !
       ! output fluxes
       CALL isf_diags_flx( Kmm, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc)
-      !
-      ! set temperature content
-      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
       !
       ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
