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 11275)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90	(revision 11304)
@@ -71,5 +71,8 @@
       SELECT CASE( ksbc )
          CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
-         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk         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                              ! Bulk         formulation
          CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation
       END SELECT
@@ -140,5 +143,6 @@
                                   CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
       CASE( jp_blk )              !--- bulk formulation
-                                  CALL blk_ice_flx    ( t_su, h_s, h_i, alb_ice )    ! 
+                                  CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),    &
+            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )    ! 
          IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
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 11275)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11304)
@@ -29,5 +29,5 @@
    USE lib_mpp           ! MPP library
 #if defined key_iomput
-   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3w_abl, jpka, jpkam1
+   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1
    USE trc_oce  , ONLY :   nn_dttrc                  ! frequency of step on passive tracers
    USE icb_oce  , ONLY :   nclasses, class_num       ! iceberg classes
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 11275)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11304)
@@ -31,6 +31,6 @@
    !!   L_vap         : latent heat of vaporization of water as a function of temperature
    !!             sea-ice case only : 
-   !!   blk_ice_tau   : provide the air-ice stress
-   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_1   : provide the air-ice stress
+   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
    !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
    !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
@@ -48,5 +48,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_si3
-   USE ice     , ONLY :   u_ice, v_ice, 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, t_su, rn_cnd_s, hfx_err_dif
    USE icethd_dh      ! for CALL ice_thd_snwblow
 #endif
@@ -72,6 +72,6 @@
    PUBLIC   cp_air        ! called in ablmod
 #if defined key_si3
-   PUBLIC   blk_ice_tau   ! routine called in icesbc
-   PUBLIC   blk_ice_flx   ! routine called in icesbc
+   PUBLIC   blk_ice_1     ! routine called in icesbc
+   PUBLIC   blk_ice_2     ! routine called in icesbc
    PUBLIC   blk_ice_qcn   ! routine called in icesbc
 #endif 
@@ -83,10 +83,10 @@
   !!Lolo: should ultimately be moved in the module with all physical constants ?
 !!gm  : In principle, yes.
-   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), 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), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   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)        , 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
 
    INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
@@ -118,7 +118,7 @@
    LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
    !
-   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
-   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation
-   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress
+   REAL(wp)         ::   rn_pfac        ! multiplication factor for precipitation
+   REAL(wp), PUBLIC ::   rn_efac        !: multiplication factor for evaporation
+   REAL(wp), PUBLIC ::   rn_vfac        !: multiplication factor for ice/ocean velocity in the calculation of wind stress
    REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
@@ -224,5 +224,5 @@
       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
-      slf_i(jp_slp)  = sn_slp
+      slf_i(jp_slp ) = sn_slp
       slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj  
       !
@@ -506,5 +506,5 @@
 !!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) ) )
+               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)
                pcd_du(ji,jj) = zztmp
@@ -615,4 +615,5 @@
       zz2 = rn_pfac * rcp
       zz3 = rn_pfac * rcpi
+      zztmp = 1._wp - albo
       qns(:,:) = (  zqlw(:,:) - psen(:,:) - zqla(:,:)                          &   ! Downward Non Solar
          &        - psnow(:,:) * zztmp                                         &   ! remove latent melting heat for solid precip
@@ -788,6 +789,6 @@
    !!   'key_si3'                                       SI3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   blk_ice_tau : provide the air-ice stress
-   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_1   : provide the air-ice stress
+   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
    !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
    !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
@@ -795,7 +796,8 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE blk_ice_tau
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE blk_ice_tau  ***
+   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice,       &    ! inputs
+      &                  putaui, pvtaui, pseni, pevpi, ptsui, pssqi, pcd_dui   )    ! outputs
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_1  ***
       !!
       !! ** Purpose :   provide the surface boundary condition over sea-ice
@@ -805,8 +807,23 @@
       !!                NB: ice drag coefficient is assumed to be a constant
       !!---------------------------------------------------------------------
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndi   ! atmospheric wind at T-point [m/s]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndj   ! atmospheric wind at T-point [m/s]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric wind at T-point [m/s]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   phumi   ! atmospheric wind at T-point [m/s]      
+      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(  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_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)
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui     ! transfer coefficient for momentum      (tau)
       !!---------------------------------------------------------------------
       !
@@ -822,6 +839,6 @@
       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) )  )
+            zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  )
+            zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  )
             wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
          END DO
@@ -842,32 +859,45 @@
       ! 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))
-
-      ! ------------------------------------------------------------ !
-      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
-      ! ------------------------------------------------------------ !
-      ! 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) * 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) * 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) )
+      zrhoa  (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) )
+      zcd_dui(:,:) = wndm_ice(:,:) * Cd_atm(:,:)
+      
+      IF( ln_blk ) THEN 
+         ! ------------------------------------------------------------ !
+         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
+         ! ------------------------------------------------------------ !
+         ! C-grid ice dynamics :   U & V-points (same as ocean)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               putaui(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * zcd_dui(ji+1,jj)             &
+                  &                      + zrhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          &
+                  &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) )
+               pvtaui(ji,jj) = 0.5_wp * (  zrhoa(ji,jj+1) * zcd_dui(ji,jj+1)             &
+                  &                      + zrhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          &
+                  &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) )
+            END DO
          END DO
-      END DO
-      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
-      !
-      !
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
-         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
-      ENDIF
-      !
-   END SUBROUTINE blk_ice_tau
-
-
-   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE blk_ice_flx  ***
+         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
+         !
+         IF(ln_ctl)   CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
+            &                     , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
+      ELSE
+         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)
+               pssqi  (ji,jj) = 11637800.0_wp * EXP( -5897.8_wp / (ptsui(ji,jj)+rt0) ) / zrhoa(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      IF(ln_ctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
+      !
+   END SUBROUTINE blk_ice_1
+
+
+   SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow  )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_2  ***
       !!
       !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
@@ -883,9 +913,15 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   ptair
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   phumi
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pslp
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pqlw
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pprec
+      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   psnow
       !!
       INTEGER  ::   ji, jj, jl               ! dummy loop indices
       REAL(wp) ::   zst3                     ! local variable
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
-      REAL(wp) ::   zztmp, z1_rLsub           !   -      -
+      REAL(wp) ::   zztmp, zztmp2, z1_rLsub          !   -      -
       REAL(wp) ::   zfr1, zfr2               ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
@@ -901,5 +937,5 @@
       zcoef_dqla = -Ls * 11637800. * (-5897.8)
       !
-      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) )
       !
       zztmp = 1. / ( 1. - albo )
@@ -919,5 +955,5 @@
                qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
                ! Long  Wave (lw)
-               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
+               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
@@ -927,14 +963,15 @@
                ! ----------------------------!
 
-               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
+               ! ... turbulent heat fluxes with Ch_atm 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) - 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) - ptair(ji,jj))
                ! Latent Heat
-               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 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
+               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) *  &
+                  &                ( 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 * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
-                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
+                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_atm(ji,jj) * wndm_ice(ji,jj) *  &
+                     &                 z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2
                ELSE
                   dqla_ice(ji,jj,jl) = 0._wp
@@ -943,5 +980,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)
-
+               
                ! ----------------------------!
                !     III    Total FLUXES     !
@@ -957,8 +994,8 @@
       END DO
       !
-      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
-      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
-      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
-      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
+      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
+      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
+      CALL iom_put( 'snowpre', sprecip )                  ! Snow precipitation
+      CALL iom_put( 'precip' , tprecip )                  ! Total precipitation
 
       ! --- evaporation --- !
@@ -977,9 +1014,9 @@
       ! --- heat flux associated with emp --- !
       qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
-         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
+         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair
          &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
       qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
 
       ! --- total solar and non solar fluxes --- !
@@ -989,6 +1026,6 @@
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
-
+      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+      
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
       DO jl = 1, jpl
@@ -1016,7 +1053,7 @@
          CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
          CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
-      ENDIF
-      !
-   END SUBROUTINE blk_ice_flx
+      END IF
+      !
+   END SUBROUTINE blk_ice_2
    
 
