Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11803)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11804)
@@ -520,5 +520,5 @@
       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
-      !!    (since reanalysis products provide T at z, not theta !)
+      !!    (since reanalysis products provide absolute temperature "T" at z, not theta !)
       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
 
@@ -594,26 +594,34 @@
       END IF ! IF ( ln_skin_cs .OR. ln_skin_wl )
 
-
-
-
-      !                          ! Compute true air density :
-      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove rhoa*grav*rn_zu from SLP...)
-         rhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
-      ELSE                                      ! At zt:
-         rhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), zqair(:,:), sf(jp_slp)%fnow(:,:,1) )
-      END IF
-
       !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
       !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
-
-      DO jj = 1, jpj             ! tau module, i and j component
+      
+      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
+         !! If zu == zt, then ensuring once for all that:
+         t_zu(:,:) = ztpot(:,:)
+         q_zu(:,:) = zqair(:,:)
+      END IF
+      
+
+      !  Turbulent fluxes over ocean  => TURB_FLUXES @ sbcblk_phy.F90
+      ! -------------------------------------------------------------
+      
+      CALL TURB_FLUXES( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:), &
+         &                 wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1), &
+         &                 taum(:,:), zqsb(:,:), zqla(:,:),               &
+         &                 pEvap=zevap(:,:), prhoa=rhoa(:,:) )
+      
+      zqla(:,:)  =  zqla(:,:) * tmask(:,:,1)
+      zqsb(:,:)  =  zqsb(:,:) * tmask(:,:,1)
+      taum(:,:)  =  taum(:,:) * tmask(:,:,1)
+      zevap(:,:) = zevap(:,:) * tmask(:,:,1)
+      
+      DO jj = 1, jpj             ! tau i and j component on T-grid points
          DO ji = 1, jpi
-            zztmp = rhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
-            taum  (ji,jj) = zztmp * wndm  (ji,jj)
+            zztmp = taum(ji,jj) / wndm(ji,jj)
             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
             zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
          END DO
       END DO
-
       !                          ! add the HF tau contribution to the wind stress module
       IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
@@ -634,23 +642,4 @@
       CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
 
-      !  Turbulent fluxes over ocean
-      ! -----------------------------
-
-      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
-      zqla(:,:) = rhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
-
-      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
-         !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,  zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - zqair(:,:)) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)  ) ! Sensible Heat, using bulk wind speed
-      ELSE
-         !! q_air and t_air are not given at 10m (wind reference height)
-         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
-         zevap(:,:) = rn_efac*MAX( 0._wp,  zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:)) )  ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
-      ENDIF
-
-      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
-
 
       ! ----------------------------------------------------------------------------- !
@@ -658,5 +647,5 @@
       ! ----------------------------------------------------------------------------- !
 
-      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
+      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
       zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
 
@@ -681,5 +670,5 @@
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
       !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
+      qns(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
@@ -691,5 +680,5 @@
       !
 #if defined key_si3
-      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
+      qns_oce(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                ! non solar without emp (only needed by SI3)
       qsr_oce(:,:) = qsr(:,:)
 #endif
@@ -698,8 +687,8 @@
       CALL iom_put( "rho_air"  ,   rhoa )                 ! output air density (kg/m^3) !#LB
       CALL iom_put( "qlw_oce"  ,   zqlw )                 ! output downward longwave heat over the ocean
-      CALL iom_put( "qsb_oce"  , - zqsb )                 ! output downward sensible heat over the ocean
-      CALL iom_put( "qla_oce"  , - zqla )                 ! output downward latent   heat over the ocean
+      CALL iom_put( "qsb_oce"  ,   zqsb )                 ! output downward sensible heat over the ocean
+      CALL iom_put( "qla_oce"  ,   zqla )                 ! output downward latent   heat over the ocean
       CALL iom_put( "evap_oce" ,  zevap )                 ! evaporation
-      CALL iom_put( "qemp_oce" ,   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
+      CALL iom_put( "qemp_oce" ,   qns-zqlw-zqsb-zqla )   ! output downward heat content of E-P over the ocean
       CALL iom_put( "qns_oce"  ,   qns  )                 ! output downward non solar heat over the ocean
       CALL iom_put( "qsr_oce"  ,   qsr  )                 ! output downward solar heat over the ocean
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11803)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11804)
@@ -495,5 +495,5 @@
       REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
       !
-      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zUrho, zTs2, zz0, &
+      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, &
          &        zQlat, zQsen, zQlw
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -508,19 +508,10 @@
             zCh = zz0*ptst(ji,jj)/zdt
             zCe = zz0*pqst(ji,jj)/zdq
-
-            !zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10
-            zTs2  = pTs(ji,jj)*pTs(ji,jj)
-
+            
             CALL TURB_FLUXES( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
                &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), &
                &              pTau(ji,jj), zQsen, zQlat )
-
-
-            ! Wind stress module:
-            !pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo?
-
-            ! Non-Solar heat flux to the ocean:
-            !zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp )  ! we do not want a Qlat > 0 !
-            !zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt
+            
+            zTs2  = pTs(ji,jj)*pTs(ji,jj)            
             zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
 
@@ -533,9 +524,6 @@
 
 
-
-
-
    SUBROUTINE TURB_FLUXES_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
-      &                                 pTau, pQsen, pQlat,  pEvap )
+      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
       !!----------------------------------------------------------------------------------
       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
@@ -555,5 +543,6 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
       !!
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap !  [kg/m^2/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
       !!
       REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
@@ -576,9 +565,10 @@
             pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
 
-            zevap        = MIN( zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) , 0._wp )   ! we do not want condensation & Qlat > 0 !
-            pQsen(ji,jj) =      zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
-            pQlat(ji,jj) =  L_vap(pTs(ji,jj)) * zevap
+            zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
+            pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
+            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
 
             IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
+            IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
 
          END DO
@@ -588,5 +578,5 @@
 
    SUBROUTINE TURB_FLUXES_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
-      &                                 pTau, pQsen, pQlat,  pEvap )
+      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
       !!----------------------------------------------------------------------------------
       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
@@ -606,5 +596,6 @@
       REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
       !!
-      REAL(wp), INTENT(out), OPTIONAL :: pEvap !  [kg/m^2/s]
+      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
+      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
       !!
       REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
@@ -625,12 +616,12 @@
       pTau = zUrho * pCd * pwnd ! Wind stress module
 
-      zevap        = MIN( zUrho * pCe * (pqa - pqs) , 0._wp )   ! we do not want condensation & Qlat > 0 !
-      pQsen =      zUrho * pCh * (pTa - pTs) * cp_air(pqa)
-      pQlat =  L_vap(pTs) * zevap
+      zevap = zUrho * pCe * (pqa - pqs)
+      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
+      pQlat = L_vap(pTs) * zevap
 
       IF ( PRESENT(pEvap) ) pEvap = - zevap
+      IF ( PRESENT(prhoa) ) prhoa = zrho
 
    END SUBROUTINE TURB_FLUXES_SCLR
-
 
 
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 11803)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 11804)
@@ -63,5 +63,5 @@
       !!---------------------------------------------------------------------
       !!
-      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
+      !! Cool-skin parameterization, based on Fairall et al., 1996, revisited for COARE 3.6 (Fairall et al., 2019)
       !!
       !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
@@ -78,11 +78,10 @@
       !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
       !!     *pQlat*      surface latent heat flux                       [K]
-      !!
       !!------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST ! bulk SST [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQlat  ! latent heat flux [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jc
@@ -92,16 +91,16 @@
          DO ji = 1, jpi
 
-            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$
-            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$
+            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
+            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
 
             zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
 
             DO jc = 1, 4 ! because implicit in terms of zdelta...
-               zfr    = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
-               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface
+               zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
+               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
                zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
             END DO
 
-            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment
+            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
 
          END DO
@@ -189,6 +188,4 @@
                   &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
                zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!!
-               !PRINT *, '#LBD:  Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj))
-               !PRINT *, '#LBD:       =>Qabs:', zQabs,' zfr=', zfr
 
                IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
@@ -196,7 +193,4 @@
                   ! since zQabs <= 0._wp
                   ! => no need to go further
-                  !PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)'
-                  !PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs
-                  !PRINT *, '#LBD: => leaving without changing anything...'
                   l_exit = .TRUE.
                END IF
@@ -207,6 +201,4 @@
             !LOLO: remove??? has a strong influence !!!
             IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
-               !PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt
-               !PRINT *, '#LBD:  => time to destroy the warm-layer!'
                l_exit       = .TRUE.
                l_destroy_wl = .TRUE.
@@ -219,9 +211,4 @@
                ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
                ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !
-
-               !PRINT *, '#LBD:======================================================'
-               !PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4)
-               !PRINT *, '#LBD:======================================================'
-               !PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
 
                ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
@@ -238,7 +225,4 @@
                   zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
                END DO
-               !PRINT *, '#LBD: updated absorption and WL depth=',  REAL(zfr,4), REAL(zHwl,4)
-               !PRINT *, '#LBD: updated value for Qabs=',  REAL(zQabs,4), 'W/m2'
-               !PRINT *, '#LBD: updated value for Qac =',  REAL(zqac,4), 'J'
 
                IF ( zqac <= 0._wp ) THEN
@@ -263,15 +247,10 @@
             END IF
 
-            !PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4)
-
             IF ( iwait == 0 ) THEN
                !! Iteration loop within bulk algo is over, time to update what needs to be updated:
                dT_wl(ji,jj)  = zdTwl
                Hz_wl(ji,jj)  = zHwl
-               !PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4)
                Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
                Tau_ac(ji,jj) = ztac
-               !PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
-               !PRINT *, '#LBD'
             END IF
 
@@ -284,5 +263,5 @@
 
 
-   FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a )
+   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a )
       !!---------------------------------------------------------------------
       !! Computes the thickness (m) of the viscous skin layer.
@@ -298,20 +277,25 @@
       REAL(wp)                :: delta_skin_layer
       REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
-      REAL(wp), INTENT(in)    :: pQabs    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
+      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
       REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
       REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
       !!---------------------------------------------------------------------
-      REAL(wp) :: zusw, zusw2, zlamb, zQb
-      !!---------------------------------------------------------------------
-
-      zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
-
+      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
+      !!---------------------------------------------------------------------
+      
+      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
+
+      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
+      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
+      !
       zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
       zusw2 = zusw*zusw
-
-      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996
-
-      delta_skin_layer = zlamb*rnu0_w/zusw
-
+      !
+      zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
+      !  => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
+      !
+      ztmp = rnu0_w/zusw
+      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
+         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
    END FUNCTION delta_skin_layer
 
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11803)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11804)
@@ -88,19 +88,17 @@
          DO ji = 1, jpi
 
-            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$
-            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$
-            !zQabs  = pQnsol(ji,jj)
+            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
+            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
 
             zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
 
             DO jc = 1, 4 ! because implicit in terms of zdelta...
-               zfr    = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
+               zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
                !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
-               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface
-               !zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
+               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
                zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
             END DO
 
-            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment
+            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
 
          END DO
@@ -142,5 +140,5 @@
          & zalpha, & !: thermal expansion coefficient of sea-water [1/K]
          & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
-         & zfr, zeta, ztmp, &
+         & zfr, zeta, &
          & zusw, zusw2, &
          & zLa, zfLa, &
@@ -239,5 +237,5 @@
 
 
-   FUNCTION delta_skin_layer( palpha, pQabs, pustar_a )
+   FUNCTION delta_skin_layer( palpha, pQd, pustar_a )
       !!---------------------------------------------------------------------
       !! Computes the thickness (m) of the viscous skin layer.
@@ -253,28 +251,22 @@
       REAL(wp)                :: delta_skin_layer
       REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
-      REAL(wp), INTENT(in)    :: pQabs    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
+      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
       REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
       !!---------------------------------------------------------------------
-      REAL(wp) :: zusw, zusw2, zlamb, zQb
-      !!---------------------------------------------------------------------
-
-      zQb = pQabs
-
-      !zQ = MIN( -0.1_wp , pQabs )
-
-      !ztf = 0.5_wp + SIGN(0.5_wp, zQ)  ! Qabs < 0 => cooling of the layer => ztf = 0 (normal case)
-      !                                   ! Qabs > 0 => warming of the layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
+      REAL(wp) :: zusw, zusw2, zlamb, ztf, ztmp
+      !!---------------------------------------------------------------------
+      ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
+      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
+      !
       zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
       zusw2 = zusw*zusw
-
-      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996
-      !zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQ, 0._wp)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996
-
-      delta_skin_layer = zlamb*rnu0_w/zusw
-
-      !delta_skin_layer =  (1._wp - ztf) * zlamb*rnu0_w/zusw    &         ! see eq.(12) in Fairall et al., 1996
-      !   &               +     ztf  * MIN(6._wp*rnu0_w/zusw , 0.007_wp)
+      !
+      zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*pQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
+      !  => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
+      !
+      ztmp = rnu0_w/zusw
+      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
+         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
    END FUNCTION delta_skin_layer
-
 
 
