Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8637)
@@ -16,7 +16,8 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
-   USE ice             ! LIM_3 ice variables
-   USE icevar
-   USE icectl
+   USE ice             ! sea-ice: variables
+   USE icevar          ! sea-ice: operations
+   USE iceitd          ! sea-ice: rebining
+   USE icectl          ! sea-ice: control prints
    USE phycst          ! physical constant
    USE eosbn2          ! equation of state
@@ -39,5 +40,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id: bdyice.F90 8306 2017-07-10 10:18:03Z clem $
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -249,4 +250,8 @@
       END DO !jl
       !
+      ! --- In case categories are out of bounds, do a remapping --- !
+      !     i.e. inputs have not the same ice thickness distribution 
+      !          (set by rn_himean) than the regional simulation
+      IF( jpl > 1 )   CALL ice_itd_reb( kt )
       !      
       IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs')
@@ -273,5 +278,5 @@
       !!------------------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_dyn')
+      IF( ln_timing )   CALL timing_start('bdy_ice_dyn')
       !
       DO ib_bdy=1, nb_bdy
@@ -350,5 +355,5 @@
       END DO
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_dyn')
+      IF( ln_timing )   CALL timing_stop('bdy_ice_dyn')
       !
     END SUBROUTINE bdy_ice_dyn
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/albedooce.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/albedooce.F90	(revision 8586)
+++ 	(revision )
@@ -1,49 +1,0 @@
-MODULE albedooce
-   !!======================================================================
-   !!                       ***  MODULE  albedooce  ***
-   !! Ocean forcing:  bulk thermohaline forcing of the ocean
-   !!=====================================================================
-   !! History :
-   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ocean and ice albedos
-   !!----------------------------------------------------------------------
-   !!   albedo_oce    : albedo for ocean (clear and overcast skies)
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constants
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   albedo_oce   ! routine called by sbccpl
-  
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
-   !! $Id: albedooce.F90 8268 2017-07-03 15:01:04Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE albedo_oce( pa_oce_os , pa_oce_cs )
-      !!----------------------------------------------------------------------
-      !!               ***  ROUTINE albedo_oce  ***
-      !! 
-      !! ** Purpose :   Computation of the albedo of the ocean
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_os   !  albedo of ocean under overcast sky
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
-      !!
-      REAL(wp) ::   zcoef 
-      REAL(wp) ::   rmue = 0.40    !  cosine of local solar altitude
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
-      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
-      !
-   END SUBROUTINE albedo_oce
-
-   !!======================================================================
-END MODULE albedooce
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90	(revision 8637)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90	(revision 8637)
@@ -0,0 +1,49 @@
+MODULE ocealb
+   !!======================================================================
+   !!                       ***  MODULE  ocealb  ***
+   !! Ocean forcing:  bulk thermohaline forcing of the ocean
+   !!=====================================================================
+   !! History :
+   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ocean and ice albedos
+   !!----------------------------------------------------------------------
+   !!   oce_alb    : albedo for ocean (clear and overcast skies)
+   !!----------------------------------------------------------------------
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   oce_alb   ! routine called by sbccpl
+  
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
+   !! $Id: ocealb.F90 8268 2017-07-03 15:01:04Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE oce_alb( palb_os , palb_cs )
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE oce_alb  ***
+      !! 
+      !! ** Purpose :   Computation of the albedo of the ocean
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   palb_os   !  albedo of ocean under overcast sky
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   palb_cs   !  albedo of ocean under clear sky
+      !!
+      REAL(wp) ::   zcoef 
+      REAL(wp) ::   rmue = 0.40    !  cosine of local solar altitude
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      palb_cs(:,:) = zcoef 
+      palb_os(:,:) = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      !
+   END SUBROUTINE oce_alb
+
+   !!======================================================================
+END MODULE ocealb
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8637)
@@ -40,5 +40,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_lim3
-   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
    USE icethd_dh      ! for CALL ice_thd_snwblow
 #endif
@@ -92,5 +92,5 @@
    REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
    REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
-   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice
+   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
    REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
    !
@@ -108,7 +108,14 @@
    REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
 !!gm ref namelist initialize it so remove the setting to false below
-   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012)
+   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
+   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
    !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem)
+!!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
 
    INTEGER  ::   nblk           ! choice of the bulk algorithm
@@ -132,5 +139,7 @@
       !!             ***  ROUTINE sbc_blk_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc )
+!!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc )
+      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
+         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc )
@@ -164,5 +173,5 @@
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
-         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12
+         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
       !!---------------------------------------------------------------------
       !
@@ -255,4 +264,6 @@
          WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
          WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
+         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
+         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
          !
          WRITE(numout,*)
@@ -361,10 +372,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
       REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCd               ! transfer coefficient for momentum      (tau)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCh               ! transfer coefficient for sensible heat (Q_sens)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCe               ! tansfert coefficient for evaporation   (Q_lat)
       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
-      REAL(wp), DIMENSION(jpi,jpj) ::   zt_zu             ! air temperature at wind speed height
-      REAL(wp), DIMENSION(jpi,jpj) ::   zq_zu             ! air spec. hum.  at wind speed height
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
@@ -418,6 +424,4 @@
       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
 
-
-
       ! ----------------------------------------------------------------------------- !
       !     II    Turbulent FLUXES                                                    !
@@ -435,11 +439,11 @@
       !
       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
-         &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
-         &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
-         &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
-         &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
       CASE DEFAULT
          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
@@ -448,14 +452,15 @@
       !                          ! Compute true air density :
       IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
-         zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) )
+         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
       ELSE                                      ! At zt:
          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       END IF
 
-      Cd_oce(:,:) = zCd(:,:)     ! record value of pure ocean-atm. drag (clem)
+!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
+!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
 
       DO jj = 1, jpj             ! tau module, i and j component
          DO ji = 1, jpi
-            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * zCd(ji,jj)   ! using bulk wind speed
+            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
             taum  (ji,jj) = zztmp * wndm  (ji,jj)
             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
@@ -492,11 +497,11 @@
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) )   ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - ztpot(:,:)               )   ! Sensible Heat, using bulk wind speed
+         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
+         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
       ELSE
          !! q_air and t_air are not given at 10m (wind reference height)
          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed
+         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
+         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
       ENDIF
 
@@ -505,6 +510,6 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=zCe , clinfo2=' Ce  : ' )
-         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zCh , clinfo2=' Ch  : ' )
+         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
+         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
          CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
          CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
@@ -574,31 +579,21 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f      ! relative wind module and components at F-point
-      REAL(wp) ::   zwndi_t , zwndj_t                ! relative wind components at T-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCd, zrhoa   ! transfer coefficient for momentum      (tau)
+      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
+      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau)
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('blk_ice_tau')
       !
-      IF( ln_Cd_L12 ) THEN
-         CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012)
-      ELSE
-         zCd(:,:) = Cd_ice             ! constant air-ice drag
-      ENDIF
-
-      ! local scalars ( place there for vector optimisation purposes)
-      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
-      !!
-      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
-
-      !!gm brutal....
-      utau_ice  (:,:) = 0._wp
-      vtau_ice  (:,:) = 0._wp
-      wndm_ice  (:,:) = 0._wp
-      !!gm end
-
-      ! ----------------------------------------------------------------------------- !
-      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   !
-      ! ----------------------------------------------------------------------------- !
+      ! set transfer coefficients to default sea-ice values
+      Cd_atm(:,:) = Cd_ice
+      Ch_atm(:,:) = Cd_ice
+      Ce_atm(:,:) = Cd_ice
+
+      wndm_ice(:,:) = 0._wp      !!gm brutal....
+
+      ! ------------------------------------------------------------ !
+      !    Wind module relative to the moving ice ( U10m - U_ice )   !
+      ! ------------------------------------------------------------ !
       SELECT CASE( cp_ice_msh )
       CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
@@ -606,13 +601,4 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! B grid : NO vector opt
-               ! ... scalar wind at I-point (fld being at T-point)
-               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
-               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
-               zwnorm_f = zrhoa(ji,jj) * zCd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
-               ! ... ice stress at I-point
-               utau_ice(ji,jj) = zwnorm_f * zwndi_f
-               vtau_ice(ji,jj) = zwnorm_f * zwndj_f
                ! ... scalar wind at T-point (fld being at T-point)
                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
@@ -623,11 +609,9 @@
             END DO
          END DO
-         CALL lbc_lnk( utau_ice, 'I', -1. )
-         CALL lbc_lnk( vtau_ice, 'I', -1. )
          CALL lbc_lnk( wndm_ice, 'T',  1. )
          !
       CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vect. opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
@@ -635,9 +619,55 @@
             END DO
          END DO
+         CALL lbc_lnk( wndm_ice, 'T',  1. )
+         !
+      END SELECT
+
+      ! Make ice-atm. drag dependent on ice concentration
+      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
+         CALL Cdn10_Lupkes2012( Cd_atm )
+         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
+      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
+         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
+      ENDIF
+
+!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
+!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
+
+      ! local scalars ( place there for vector optimisation purposes)
+      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
+      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
+
+      !!gm brutal....
+      utau_ice  (:,:) = 0._wp
+      vtau_ice  (:,:) = 0._wp
+      !!gm end
+
+      ! ------------------------------------------------------------ !
+      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
+      ! ------------------------------------------------------------ !
+      SELECT CASE( cp_ice_msh )
+      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! B grid : NO vector opt
+               ! ... scalar wind at I-point (fld being at T-point)
+               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
+                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
+               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
+                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
+               ! ... ice stress at I-point
+               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
+               utau_ice(ji,jj) = zwnorm_f * zwndi_f
+               vtau_ice(ji,jj) = zwnorm_f * zwndj_f
+            END DO
+         END DO
+         CALL lbc_lnk( utau_ice, 'I', -1. )
+         CALL lbc_lnk( vtau_ice, 'I', -1. )
+         !
+      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
-               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
+               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
-               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
+               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
             END DO
@@ -645,5 +675,4 @@
          CALL lbc_lnk( utau_ice, 'U', -1. )
          CALL lbc_lnk( vtau_ice, 'V', -1. )
-         CALL lbc_lnk( wndm_ice, 'T',  1. )
          !
       END SELECT
@@ -684,17 +713,10 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zCd            ! transfer coefficient for momentum      (tau)
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('blk_ice_flx')
       !
-      IF( ln_Cd_L12 ) THEN
-         CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012)
-      ELSE
-         zCd(:,:) = Cd_ice             ! constant air-ice drag
-      ENDIF
-
-      !
-      ! local scalars ( place there for vector optimisation purposes)
+      !
+      ! local scalars
       zcoef_dqlw   = 4.0 * 0.95 * Stef
       zcoef_dqla   = -Ls * 11637800. * (-5897.8)
@@ -724,13 +746,13 @@
                ! ----------------------------!
 
-               ! ... turbulent heat fluxes
+               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
+               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * zCd(ji,jj) * wndm_ice(ji,jj)   &
-                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)  ) )
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
+                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
                ! Latent heat sensitivity for ice (Dqla/Dt)
                IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
-                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * zCd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) )
+                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl))
                ELSE
                   dqla_ice(ji,jj,jl) = 0._wp
@@ -738,5 +760,5 @@
 
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
-               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj)
+               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
 
                ! ----------------------------!
@@ -935,5 +957,5 @@
 #if defined key_lim3
 
-   SUBROUTINE Cdn10_Lupkes2012( pCd )
+   SUBROUTINE Cdn10_Lupkes2012( Cd )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
@@ -965,5 +987,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd   ! air-ice drag coefficient
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
       REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
       REAL(wp), PARAMETER ::   znu   = 1._wp
@@ -973,17 +995,132 @@
       !!----------------------------------------------------------------------
       zcoef = znu + 1._wp / ( 10._wp * zbeta )
-      !
+
       ! generic drag over a cell partly covered by ice
-      !!pCd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
-      !!   &       Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
-      !!   &       zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
+      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
+      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
+      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
 
       ! ice-atm drag
-      pCd(:,:) = Cd_ice +  &                                                         ! pure ice drag
-         &       zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
-      !
+      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag
+         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
+      
    END SUBROUTINE Cdn10_Lupkes2012
+
+
+   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
+      !!
+      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation
+      !!                 between sea-ice and atmosphere with distinct momentum 
+      !!                 and heat coefficients depending on sea-ice concentration 
+      !!                 and atmospheric stability (no meltponds effect for now).
+      !!                
+      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
+      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
+      !!                 it considers specific skin and form drags (Andreas et al. 2010)
+      !!                 to compute neutral transfert coefficients for both heat and 
+      !!                 momemtum fluxes. Atmospheric stability effect on transfert
+      !!                 coefficient is also taken into account following Louis (1979).
+      !!
+      !! ** References : Lupkes et al. JGR 2015 (theory)
+      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
+      !!
+      !!----------------------------------------------------------------------
+      !
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
+      !
+      ! ECHAM6 constants
+      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
+      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
+      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
+      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
+      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
+      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
+      REAL(wp), PARAMETER ::   zc2          = zc * zc
+      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
+      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
+      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
+      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
+      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
+      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
+      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
+      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
+      REAL(wp) ::   zrib_o, zrib_i
+      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
+      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
+      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
+      REAL(wp) ::   zCdn_form_tmp
+      !!----------------------------------------------------------------------
+
+      ! Momentum Neutral Transfert Coefficients (should be a constant)
+      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
+      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
+      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
+      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
+
+      ! Heat Neutral Transfert Coefficients
+      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details)
+     
+      ! Atmospheric and Surface Variables
+      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin
+      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
+      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
+      !
+      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
+         DO ji = fs_2, fs_jpim1
+            ! Virtual potential temperature [K]
+            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
+            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
+            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
+            
+            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
+            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
+            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
+            
+            ! Momentum and Heat Neutral Transfert Coefficients
+            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
+            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
+                       
+            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
+            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
+            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
+            IF( zrib_o <= 0._wp ) THEN
+               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10
+               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
+                  &             )**zgamma )**z1_gamma
+            ELSE
+               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
+               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
+            ENDIF
+            
+            IF( zrib_i <= 0._wp ) THEN
+               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
+               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
+            ELSE
+               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
+               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
+            ENDIF
+            
+            ! Momentum Transfert Coefficients (Eq. 38)
+            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
+               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
+            
+            ! Heat Transfert Coefficients (Eq. 49)
+            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
+               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
+            !
+         END DO
+      END DO
+      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
+      !
+   END SUBROUTINE Cdn10_Lupkes2015
+
 #endif
-   
 
    !!======================================================================
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90	(revision 8637)
@@ -44,5 +44,4 @@
    USE lib_fortran    ! to use key_nosignedzero
 
-
    IMPLICIT NONE
    PRIVATE
@@ -50,15 +49,15 @@
    PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
 
-   !! COARE own values for given constants:
-   REAL(wp), PARAMETER :: &
-      &   zi0     = 600.,  &   !: scale height of the atmospheric boundary layer...1
-      &  Beta0    = 1.25, &   !: gustiness parameter
-      &  rctv0    = 0.608     !: constant to obtain virtual temperature...
-
-
+   !                                           !! COARE own values for given constants:
+   REAL(wp), PARAMETER ::   zi0     = 600.      ! scale height of the atmospheric boundary layer...1
+   REAL(wp), PARAMETER ::   Beta0   =   1.25    ! gustiness parameter
+   REAL(wp), PARAMETER ::   rctv0   =   0.608   ! constant to obtain virtual temperature...
+
+   !!----------------------------------------------------------------------
 CONTAINS
 
    SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
-      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk )
+      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+      &                   Cdn, Chn, Cen                       )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare  ***
@@ -106,4 +105,5 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
       INTEGER :: j_itt
@@ -246,4 +246,9 @@
       Ce   = ztmp0*q_star/dq_zu
       !
+      ztmp1 = zu + z0
+      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
+      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
+      Cen = Chn
+      !
       CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu )
       CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90	(revision 8637)
@@ -58,6 +58,7 @@
 CONTAINS
 
-   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
-      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk )
+   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       &
+      &                      Cdn, Chn, Cen                        )
       !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare3p5  ***
@@ -105,4 +106,5 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
       INTEGER :: j_itt
@@ -252,4 +254,9 @@
       Ce   = ztmp0*q_star/dq_zu
       !
+      ztmp1 = zu + z0
+      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
+      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
+      Cen = Chn
+      !
       CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu )
       CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90	(revision 8637)
@@ -64,5 +64,6 @@
 
    SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   &
-      &                   Cd, Ch, Ce , t_zu, q_zu, U_blk )
+      &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         &
+      &                   Cdn, Chn, Cen                           )
       !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
@@ -112,4 +113,5 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
       INTEGER :: j_itt
@@ -266,4 +268,5 @@
             dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
             dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
+
          END IF
 
@@ -271,5 +274,5 @@
          ztmp1 = zu + z0
          ztmp0 = ztmp1*Linv
-         func_m = log(ztmp1) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv)
+         func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
          func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
 
@@ -280,4 +283,9 @@
       ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
       Ce = vkarmn*vkarmn/(func_m*ztmp1)
+
+      ztmp1 = zu + z0
+      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
+      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
+      Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q))
 
       CALL wrk_dealloc( jpi,jpj,   u_star, t_star, q_star, func_m, func_h, dt_zu, dq_zu, Linv )
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 8637)
@@ -53,5 +53,6 @@
 
    SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
-      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk )
+      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+      &                  Cdn, Chn, Cen                       )
       !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ncar  ***
@@ -112,4 +113,5 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
       INTEGER ::   j_itt
@@ -199,8 +201,10 @@
             ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))
             ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
+            Cdn(:,:) = ztmp0
             sqrt_Cd_n10 = sqrt(ztmp0)
 
             stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
             Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
+            Chn(:,:) = Cx_n10
 
             !! Update of transfer coefficients:
@@ -216,14 +220,15 @@
 
          Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
+         Cen(:,:) = Cx_n10
          ztmp1 = 1. + Cx_n10*ztmp0
          Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
-
+         !
       END DO
-
+      !
       CALL wrk_dealloc( jpi,jpj,   Cx_n10, sqrt_Cd_n10, zeta_u, stab )
       CALL wrk_dealloc( jpi,jpj,   zpsi_h_u, ztmp0, ztmp1, ztmp2 )
-
+      !
       IF( nn_timing == 1 )   CALL timing_stop('turb_ncar')
-
+      !
    END SUBROUTINE turb_ncar
 
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8637)
@@ -33,5 +33,5 @@
    USE geo2ocean      ! 
    USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
-   USE albedooce      ! 
+   USE ocealb         ! 
    USE eosbn2         ! 
    USE sbcrnf, ONLY : l_rnfcpl
@@ -173,10 +173,10 @@
                                          !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
    TYPE ::   DYNARR     
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3   
    END TYPE DYNARR
 
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -202,5 +202,5 @@
       ierr(:) = 0
       !
-      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
+      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
       
 #if ! defined key_lim3 && ! defined key_cice
@@ -736,7 +736,7 @@
       !     2. receiving mixed oce-ice solar radiation 
       IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
-         CALL albedo_oce( zaos, zacs )
+         CALL oce_alb( zaos, zacs )
          ! Due to lack of information on nebulosity : mean clear/overcast sky
-         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
+         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
       ENDIF
 
@@ -1885,6 +1885,6 @@
 !       ( see OASIS3 user guide, 5th edition, p39 )
          zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( albedo_oce_mix(:,:  ) * ziceld(:,:)       &
-            &                     + palbi         (:,:,1) * picefr(:,:) ) )
+            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
+            &                     + palbi      (:,:,1) * picefr(:,:) ) )
       END SELECT
       IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
@@ -2052,5 +2052,5 @@
                    ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
                 ELSEWHERE
-                   ztmp1(:,:) = albedo_oce_mix(:,:)
+                   ztmp1(:,:) = alb_oce_mix(:,:)
                 END WHERE
              CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
@@ -2080,8 +2080,8 @@
 
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
-         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
+         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
          DO jl=1,jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
-         ENDDO
+         END DO
          CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
       ENDIF
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 8637)
@@ -236,5 +236,5 @@
       ENDIF
       !
-      IF( .NOT. l_ssm_mean ) THEN   ! default initialisation. needed by lim_istate
+      IF( .NOT.l_ssm_mean ) THEN   ! default initialisation. needed by iceistate
          !
          IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays'
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 8637)
@@ -33,24 +33,19 @@
    PRIVATE
 
-   PUBLIC   zdf_iwm         ! called in step module 
-   PUBLIC   zdf_iwm_init    ! called in nemogcm module 
-
-   !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing *
-   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
-   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
-   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
-
-   REAL(wp) ::  r1_6 = 1._wp / 6._wp
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_iwm     ! power available from high-mode wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_iwm     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_iwm     ! power available from low-mode, critical slope wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_iwm     ! WKB decay scale for high-mode energy dissipation (m)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_iwm     ! decay scale for low-mode critical slope dissipation (m)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_iwm     ! local energy density available for mixing (W/kg)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_iwm     ! buoyancy flux Kz * N^2 (W/kg)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_iwm    ! vertically integrated buoyancy flux (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity
+   PUBLIC   zdf_iwm        ! called in step module 
+   PUBLIC   zdf_iwm_init   ! called in nemogcm module 
+
+   !                      !!* Namelist  namzdf_iwm : internal wave-driven mixing *
+   INTEGER ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
+   LOGICAL ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
+   LOGICAL ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
+
+   REAL(wp)::  r1_6 = 1._wp / 6._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ebot_iwm   ! power available from high-mode wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   epyc_iwm   ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ecri_iwm   ! power available from low-mode, critical slope wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbot_iwm   ! WKB decay scale for high-mode energy dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hcri_iwm   ! decay scale for low-mode critical slope dissipation (m)
 
    !! * Substitutions
@@ -67,8 +62,6 @@
       !!                ***  FUNCTION zdf_iwm_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE(     ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj)    ,   &
-      &             hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj),  emix_iwm(jpi,jpj,jpk),   &
-      &         bflx_iwm(jpi,jpj,jpk), pcmap_iwm(jpi,jpj), zav_ratio(jpi,jpj,jpk),   & 
-      &         zav_wave(jpi,jpj,jpk), STAT=zdf_iwm_alloc     )
+      ALLOCATE( ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj) ,     &
+      &         hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj)                     , STAT=zdf_iwm_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_iwm_alloc )
@@ -85,20 +78,20 @@
       !!
       !! ** Method  : - internal wave-driven vertical mixing is given by:
-      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_iwm /( Nu * N^2 )  )
-      !!              where emix_iwm is the 3D space distribution of the wave-breaking 
+      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = zemx_iwm /( Nu * N^2 )  )
+      !!              where zemx_iwm is the 3D space distribution of the wave-breaking 
       !!              energy and Nu the molecular kinematic viscosity.
       !!              The function f(Reb) is linear (constant mixing efficiency)
       !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
       !!
-      !!              - Compute emix_iwm, the 3D power density that allows to compute
+      !!              - Compute zemx_iwm, the 3D power density that allows to compute
       !!              Reb and therefrom the wave-induced vertical diffusivity.
       !!              This is divided into three components:
       !!                 1. Bottom-intensified low-mode dissipation at critical slopes
-      !!                     emix_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm )
+      !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm )
       !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm
       !!              where hcri_iwm is the characteristic length scale of the bottom 
       !!              intensification, ecri_iwm a map of available power, and H the ocean depth.
       !!                 2. Pycnocline-intensified low-mode dissipation
-      !!                     emix_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
+      !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
       !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
       !!              where epyc_iwm is a map of available power, and nn_zpyc
@@ -106,5 +99,5 @@
       !!              energy dissipation.
       !!                 3. WKB-height dependent high mode dissipation
-      !!                     emix_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
+      !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
       !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) )
       !!              where hbot_iwm is the characteristic length scale of the WKB bottom 
@@ -121,6 +114,5 @@
       !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
       !!
-      !! ** Action  : - Define emix_iwm used to compute internal wave-induced mixing
-      !!              - avt, avs, avm, increased by internal wave-driven mixing    
+      !! ** Action  : - avt, avs, avm, increased by tide internal wave-driven mixing    
       !!
       !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
@@ -132,20 +124,30 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zztmp        ! scalar workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zfact     ! Used for vertical structure
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zhdep     ! Ocean depth
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwkb      ! WKB-stretched height above bottom
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zweight   ! Weight for high mode vertical distribution
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_t     ! Molecular kinematic viscosity (T grid)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_w     ! Molecular kinematic viscosity (W grid)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zReb      ! Turbulence intensity parameter
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwkb        ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zweight     ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_t       ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_w       ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zReb        ! Turbulence intensity parameter
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zemx_iwm    ! local energy density available for mixing (W/kg)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_ratio   ! S/T diffusivity ratio (only for ln_tsdiff=T)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_wave    ! Internal wave-induced diffusivity
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3d  ! 3D workspace used for iom_put 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2d  ! 2D     -      -    -     -
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('zdf_iwm')
       !
-      !                      ! ----------------------------- !
-      !                      !  Internal wave-driven mixing  !  (compute zav_wave)
-      !                      ! ----------------------------- !
+      !                       !* Set to zero the 1st and last vertical levels of appropriate variables
+      zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp
+      zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp
+      zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp
+      !
+      !                       ! ----------------------------- !
+      !                       !  Internal wave-driven mixing  !  (compute zav_wave)
+      !                       ! ----------------------------- !
       !                             
-      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
+      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
       DO jj = 1, jpj                ! part independent of the level
@@ -156,7 +158,7 @@
          END DO
       END DO
-
+!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept_n - sshn
       DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
+         zemx_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
             &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
             &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
@@ -186,5 +188,5 @@
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
          END DO
          !
@@ -203,5 +205,5 @@
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
          END DO
          !
@@ -253,5 +255,5 @@
       !
       DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
+         zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
             &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
 !!gm  use of e3t_n just above?
@@ -269,5 +271,5 @@
       ! Calculate turbulence intensity parameter Reb
       DO jk = 2, jpkm1
-         zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
+         zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
       END DO
       !
@@ -349,18 +351,20 @@
       !                             !* output internal wave-driven mixing coefficient
       CALL iom_put( "av_wave", zav_wave )
-                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_iwm), 
-                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm)
+                                    !* output useful diagnostics: Kz*N^2 , 
+!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5)
+                                    !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm)
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
-         bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
-         pcmap_iwm(:,:) = 0._wp
+         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) )
+         z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
+         z2d(:,:) = 0._wp
          DO jk = 2, jpkm1
-            pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk)
-         END DO
-         pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:)
-         CALL iom_put( "bflx_iwm", bflx_iwm )
-         CALL iom_put( "pcmap_iwm", pcmap_iwm )
-      ENDIF
-      CALL iom_put( "bn2", rn2 )
-      CALL iom_put( "emix_iwm", emix_iwm )
+            z2d(:,:) = z2d(:,:) + e3w_n(:,:,jk) * z3d(:,:,jk) * wmask(:,:,jk)
+         END DO
+         z2d(:,:) = rau0 * z2d(:,:)
+         CALL iom_put( "bflx_iwm", z3d )
+         CALL iom_put( "pcmap_iwm", z2d )
+         DEALLOCATE( z2d , z3d )
+      ENDIF
+      CALL iom_put( "emix_iwm", zemx_iwm )
       
       IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
@@ -466,12 +470,4 @@
       ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:)
 
-      ! Set once for all to zero the first and last vertical levels of appropriate variables
-      emix_iwm (:,:, 1 ) = 0._wp
-      emix_iwm (:,:,jpk) = 0._wp
-      zav_ratio(:,:, 1 ) = 0._wp
-      zav_ratio(:,:,jpk) = 0._wp
-      zav_wave (:,:, 1 ) = 0._wp
-      zav_wave (:,:,jpk) = 0._wp
-
       zbot = glob_sum( e1e2t(:,:) * ebot_iwm(:,:) )
       zpyc = glob_sum( e1e2t(:,:) * epyc_iwm(:,:) )
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8637)
@@ -62,7 +62,7 @@
       !!                     ***  ROUTINE stp  ***
       !!
-      !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
-      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
-      !!              - Tme stepping  of TRC (passive tracer eqs.)
+      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
+      !!              - Time stepping of ESIM (dynamic and thermodynamic eqs.)
+      !!              - Time stepping of TRC  (passive tracer eqs.)
       !!
       !! ** Method  : -1- Update forcings and data
