Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref	(revision 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref	(revision 11334)
@@ -301,5 +301,5 @@
    rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented)
    rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities)
-   ln_smth_pblh  = .false.    !  Smoothing of PBL height with a 2d Hanning filter
+   ln_smth_pblh  = .true.    !  Smoothing of PBL height with a 2d Hanning filter
 /
 !-----------------------------------------------------------------------
@@ -589,13 +589,12 @@
 !-----------------------------------------------------------------------
    ln_tide     = .false.      ! Activate tides
-      ln_tide_pot   = .true.                !  use tidal potential forcing
-         ln_scal_load  = .false.               ! Use scalar approximation for
-            rn_scal_load = 0.094               !     load potential
-         ln_read_load  = .false.               ! Or read load potential from file
-            cn_tide_load = 'tide_LOAD_grid_T.nc'  ! filename for load potential
-            !      
-      ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup
-         rdttideramp   =    0.                 !  ramp duration in days
-      clname(1)     = 'DUMMY'               !  name of constituent - all tidal components must be set in namelist_cfg
+   ln_tide_pot   = .true.                !  use tidal potential forcing
+   ln_scal_load  = .false.               ! Use scalar approximation for 
+   ln_read_load  = .false.               ! Or read load potential from file
+   cn_tide_load  = 'tide_LOAD_grid_T.nc'  ! filename for load potential     
+   ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup
+   rn_scal_load  = 0.094               !     load potential 
+   rdttideramp   =    0.                 !  ramp duration in days
+   clname        = 'FOUR'               !  name of constituent - all tidal components must be set in namelist_cfg
 /
 !-----------------------------------------------------------------------
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
 
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 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90	(revision 11334)
@@ -75,8 +75,5 @@
             &                                      u_ice, v_ice,   &    ! inputs
             &                                      putaui = utau_ice, pvtaui = vtau_ice )    ! outputs                             
-         CASE( jp_abl     )   ;    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  
+ !        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
       END SELECT
@@ -146,5 +143,5 @@
       CASE( jp_usr )              !--- user defined formulation
                                   CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
-      CASE( jp_blk, jp_abl )  !--- bulk formulation (temporary version for ABL option)
+      CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation
                                   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) )    ! 
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90	(revision 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90	(revision 11334)
@@ -114,13 +114,18 @@
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
-
+      !
+      ALLOCATE( z1_at_i(jpi,jpj) )
+      WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
+      ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
+      END WHERE
+      !
+      tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
+      WHERE( at_i(:,:)<=epsi20 ); tm_su(:,:) = rt0; END WHERE      
+      !
       ! The following fields are calculated for diagnostics and outputs only
       ! ==> Do not use them for other purposes
       IF( kn > 1 ) THEN
          !
-         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
-         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
-         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
-         END WHERE
+         ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
          WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
          ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
@@ -135,5 +140,4 @@
          !         
          !                          ! mean temperature (K), salinity and age
-         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
@@ -153,5 +157,4 @@
          !                           ! put rt0 where there is no ice
          WHERE( at_i(:,:)<=epsi20 )
-            tm_su(:,:) = rt0
             tm_si(:,:) = rt0
             tm_i (:,:) = rt0
@@ -159,6 +162,8 @@
          END WHERE
 
-         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
+         DEALLOCATE( z1_vt_i , z1_vt_s )
       ENDIF
+      !
+      DEALLOCATE( z1_at_i )
       !
    END SUBROUTINE ice_var_agg
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 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11334)
@@ -26,4 +26,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers 
+   USE abl            ! abl variables in case ln_abl = .true.
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -66,13 +67,17 @@
    PUBLIC   dia_wri_state
    PUBLIC   dia_wri_alloc           ! Called by nemogcm module
-
+#if ! defined key_iomput   
+   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.)
+#endif
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
-   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
+   INTEGER ::          nb_T, ndim_bT                 ! grid_T file
    INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
    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 ::   nid_A, nz_A, nh_A                    ! grid_ABL file   
    INTEGER ::   ndex(1)                              ! ???
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
+   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_bT
 
@@ -409,11 +414,16 @@
          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
-         !
-      dia_wri_alloc = MAXVAL(ierr)
+	  dia_wri_alloc = MAXVAL(ierr)
       CALL mpp_sum( 'diawri', dia_wri_alloc )
       !
    END FUNCTION dia_wri_alloc
-
-   
+ 
+   INTEGER FUNCTION dia_wri_alloc_abl()
+      !!----------------------------------------------------------------------
+	  ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl)
+      CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
+      !
+   END FUNCTION dia_wri_alloc_abl
+ 
    SUBROUTINE dia_wri( kt )
       !!---------------------------------------------------------------------
@@ -435,5 +445,5 @@
       INTEGER  ::   ji, jj, jk                               ! dummy loop indices
       INTEGER  ::   ierr                                     ! error code return from allocation
-      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma, ipka   ! local integers
       INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
@@ -441,4 +451,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
       !!----------------------------------------------------------------------
       ! 
@@ -449,4 +460,5 @@
          ninist = 0
       ENDIF
+	  
       !
       ! 0. Initialisation
@@ -472,4 +484,5 @@
       ijmi = 1      ;      ijma = jpj
       ipk = jpk
+	  IF(ln_abl) ipka = jpkam1
 
       ! define time axis
@@ -574,4 +587,20 @@
             &          "m", ipk, gdepw_1d, nz_W, "down" )
 
+         IF( ln_abl ) THEN 
+         ! Define the ABL grid FILE ( nid_A )
+            CALL dia_nam( clhstnam, nwrite, 'grid_ABL' )
+            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
+               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+               &          nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
+            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
+               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
+            !                                                            ! Index of ocean points
+			ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
+			zw3d_abl(:,:,:) = 1._wp 
+			CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
+            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
+			DEALLOCATE(zw3d_abl)
+         ENDIF
 
          ! Declare all the output fields as NETCDF variables
@@ -623,4 +652,30 @@
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+!
+         IF( ln_abl ) THEN
+         !                                                                                      !!! nid_A : 3D
+			CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )		 			   
+#if defined key_si3
+            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
+#endif
+		    CALL histend( nid_A, snc4chunks=snc4set )
+		 !
+		 ENDIF
 !
          IF( ln_icebergs ) THEN
@@ -676,6 +731,6 @@
          
          clmx ="l_max(only(x))"    ! max index on a period
-!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
-!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
+         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
 #if defined key_diahth
          CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
@@ -790,4 +845,27 @@
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
+!
+      IF( ln_abl ) THEN 
+	     ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
+	     IF( ln_mskland )   THEN 
+		    DO jk=1,jpka
+		       zw3d_abl(:,:,jk) = tmask(:,:,1)
+            END DO
+		 ELSE
+            zw3d_abl(:,:,:) = 1._wp		 
+         ENDIF			
+		 CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh 
+	     CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
+	     CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
+	     CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
+	     CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl		 
+	     CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
+	     CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
+	     CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl	
+#if defined key_si3
+         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
+#endif
+		 DEALLOCATE(zw3d_abl)
+	  ENDIF
 !
       IF( ln_icebergs ) THEN
@@ -826,6 +904,6 @@
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
-!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
-!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
+      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
+      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
 
 #if defined key_diahth
@@ -862,4 +940,5 @@
          CALL histclo( nid_V )
          CALL histclo( nid_W )
+         IF(ln_abl) CALL histclo( nid_A )
       ENDIF
       !
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/abl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/abl.F90	(revision 11334)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/abl.F90	(revision 11334)
@@ -0,0 +1,31 @@
+MODULE abl
+   !!======================================================================
+   !!                      ***  MODULE  abl  ***
+   !! Abl        :  ABL dynamics and active tracers defined in memory 
+   !!======================================================================
+   USE par_kind        ! abl parameters
+  
+   IMPLICIT NONE
+   PRIVATE
+   !! --------------------------                            ! 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   u_abl        !: i-horizontal velocity   [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   v_abl        !: j-horizontal velocity   [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:,:) ::   tq_abl       !: 4D T-q fields           [Kelvin,kg/kg] 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avm_abl      !: turbulent viscosity   [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avt_abl      !: turbulent diffusivity [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxl_abl      !: mixing length         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   tke_abl      !: turbulent kinetic energy [m2/s2]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   fft_abl      !: Coriolis parameter    [1/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   pblh         !: PBL height            [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   rest_eq
+   !
+   INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2)
+   ! 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+END MODULE abl
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90	(revision 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90	(revision 11334)
@@ -8,4 +8,5 @@
    !!----------------------------------------------------------------------
    USE sbc_oce, ONLY :   ght_abl, ghw_abl, e3t_abl, e3w_abl
+   USE lib_mpp, ONLY :   ctl_stop
 
    IMPLICIT NONE
@@ -31,5 +32,7 @@
       !!
       !!----------------------------------------------------------------------
-      
+      CALL ctl_stop( 'STOP', 'ln_abl = .true. but ABL source directory was not included',   &
+	     & '(Either switch to ln_abl = .false. or modify your cfg.txt file and recompile)' )
+	  !!
    END SUBROUTINE sbc_abl_init
 
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 11322)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11334)
@@ -816,5 +816,5 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
+      REAL(wp) ::   zwndi_t , zwndj_t, zootm_su, zztmp1, zztmp2 ! 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)
@@ -855,5 +855,5 @@
       zcd_dui(:,:) = wndm_ice(:,:) * Cd_atm(:,:)
       
- !     IF( ln_blk ) THEN 
+      IF( ln_blk ) THEN 
          ! ------------------------------------------------------------ !
          !    Wind stress relative to the moving ice ( U10m - U_ice )   !
@@ -874,14 +874,17 @@
          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
+      ELSE
+         zztmp1 = 11637800.0_wp
+		 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)
+            END DO
+         END DO
+      ENDIF
       !
       IF(ln_ctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
@@ -1248,6 +1251,7 @@
       REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
       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(:,:)
