Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11359)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11360)
@@ -335,9 +335,8 @@
   
 #if defined key_si3
-      CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in
-         &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in
-         &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice           ,   &   !   <<= in
-         &            pseni=zseni, pevpi=zevpi, ptsui=tm_su             ,   &   !   <<= out
-         &            pssqi=zssqi, pcd_dui=zcd_dui                      )       !   =>> out
+      CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
+         &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),    &   !   <<= in
+         &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
+         &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
 #endif  
 
@@ -364,8 +363,8 @@
       !!-------------------------------------------------------------------------------------------   
 
-      CALL blk_oce_2( kt, tq_abl(:,:,2,nt_n,jp_ta),                            &
-         &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
-         &                sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
-         &                sst_m, zsen, zevp                                )
+      CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            &
+         &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
+         &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
+         &            sst_m, zsen, zevp                                )
 
 #if defined key_si3
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90	(revision 11359)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90	(revision 11360)
@@ -72,7 +72,7 @@
          CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
          CASE( jp_blk     )   ;    CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &
-            &                                      sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1),   &
-            &                                      u_ice, v_ice,   &    ! inputs
-            &                                      putaui = utau_ice, pvtaui = vtau_ice )    ! outputs                             
+            &                                      sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &
+            &                                      sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su    ,   &   ! inputs
+            &                                      putaui = utau_ice, pvtaui = vtau_ice            )       ! outputs                             
  !        CASE( jp_abl     )    utau_ice & vtau_ice are computed in ablmod
          CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11359)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11360)
@@ -75,9 +75,10 @@
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
    INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
-   INTEGER ::   ndim_A, ndim_hA   ! grid_T file   
+   INTEGER ::   ndim_A, ndim_hA                      ! ABL file   
    INTEGER ::   nid_A, nz_A, nh_A                    ! grid_ABL file   
    INTEGER ::   ndex(1)                              ! ???
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV, ndex_hA
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V, ndex_A
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
 
@@ -445,5 +446,6 @@
       INTEGER  ::   ji, jj, jk                               ! dummy loop indices
       INTEGER  ::   ierr                                     ! error code return from allocation
-      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma, ipka   ! local integers
+      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      INTEGER  ::   ipka                                     ! ABL
       INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
@@ -451,5 +453,5 @@
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
       ! 
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11359)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11360)
@@ -112,5 +112,5 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds
-      REAL(wp), DIMENSION(2   ,jpkam1)      ::   za_bnds     ! ABL vertical boundaries
+      REAL(wp), DIMENSION(2,jpkam1)         :: za_bnds   ! ABL vertical boundaries
       LOGICAL ::   ll_tmppatch = .TRUE.    !: seb: patch before we remove periodicity
       INTEGER ::   nldi_save, nlei_save    !:      and close boundaries in output files
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11359)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11360)
@@ -48,5 +48,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_si3
-   USE ice     , ONLY :   jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
+   USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif
    USE icethd_dh      ! for CALL ice_thd_snwblow
 #endif
@@ -77,7 +77,7 @@
 #endif 
 
-   INTERFACE cp_air
-      MODULE PROCEDURE cp_air_0d, cp_air_2d
-   END INTERFACE
+  INTERFACE cp_air
+    MODULE PROCEDURE cp_air_0d, cp_air_2d
+  END INTERFACE
 
    !                                   !!* Namelist namsbc_blk : bulk parameters
@@ -106,18 +106,18 @@
    REAL(wp)        , PARAMETER ::   Cp_dry = 1005.0        !  Specic heat of dry air, constant pressure      [J/K/kg]
    REAL(wp)        , PARAMETER ::   Cp_vap = 1860.0        !  Specic heat of water vapor, constant pressure  [J/K/kg]
-   REAL(wp), PUBLIC, PARAMETER ::   R_dry  = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
+   REAL(wp), PUBLIC, PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
    REAL(wp)        , PARAMETER ::   R_vap  = 461.495_wp    !  Specific gas constant for water vapor          [J/K/kg]
    REAL(wp)        , PARAMETER ::   reps0  = R_dry/R_vap   !  ratio of gas constant for dry air and water vapor => ~ 0.622
-   REAL(wp), PUBLIC, PARAMETER ::   rctv0  = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   REAL(wp), PUBLIC, PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
    !                                                      !!! Bulk parameters
    REAL(wp)        , PARAMETER ::   cpa    = 1000.5        ! specific heat of air (only used for ice fluxes now...)
    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     ! transfer coefficient over ice
+   REAL(wp)        , PARAMETER ::   rcdice =    1.4e-3     ! transfer coefficient over ice
    REAL(wp)        , PARAMETER ::   albo   =    0.066      ! ocean albedo assumed to be constant
    !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm , Ch_atm , Ce_atm   ! transfer coeffs. for momentum, sensible heat, and evaporation
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral  coeffs (needed for Lupkes et al. 2015 bulk scheme)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (needed for Lupkes 2015 bulk scheme)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
 
    INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
@@ -149,6 +149,7 @@
       !!             ***  ROUTINE 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 )
+      ALLOCATE( t_zu(jpi,jpj), q_zu(jpi,jpj),                                            &
+         &      Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj),                    &
+         &      Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc )
       !
       CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
@@ -268,7 +269,12 @@
          rn_zu  = ght_abl(2)
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '  ABL formulation:   overwrite rn_zqt & rn_zu with ABL first level altitude'
-      ENDIF
-      !           
+         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
+      ENDIF
+      !
+      ! set transfer coefficients to default sea-ice values
+      Cd_ice(:,:) = rcdice
+      Ch_ice(:,:) = rcdice
+      Ce_ice(:,:) = rcdice
+      !
       IF(lwp) THEN                     !** Control print
          !
@@ -349,5 +355,5 @@
             &                zssq, zcd_du, zsen, zevp )                              !   =>> out
          
-         CALL blk_oce_2( kt, sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in
+         CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in
             &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in
             &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in
@@ -410,4 +416,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean
+      REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean
+      REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean
       !!---------------------------------------------------------------------
       !
@@ -461,11 +470,11 @@
       !
       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! NCAR-COREv2
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
+         &                                           zcd_oce, zch_oce, zce_oce, 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, pssq, phumi, wndm,   &  ! COARE v3.0
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
+         &                                           zcd_oce, zch_oce, zce_oce, 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, pssq, phumi, wndm,   &  ! COARE v3.5
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
+         &                                           zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! ECMWF
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
+         &                                           zcd_oce, zch_oce, zce_oce, 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' )
@@ -479,7 +488,7 @@
                zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1)
                wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
-               pcd_du(ji,jj) = zztmp * Cd_atm(ji,jj)
-               psen(ji,jj)   = zztmp * Ch_atm(ji,jj)
-               pevp(ji,jj)   = zztmp * Ce_atm(ji,jj)       
+               pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
+               psen(ji,jj)   = zztmp * zch_oce(ji,jj)
+               pevp(ji,jj)   = zztmp * zce_oce(ji,jj)       
             END DO
          END DO
@@ -498,7 +507,7 @@
                zztmp         = zrhoa(ji,jj) * zU_zu(ji,jj)
 !!gm to be done by someone: check the efficiency of the call of cp_air within do loops 
-               psen  (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * Ch_atm(ji,jj) * (  zst(ji,jj) - t_zu(ji,jj) )
-               pevp  (ji,jj) = rn_efac*MAX( 0._wp,     zztmp * Ce_atm(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) )
-               zztmp         = zztmp * Cd_atm(ji,jj)
+               psen  (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * zch_oce(ji,jj) * (  zst(ji,jj) - t_zu(ji,jj) )
+               pevp  (ji,jj) = rn_efac*MAX( 0._wp,     zztmp * zce_oce(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) )
+               zztmp         = zztmp * zcd_oce(ji,jj)
                pcd_du(ji,jj) = zztmp
                taum  (ji,jj) = zztmp * wndm  (ji,jj) 
@@ -540,6 +549,6 @@
 
 
-   SUBROUTINE blk_oce_2( kt, ptair, pqsr, pqlw, pprec,   &   ! <<= in
-              &              psnow, pst , psen, pevp     )   ! <<= in
+   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in
+              &          psnow, pst , psen, pevp     )   ! <<= in
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_oce_2  ***
@@ -557,5 +566,4 @@
       !!              - emp     : evaporation minus precipitation       (kg/m2/s)
       !!---------------------------------------------------------------------
-      INTEGER , INTENT(in)                 ::   kt    ! time step index
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr
@@ -639,6 +647,6 @@
          CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
          CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ')
-      ENDIF 
-	  !
+      ENDIF
+      !
    END SUBROUTINE blk_oce_2
 
@@ -789,6 +797,6 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice,       &    ! inputs
-      &                  putaui, pvtaui, pseni, pevpi, ptsui, pssqi, pcd_dui   )    ! outputs
+   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui,  &   ! inputs
+      &                  putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui             )   ! optional outputs
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_ice_1  ***
@@ -807,22 +815,19 @@
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   puice   ! sea-ice velocity on I or C grid [m/s]
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pvice   ! "
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptsui   ! sea-ice surface temperature [K]
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   putaui  ! if ln_blk
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pvtaui  ! if ln_blk
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pseni   ! if ln_abl
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pevpi   ! if ln_abl
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   ptsui   ! if ln_abl 
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pssqi   ! if ln_abl
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pcd_dui ! if ln_abl 
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zwndi_t , zwndj_t, zootm_su, zztmp1, zztmp2 ! relative wind components at T-point
+      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
+      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature
+      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui     ! transfer coefficient for momentum      (tau)
-      !!---------------------------------------------------------------------
-      !
-      ! set transfer coefficients to default sea-ice values
-      Cd_atm(:,:) = Cd_ice
-      Ch_atm(:,:) = Cd_ice
-      Ce_atm(:,:) = Cd_ice
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui   ! transfer coefficient for momentum      (tau)
+      !!---------------------------------------------------------------------
 
       ! ------------------------------------------------------------ !
@@ -841,17 +846,19 @@
       ! 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
+         CALL Cdn10_Lupkes2012( Cd_ice )
+         Ch_ice(:,:) = Cd_ice(:,:)       ! momentum and heat transfer coef. are considered identical
+         Ce_ice(:,:) = Cd_ice(:,:)
       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.
+         CALL Cdn10_Lupkes2015( ptsui, Cd_ice, Ch_ice ) 
+         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical
+      ENDIF
+
+      IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef.
+      IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! 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( ptair(:,:), phumi(:,:), pslp(:,:) )
-      zcd_dui(:,:) = wndm_ice(:,:) * Cd_atm(:,:)
+      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
       
       IF( ln_blk ) THEN 
@@ -876,12 +883,12 @@
       ELSE
          zztmp1 = 11637800.0_wp
-		 zztmp2 =    -5897.8_wp
-		 DO jj = 1, jpj
+	 zztmp2 =    -5897.8_wp
+         DO jj = 1, jpj
             DO ji = 1, jpi
                pcd_dui(ji,jj) = zcd_dui (ji,jj)
-               pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_atm(ji,jj)
-               pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_atm(ji,jj)
-			   zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??)
-			   pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj)
+               pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
+               pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
+               zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??)
+               pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj)
             END DO
          END DO
@@ -905,5 +912,5 @@
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature [K]
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
@@ -959,14 +966,14 @@
                ! ----------------------------!
 
-               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_1
+               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))
+               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))
                ! Latent Heat
                zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) )
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ce_atm(ji,jj) * wndm_ice(ji,jj) *  &
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  &
                   &                ( 11637800. * zztmp2 / zrhoa(ji,jj) - phumi(ji,jj) ) )
                ! 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 * Ce_atm(ji,jj) * wndm_ice(ji,jj) *  &
+                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  &
                      &                 z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2
                ELSE
@@ -975,5 +982,5 @@
 
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
-               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
+               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj)
                
                ! ----------------------------!
@@ -1151,5 +1158,5 @@
    
 
-   SUBROUTINE Cdn10_Lupkes2012( Cd )
+   SUBROUTINE Cdn10_Lupkes2012( pcd )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
@@ -1181,5 +1188,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd
       REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
       REAL(wp), PARAMETER ::   znu   = 1._wp
@@ -1196,11 +1203,11 @@
 
       ! ice-atm drag
-      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
+      pcd(:,:) = rcdice +  &                                                         ! 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 )
+   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pcd, pch )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
@@ -1223,7 +1230,8 @@
       !!----------------------------------------------------------------------
       !
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
-      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ptm_su ! sea-ice surface temperature [K]
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd    ! momentum transfert coefficient
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pch    ! heat transfert coefficient
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
       !
       ! ECHAM6 constants
@@ -1252,30 +1260,24 @@
       REAL(wp) ::   zCdn_form_tmp
       !!---------------------------------------------------------------------------
-      ! fl: global variable tm_su could be used directly instead of recomputing it 
-	  !
-      ! mean temperature
-      WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
-      ELSEWHERE                       ;   ztm_su(:,:) = rt0
-      ENDWHERE
-      
+      !
       ! 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      = zCdn_skin_ice   ! Eq. 7
       !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)
+      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
      
       ! 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( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
-      !
-      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
+      zqi_sat(:,:) = 0.98_wp * q_sat( ptm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
+      !
+      DO jj = 2, jpjm1           ! reduced loop is necessary for reproductibility
          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 = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
+            zthetav_is = ptm_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
             
@@ -1290,5 +1292,5 @@
             ! 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)
+            z0i = z0_skin_ice                                             ! over ice
             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
@@ -1309,14 +1311,15 @@
             
             ! 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) )
+            pcd(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) )
+            pch(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( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
+      !
+      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1., pch, 'T', 1. )
       !
    END SUBROUTINE Cdn10_Lupkes2015
