Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90	(revision 8586)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_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
-   
 
    !!======================================================================
