Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90	(revision 12015)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90	(revision 12021)
@@ -140,5 +140,5 @@
       !!             ***  ROUTINE sbc_blk_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( t_zu(jpi,jpj), q_zu(jpi,jpj),                                            &
+      ALLOCATE( t_zu(jpi,jpj)   , q_zu(jpi,jpj)   , tsk(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 )
@@ -206,7 +206,12 @@
 
       !                             !** initialization of the cool-skin / warm-layer parametrization
-      IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) &
-         & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
-      !
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
+         !! Some namelist sanity tests:
+         IF( ln_NCAR )      &
+            & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )         
+         IF( nn_fsbc /= 1 ) &
+            & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')                  
+      END IF
+      
       ioptio = 0
       IF( ln_humi_sph ) THEN
@@ -351,13 +356,14 @@
       !!              (momentum, heat, freshwater and runoff)
       !!
-      !! ** Method  : (1) READ each fluxes in NetCDF files:
-      !!      the 10m wind velocity (i-component) (m/s)    at T-point
-      !!      the 10m wind velocity (j-component) (m/s)    at T-point
-      !!      the 10m or 2m specific humidity     (kg/kg)
-      !!      the solar heat                      (W/m2)
-      !!      the Long wave                       (W/m2)
-      !!      the 10m or 2m air temperature       (Kelvin)
-      !!      the total precipitation (rain+snow) (Kg/m2/s)
-      !!      the snow (solid precipitation)      (kg/m2/s)
+      !! ** Method  : 
+      !!              (1) READ each fluxes in NetCDF files:
+      !!      the wind velocity (i-component) at z=rn_zu  (m/s) at T-point
+      !!      the wind velocity (j-component) at z=rn_zu  (m/s) at T-point
+      !!      the specific humidity           at z=rn_zqt (kg/kg)
+      !!      the solar heat                  at z=rn_zqt (W/m2)
+      !!      the Long wave                               (W/m2)
+      !!      the air temperature                         (Kelvin)
+      !!      the total precipitation (rain+snow)         (Kg/m2/s)
+      !!      the snow (solid precipitation)              (kg/m2/s)
       !!      ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 
       !!              (2) CALL blk_oce_1 and blk_oce_2
@@ -556,7 +562,11 @@
 
 
-      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
-
-         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
+      !! Time to call the user-selected bulk parameterization for
+      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more...
+      SELECT CASE( nblk )        
+
+         CASE( np_NCAR      )
+            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                              &
+               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_COARE_3p0 )
@@ -576,7 +586,9 @@
 
          CASE DEFAULT
-            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' )
-         END SELECT
-
+         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
+
+      END SELECT
+
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq:
          WHERE ( fr_i < 0.001_wp )
@@ -589,33 +601,6 @@
             pssq(:,:) = zztmp2(:,:) !LB:  "   "   "
          END WHERE
-
-         !LB: Update of tsk, the "official" array for skin temperature
-         tsk(:,:) = zst(:,:)
-
-      ELSE !IF( ln_skin_cs .OR. ln_skin_wl )
-
-         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
-            !
-         CASE( np_NCAR      )
-            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                                &
-               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-
-         CASE( np_COARE_3p0 )
-            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,   &
-               &                 zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-
-         CASE( np_COARE_3p6 )
-            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,   &
-               &                 zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-
-         CASE( np_ECMWF     )
-            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &
-               &                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' )
-         END SELECT
-
-      ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl )
+         tsk(:,:) = zst(:,:) !#LB: Update of tsk, the "official" array for skin temperature
+      END IF
 
       !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef.
@@ -632,14 +617,16 @@
       ! -------------------------------------------------------------
       
-      IF (ln_abl) THEN
-
-         CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
-            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         &
-            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 &
-            &               pcd_du(:,:), psen(:,:), pevp(:,:),                &
-            &               prhoa=rhoa(:,:) )
-
-      ELSE
-
+      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
+!! FL do we need this multiplication by tmask ... ???
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               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 * zcd_oce(ji,jj)
+               psen(ji,jj)   = zztmp * zch_oce(ji,jj)
+               pevp(ji,jj)   = zztmp * zce_oce(ji,jj)       
+            END DO
+         END DO
+      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
          CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
             &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         &
@@ -673,4 +660,5 @@
          END DO
          CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+
          IF(ln_ctl) THEN
             CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ')
@@ -738,5 +726,5 @@
       ! -----------------------------
 
-      zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:)     ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
+      zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) * -1._wp    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
 
       IF(ln_ctl) THEN
@@ -866,5 +854,5 @@
 
       ! local scalars ( place there for vector optimisation purposes)
-      IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)
+      !IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)
       zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
       
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90	(revision 12015)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90	(revision 12021)
@@ -20,5 +20,4 @@
 
    USE dom_oce        ! ocean space and time domain
-   USE sbc_oce, ONLY  : ln_abl !!GS
    USE phycst         ! physical constants
 
@@ -575,21 +574,4 @@
       INTEGER  :: ji, jj, jq     ! dummy loop indices
       !!----------------------------------------------------------------------------------
-
-      !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
-      IF (ln_abl) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-   
-               zUrho = pUb(ji,jj)
-   
-               pTau(ji,jj)  = zUrho * pCd(ji,jj)   !!GS: store cd_du in wind stress
-   
-               zevap        = zUrho * pCe(ji,jj)
-               pQsen(ji,jj) = zUrho * pCh(ji,jj)
-               pQlat(ji,jj) = zevap                !!GS: store evap in latent
-   
-            END DO
-         END DO
-      ELSE
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -617,6 +599,4 @@
             END DO
          END DO
-      END IF
-
    END SUBROUTINE BULK_FORMULA_VCTR
 
@@ -650,16 +630,4 @@
       !!----------------------------------------------------------------------------------
 
-      !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
-      IF (ln_abl) THEN
-   
-         zUrho = pUb
-   
-         pTau  = zUrho * pCd
-   
-         zevap = zUrho * pCe
-         pQsen = zUrho * pCh
-         pQlat = zevap
-
-      ELSE
          !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
          ztaa = pTa ! first guess...
@@ -681,5 +649,4 @@
          IF( PRESENT(pEvap) ) pEvap = - zevap
          IF( PRESENT(prhoa) ) prhoa = zrho
-     END IF
 
    END SUBROUTINE BULK_FORMULA_SCLR
