Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfcavgam.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfcavgam.F90	(revision 14353)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfcavgam.F90	(revision 15574)
@@ -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
@@ -30,4 +30,6 @@
    PUBLIC   isfcav_gammats
 
+   !! * Substitutions   
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
@@ -42,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
@@ -55,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
@@ -92,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)')
@@ -133,13 +136,16 @@
       REAL(wp),                     INTENT(in   ) :: pke2         ! background velocity
       !!---------------------------------------------------------------------
+      INTEGER  :: ji, jj                     ! loop index
       REAL(wp), DIMENSION(jpi,jpj) :: zustar
       !!---------------------------------------------------------------------
       !
-      ! compute ustar (AD15 eq. 27)
-      zustar(:,:) = SQRT( pCd(:,:) * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) ) * mskisf_cav(:,:)
-      !
-      ! Compute gammats
-      pgt(:,:) = zustar(:,:) * rn_gammat0
-      pgs(:,:) = zustar(:,:) * rn_gammas0
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ! compute ustar (AD15 eq. 27)
+         zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj)
+         !
+         ! Compute gammats
+         pgt(ji,jj) = zustar(ji,jj) * rn_gammat0
+         pgs(ji,jj) = zustar(ji,jj) * rn_gammas0
+      END_2D
       !
       ! output ustar
@@ -148,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 
@@ -166,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
@@ -185,10 +192,4 @@
       !!---------------------------------------------------------------------
       !
-      ! compute ustar
-      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )
-      !
-      ! output ustar
-      CALL iom_put('isfustar',zustar(:,:))
-      !
       ! compute Pr and Sc number (eq ??)
       zPr =   13.8_wp
@@ -200,57 +201,51 @@
       !
       ! compute gamma
-      DO ji = 2, jpi
-         DO jj = 2, jpj
-            ikt = mikt(ji,jj)
-
-            IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
-               pgt = rn_gammat0
-               pgs = 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)
-               zts(jp_sal) = pstbl(ji,jj)
-               zdep        = gdepw(ji,jj,ikt,Kmm)
-               !
-               CALL eos_rab( zts, zdep, zab, Kmm )
-               !
-               ! compute length scale (Eq ??)
-               zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
-               !
-               ! compute Monin Obukov Length
-               ! Maximum boundary layer depth (Eq ??)
-               zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
-               !
-               ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
-               zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
-               !
-               ! 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)))
-               !
-               ! compute the sublayer thickness (Eq ??)
-               zhnu = 5 * znu / zustar(ji,jj)
-               !
-               ! compute gamma turb (Eq ??)
-               zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
+      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
+            pgt(ji,jj) = rn_gammat0
+            pgs(ji,jj) = rn_gammas0
+         ELSE
+            ! compute bouyancy 
+            zts(jp_tem) = pttbl(ji,jj)
+            zts(jp_sal) = pstbl(ji,jj)
+            zdep        = gdepw(ji,jj,ikt,Kmm)
+            !
+            CALL eos_rab( zts, zdep, zab, Kmm )
+            !
+            ! compute length scale (Eq ??)
+            zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
+            !
+            ! compute Monin Obukov Length
+            ! Maximum boundary layer depth (Eq ??)
+            zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
+            !
+            ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
+            zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+            zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+            !
+            ! compute eta* (stability parameter) (Eq ??)
+            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 / MAX( 1.e-20, zustar(ji,jj) )
+            !
+            ! compute gamma turb (Eq ??)
+            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
-               !
-               ! compute gammats
-               pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
-               pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
-            END IF
-         END DO
-      END DO
+            !
+            ! compute gammats
+            pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+            pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
+         END IF
+      END_2D
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
 
    END SUBROUTINE gammats_vel_stab
