Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11334)
@@ -38,8 +38,14 @@
 !===================================================================================================
    SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in
-              &            pu_dta, pv_dta, pt_dta, pq_dta, &
-              &            pslp_dta, pgu_dta, pgv_dta, &
-              &            pcd_du, psen, pevp,     &                            ! in/out
-              &            pwndm, ptaui, ptauj, ptaum )
+              &            pu_dta, pv_dta, pt_dta, pq_dta,    & 
+              &            pslp_dta, pgu_dta, pgv_dta,        &
+              &            pcd_du, psen, pevp,                &                            ! in/out
+              &            pwndm, ptaui, ptauj, ptaum         & 
+#if defined key_si3			  
+              &          , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice  & 
+              &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      & 
+			  &          , ptaui_ice,     ptauj_ice                      &
+#endif			  
+			  &   )
 !---------------------------------------------------------------------------------------------------
 
@@ -80,9 +86,21 @@
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau|| 
       !
+#if defined key_si3	  
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point) 	  
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity 
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point)
+	  REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point)
+	  REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
+	  REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||  
+	  REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pfrac_oce	   !
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point)
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point) 	  
+#endif	  
+	  !
       REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zrhoa, zwnd_i, zwnd_j
-!      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   zFC
       REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF    
-!      REAL(wp), DIMENSION(1:jpi,  jptq  )        ::   zBC
-      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)     ::   z_cft      !--FL--to be removed after the test phase   
+      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase   
       !
       REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
@@ -91,6 +109,6 @@
       !
       INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
-      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1
-      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv
+      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
+      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
       !
       !!---------------------------------------------------------------------      
@@ -108,5 +126,11 @@
       DO jj = 1,jpj
          DO ji = 1,jpi
-            ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj)
+		    zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
+#if defined key_si3				
+			zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
+		    ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
+#else
+	        ustar2(ji,jj) = zzoce	
+#endif						
          END DO
       END DO  
@@ -154,14 +178,24 @@
             IF(jtra == jp_ta) THEN
                DO ji = 1,jpi  ! surface boundary condition for temperature              
-                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj) 
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
-                    &                               + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )                
-                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
+                  zztmp1 = psen(ji, jj) 
+				  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )	
+#if defined key_si3				  
+                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
+                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
+#endif				  
+				  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1				  
+                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2               
+				  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
                END DO 
             ELSE   
                DO ji = 1,jpi  ! surface boundary condition for humidity              
-                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj) 
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
-                    &                               + rdt * pevp(ji, jj) * pssq(ji, jj)               
+                  zztmp1 = pevp(ji, jj) 
+				  zztmp2 = pevp(ji, jj) * pssq(ji, jj)   
+#if defined key_si3				  
+                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 
+				  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)    
+#endif	
+                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1
+                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2               
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
                END DO                                     
@@ -251,5 +285,5 @@
          ENDDO 
       END IF
-
+	  
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  4 *** Advance u,v to time n+1 
@@ -273,7 +307,15 @@
             z_elem_a( ji,     2    ) = 0._wp
             z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
-            zcff  = rdt * pcd_du(ji, jj)
-            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
-            u_abl( ji, jj,    2, nt_a ) = u_abl( ji, jj,    2, nt_a ) + 0.5_wp * zcff * ( pssu(ji-1, jj) + pssu(ji,jj) ) 
+            !
+			zztmp1  = pcd_du(ji, jj)
+			zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )			
+#if defined key_si3				  
+            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
+			zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) )
+			zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
+#endif           
+			z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1 			
+            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rdt * zztmp2
+			
             !++ Top Neumann B.C.
             !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 
@@ -331,8 +373,15 @@
             !++ Surface boundary condition
             z_elem_a( ji,     2    ) = 0._wp
-            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
-            zcff  = rdt * pcd_du(ji, jj)
-            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
-            v_abl( ji, jj,    2, nt_a ) =  v_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssv(ji, jj) + pssv(ji, jj-1) )  
+            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )        
+            !
+			zztmp1 = pcd_du(ji, jj)
+			zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )  
+#if defined key_si3				  
+            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
+			zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) )
+			zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
+#endif 			
+            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1 
+            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2
             !++ Top Neumann B.C.            
             !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
@@ -372,5 +421,5 @@
       END DO             ! end outer loop 
       !-------------   
-  
+    
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  5 *** Apply nudging on the dynamics and the tracers 
@@ -428,12 +477,12 @@
       !-------------
       END DO             ! end outer loop
-      !------------- 
+      !------------- 	  
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  6 *** MPI exchanges  
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !
-      CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), kfillmode = jpfillnothing  )
+	  !
+      CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), 'T', -1., kfillmode = jpfillnothing  )
       CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing  )   ! ++++ this should not be needed...
-      !
+	  !
       CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      ))
       CALL iom_put ( 'v_abl',  v_abl(:,:,2:jpka,nt_a      ))
@@ -504,5 +553,28 @@
          CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' )
       ENDIF
-	  
+
+#if defined key_si3
+         ! ------------------------------------------------------------ !
+         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
+         ! ------------------------------------------------------------ !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1  
+               
+			   zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
+			   zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
+			   
+			   ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
+                  &                      +    zrhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
+                  &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
+               ptauj_ice(ji,jj) = 0.5_wp * (  zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             &
+                  &                      +    zrhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          &
+                  &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
+         !
+         IF(ln_ctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
+            &                     , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
+#endif
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  8 *** Swap time indices for the next timestep
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 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11334)
@@ -30,5 +30,11 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
-
+#if defined key_si3
+   USE ice    , ONLY : u_ice, v_ice, tm_su, ato_i      ! ato_i = total open water fractional area
+   USE sbc_ice, ONLY : wndm_ice, utau_ice, vtau_ice
+#endif   
+#if ! defined key_iomput
+   USE diawri    , ONLY : dia_wri_alloc_abl
+#endif
    IMPLICIT NONE
    PRIVATE
@@ -141,4 +147,8 @@
       CALL iom_get( inum, jpdom_unknown, 'ghw_abl', ghw_abl(:) )
       CALL iom_close( inum )
+	  
+#if ! defined key_iomput
+	  IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' )
+#endif
 
       IF(lwp) THEN
@@ -303,4 +313,7 @@
       !!
       REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
+#if defined key_si3
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi
+#endif  	  
       INTEGER                      ::   jbak, jbak_dta, ji, jj
       !!---------------------------------------------------------------------
@@ -319,5 +332,11 @@
               &           sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, &
               &              zssq, zcd_du, zsen, zevp )
-  
+
+#if defined key_si3
+	  CALL blk_ice_1( u_abl(:,:,2,nt_n), v_abl(:,:,2,nt_n), tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &
+	     &            sf(jp_slp)%fnow(:,:,1), u_ice, v_ice,                                          &
+         &            pseni=zseni, pevpi=zevpi, ptsui=tm_su, pssqi=zssqi, pcd_dui=zcd_dui )   ! outputs  
+#endif  
+
       !!-------------------------------------------------------------------------------------------
       !! 3 - Advance ABL variables from now (n) to after (n+1)
@@ -329,6 +348,11 @@
               &         sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),                          &
               &         zcd_du, zsen, zevp,    &                                ! in/out
-              &         wndm, utau, vtau, taum )                                ! out
-
+              &         wndm, utau, vtau, taum &                                 ! out
+#if defined key_si3			  
+              &          , tm_su, u_ice, v_ice, zssqi, zcd_dui   & 
+              &          , zseni, zevpi, wndm_ice, ato_i,        & 
+			  &            utau_ice,     vtau_ice                &
+#endif			  
+              &                                                                )
       !!-------------------------------------------------------------------------------------------
       !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 
@@ -340,4 +364,9 @@
               &           sst_m , zsen, zevp )
 
+#if defined key_si3
+! Avoid a USE abl in icesbc module
+      sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
+#endif 
+
    END SUBROUTINE sbc_abl
 
