Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11625)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11626)
@@ -38,5 +38,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_coare  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
@@ -348,5 +348,5 @@
                &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
 
-            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
+            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
 
             T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)
@@ -360,9 +360,9 @@
             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
                &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
-            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
+            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
             IF (PRESENT(Hwl)) THEN
-               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
+               CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
             ELSE
-               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
+               CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
             END IF
             !! Updating T_s and q_s !!!
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11625)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11626)
@@ -38,5 +38,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_coare  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
@@ -348,5 +348,5 @@
                &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
 
-            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
+            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
 
             T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)
@@ -360,9 +360,9 @@
             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
                &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
-            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
+            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
             IF (PRESENT(Hwl)) THEN
-               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
+               CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
             ELSE
-               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
+               CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
             END IF
             !! Updating T_s and q_s !!!
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 11626)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 11626)
@@ -0,0 +1,301 @@
+MODULE sbcblk_skin_coare
+   !!======================================================================
+   !!                   ***  MODULE  sbcblk_skin_coare  ***
+   !! Computes:
+   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
+   !!     scheme used at ECMWF
+   !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019)
+   !!
+   !!   From routine turb_ecmwf maintained and developed in AeroBulk
+   !!            (https://github.com/brodeau/aerobulk)
+   !!
+   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!----------------------------------------------------------------------
+   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
+   !!               based on the cool-skin/warm-layer scheme used at ECMWF
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+
+   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
+
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! to use key_nosignedzero
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC :: CS_COARE, WL_COARE
+
+   !! Cool-skin related parameters:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl  !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m]
+
+   !! Warm-layer related parameters:
+   REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp    !: maximum depth of warm layer (adjustable)
+   !
+   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number
+   !
+   REAL(wp), PARAMETER :: Qabs_thr = 50._wp  !: threshold for heat flux absorbed in WL
+   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac  ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac    ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl     ! depth of warm-layer [m]
+   !!----------------------------------------------------------------------
+CONTAINS
+
+
+   SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat,  pdT )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pustar*     friction velocity u*                           [m/s]
+      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
+      !!     *pQlat*      surface latent heat flux                       [K]
+      !!
+      !!  **  INPUT/OUTPUT:
+      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
+      !!              as output =>  new estimate of dT cool-skin
+      !!
+      !!------------------------------------------------------------------
+      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(inout) :: pdT    ! dT due to cool-skin effect
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
+         &        zz1, zz2, zus, &
+         &        ztf
+      !!---------------------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
+
+            zdelta = delta_vl(ji,jj)   ! using last value of delta
+
+            !! Fraction of the shortwave flux absorbed by the cool-skin sublayer:
+            zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! 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 !
+
+            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) )  ! Total cooling at the interface
+
+            ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1
+            !                                 Qt < 0 => warming of the layer => ztf = 0
+
+            !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap):
+            zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap  ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0
+            !! LB: this terms only makes sense if > 0 i.e. in the cooling case
+            !! so similar to what's done in ECMWF:
+            zz1 = MAX(0._wp , zz1)    ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj))
+
+            zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases:
+            zz2 = zus*zus * roadrw
+            zz2 = zz2*zz2
+            zlamb =  6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders)
+
+            ! Updating molecular sublayer thickness (delta):
+            zz2    = rnu0_w/(SQRT(roadrw)*zus)
+            zdelta =      ztf    *          zlamb*zz2   &  ! Eq.12 (when alpha*Qb>0 / cooling of layer)
+               &    + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 )    ! Eq.12 (when alpha*Qb<0 / warming of layer)
+            !LB: changed 0.01 to 0.007
+
+            !! Once again with the new 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)
+            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
+
+            !! Update!
+            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
+            delta_vl(ji,jj) = zdelta
+
+         END DO
+      END DO
+
+   END SUBROUTINE CS_COARE
+
+
+
+
+   SUBROUTINE WL_COARE( kt,  pQsw, pQnsol, pTau, pSST, plon, isd, iwait,  pdT, &
+      &                    Hwl, mask_wl )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pTau*       surface wind stress                            [N/m^2]
+      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
+      !!     *plon*       longitude                                      [deg.E]
+      !!     *isd*        current UTC time, counted in second since 00h of the current day
+      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
+      !!
+      !!  **   OUTPUT:
+      !!     *pdT*   dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
+      !!---------------------------------------------------------------------
+      !!
+      !!   ** OPTIONAL OUTPUT:
+      !!     *Hwl*        depth of warm layer [m]
+      !!     *mask_wl*    mask for possible existence of a warm-layer (1) or not (0)
+      !!
+      !!------------------------------------------------------------------
+      INTEGER ,                     INTENT(in)  :: kt
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: plon     ! longitude [deg.E]
+      INTEGER ,                     INTENT(in)  :: isd      ! current UTC time, counted in second since 00h of the current day
+      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT      ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
+      !!
+      REAL(wp),   DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl     ! depth of warm layer [m]
+      INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0)
+      !
+      !
+      INTEGER :: ji,jj
+      !
+      REAL(wp) :: zdT_wl, zQabs, zfr, zdz
+      REAL(wp) :: zqac, ztac
+      REAL(wp) :: zalpha_w, zcd1, zcd2, flg
+      CHARACTER(len=128) :: cf_tmp
+      !!---------------------------------------------------------------------
+
+      REAL(wp) :: rlag_gw_h, &  ! local solar time lag in hours   / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai
+         &        rhr_sol       ! local solar time in hours since midnight
+
+      INTEGER  :: ilag_gw_s, &  ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai
+         &        isd_sol,   &  ! local solar time in number of seconds since local solar midnight
+         &        jl
+
+      !! INITIALIZATION:
+      pdT(:,:) = 0._wp    ! dT initially set to 0._wp
+      zdT_wl = 0._wp       ! total warming (amplitude) in warm layer
+      zQabs  = 0._wp       ! total heat absorped in warm layer
+      zfr    = zfr0        ! initial value of solar flux absorption
+      ztac   = 0._wp
+      zqac   = 0._wp
+      IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer
+
+            !! Need to know the local solar time from longitude and isd:
+            rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp )
+            rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp )
+            ilag_gw_s = INT( rlag_gw_h*3600._wp )
+            isd_sol = MODULO( isd + ilag_gw_s , 24*3600 )
+            rhr_sol = REAL( isd_sol , wp) / 3600._wp
+
+            !*****  variables for warm layer  ***
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w))        !mess-o-constants 1
+            zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
+
+            !********************************************************
+            !****  Compute apply warm layer  correction *************
+            !********************************************************
+            
+            !IF (isd_sol <= rdt ) THEN    !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime...
+            IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN
+               !PRINT *, '  [WL_COARE] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol
+               zdz           = H_wl_max
+               Tau_ac(ji,jj) = 0._wp
+               Qnt_ac(ji,jj) = 0._wp
+            END IF
+
+            IF ( rhr_sol > 5._wp ) THEN  ! ( 5am)
+
+               !PRINT *, '  [WL_COARE] WE DO WL !'
+               !PRINT *, '  [WL_COARE] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4)
+
+               !************************************
+               !****   set warm layer constants  ***
+               !************************************
+
+               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
+
+               !PRINT *, '  [WL_COARE] rdt,  pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4)
+
+               IF ( zQabs >= Qabs_thr ) THEN         ! Check for threshold
+
+                  !PRINT *, '  [WL_COARE] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4)
+
+                  !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! momentum integral
+                  ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
+
+                  IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN         !check threshold for warm layer existence
+                     !******************************************
+                     ! Compute the absorption profile
+                     !******************************************
+                     DO jl = 1, 5                           !loop 5 times for zfr
+                        zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) &
+                           &        + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz
+                        zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
+                        IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
+                     END DO
+
+                  ELSE
+                     !***********************
+                     ! Warm layer wiped out
+                     !***********************
+                     zfr  = 0.75
+                     zdz  = H_wl_max
+                     zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
+
+                  END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp )
+
+                  IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 !
+                  
+                  ! Warm layer correction
+                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz)
+                  pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
+                  
+               END IF ! IF ( zQabs >= Qabs_thr )
+
+            END IF ! IF ( isd_sol >= 21600 ) THEN  ! (21600 == 6am)
+
+            IF ( iwait == 0 ) THEN
+               IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN
+                  !PRINT *, '  [WL_COARE] WE UPDATE ACCUMULATED FLUXES !!!'
+                  Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
+                  Tau_ac(ji,jj) = ztac !
+                  IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1
+               END IF
+            END IF
+
+            H_wl(ji,jj) = zdz
+
+            IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj)
+
+         END DO
+      END DO
+      
+   END SUBROUTINE WL_COARE
+
+
+   !!======================================================================
+END MODULE sbcblk_skin_coare
Index: MO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare3p6.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare3p6.F90	(revision 11625)
+++ 	(revision )
@@ -1,301 +1,0 @@
-MODULE sbcblk_skin_coare3p6
-   !!======================================================================
-   !!                   ***  MODULE  sbcblk_skin_coare3p6  ***
-   !! Computes:
-   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
-   !!     scheme used at ECMWF
-   !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019)
-   !!
-   !!   From routine turb_ecmwf maintained and developed in AeroBulk
-   !!            (https://github.com/brodeau/aerobulk)
-   !!
-   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
-   !!----------------------------------------------------------------------
-   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
-   !!               based on the cool-skin/warm-layer scheme used at ECMWF
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-
-   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
-
-   USE lib_mpp         ! distribued memory computing library
-   USE in_out_manager  ! I/O manager
-   USE lib_fortran     ! to use key_nosignedzero
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC :: CS_COARE3P6, WL_COARE3P6
-
-   !! Cool-skin related parameters:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl  !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m]
-
-   !! Warm-layer related parameters:
-   REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp    !: maximum depth of warm layer (adjustable)
-   !
-   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number
-   !
-   REAL(wp), PARAMETER :: Qabs_thr = 50._wp  !: threshold for heat flux absorbed in WL
-   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac  ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac    ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl     ! depth of warm-layer [m]
-   !!----------------------------------------------------------------------
-CONTAINS
-
-
-   SUBROUTINE CS_COARE3P6( pQsw, pQnsol, pustar, pSST, pQlat,  pdT )
-      !!---------------------------------------------------------------------
-      !!
-      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
-      !!     ------------------------------------------------------------------
-      !!
-      !!  **   INPUT:
-      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
-      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
-      !!     *pustar*     friction velocity u*                           [m/s]
-      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
-      !!     *pQlat*      surface latent heat flux                       [K]
-      !!
-      !!  **  INPUT/OUTPUT:
-      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
-      !!              as output =>  new estimate of dT cool-skin
-      !!
-      !!------------------------------------------------------------------
-      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(inout) :: pdT    ! dT due to cool-skin effect
-      !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
-         &        zz1, zz2, zus, &
-         &        ztf
-      !!---------------------------------------------------------------------
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-
-            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
-
-            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
-
-            zdelta = delta_vl(ji,jj)   ! using last value of delta
-
-            !! Fraction of the shortwave flux absorbed by the cool-skin sublayer:
-            zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! 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 !
-
-            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) )  ! Total cooling at the interface
-
-            ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1
-            !                                 Qt < 0 => warming of the layer => ztf = 0
-
-            !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap):
-            zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap  ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0
-            !! LB: this terms only makes sense if > 0 i.e. in the cooling case
-            !! so similar to what's done in ECMWF:
-            zz1 = MAX(0._wp , zz1)    ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj))
-
-            zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases:
-            zz2 = zus*zus * roadrw
-            zz2 = zz2*zz2
-            zlamb =  6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders)
-
-            ! Updating molecular sublayer thickness (delta):
-            zz2    = rnu0_w/(SQRT(roadrw)*zus)
-            zdelta =      ztf    *          zlamb*zz2   &  ! Eq.12 (when alpha*Qb>0 / cooling of layer)
-               &    + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 )    ! Eq.12 (when alpha*Qb<0 / warming of layer)
-            !LB: changed 0.01 to 0.007
-
-            !! Once again with the new 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)
-            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
-
-            !! Update!
-            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
-            delta_vl(ji,jj) = zdelta
-
-         END DO
-      END DO
-
-   END SUBROUTINE CS_COARE3P6
-
-
-
-
-   SUBROUTINE WL_COARE3P6( kt,  pQsw, pQnsol, pTau, pSST, plon, isd, iwait,  pdT, &
-      &                    Hwl, mask_wl )
-      !!---------------------------------------------------------------------
-      !!
-      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
-      !!     ------------------------------------------------------------------
-      !!
-      !!  **   INPUT:
-      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
-      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
-      !!     *pTau*       surface wind stress                            [N/m^2]
-      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
-      !!     *plon*       longitude                                      [deg.E]
-      !!     *isd*        current UTC time, counted in second since 00h of the current day
-      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
-      !!
-      !!  **   OUTPUT:
-      !!     *pdT*   dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
-      !!---------------------------------------------------------------------
-      !!
-      !!   ** OPTIONAL OUTPUT:
-      !!     *Hwl*        depth of warm layer [m]
-      !!     *mask_wl*    mask for possible existence of a warm-layer (1) or not (0)
-      !!
-      !!------------------------------------------------------------------
-      INTEGER ,                     INTENT(in)  :: kt
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: plon     ! longitude [deg.E]
-      INTEGER ,                     INTENT(in)  :: isd      ! current UTC time, counted in second since 00h of the current day
-      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT      ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
-      !!
-      REAL(wp),   DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl     ! depth of warm layer [m]
-      INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0)
-      !
-      !
-      INTEGER :: ji,jj
-      !
-      REAL(wp) :: zdT_wl, zQabs, zfr, zdz
-      REAL(wp) :: zqac, ztac
-      REAL(wp) :: zalpha_w, zcd1, zcd2, flg
-      CHARACTER(len=128) :: cf_tmp
-      !!---------------------------------------------------------------------
-
-      REAL(wp) :: rlag_gw_h, &  ! local solar time lag in hours   / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai
-         &        rhr_sol       ! local solar time in hours since midnight
-
-      INTEGER  :: ilag_gw_s, &  ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai
-         &        isd_sol,   &  ! local solar time in number of seconds since local solar midnight
-         &        jl
-
-      !! INITIALIZATION:
-      pdT(:,:) = 0._wp    ! dT initially set to 0._wp
-      zdT_wl = 0._wp       ! total warming (amplitude) in warm layer
-      zQabs  = 0._wp       ! total heat absorped in warm layer
-      zfr    = zfr0        ! initial value of solar flux absorption
-      ztac   = 0._wp
-      zqac   = 0._wp
-      IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-
-            zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer
-
-            !! Need to know the local solar time from longitude and isd:
-            rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp )
-            rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp )
-            ilag_gw_s = INT( rlag_gw_h*3600._wp )
-            isd_sol = MODULO( isd + ilag_gw_s , 24*3600 )
-            rhr_sol = REAL( isd_sol , wp) / 3600._wp
-
-            !*****  variables for warm layer  ***
-            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
-
-            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w))        !mess-o-constants 1
-            zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
-
-            !********************************************************
-            !****  Compute apply warm layer  correction *************
-            !********************************************************
-            
-            !IF (isd_sol <= rdt ) THEN    !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime...
-            IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN
-               !PRINT *, '  [WL_COARE3P6] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol
-               zdz           = H_wl_max
-               Tau_ac(ji,jj) = 0._wp
-               Qnt_ac(ji,jj) = 0._wp
-            END IF
-
-            IF ( rhr_sol > 5._wp ) THEN  ! ( 5am)
-
-               !PRINT *, '  [WL_COARE3P6] WE DO WL !'
-               !PRINT *, '  [WL_COARE3P6] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4)
-
-               !************************************
-               !****   set warm layer constants  ***
-               !************************************
-
-               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
-
-               !PRINT *, '  [WL_COARE3P6] rdt,  pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4)
-
-               IF ( zQabs >= Qabs_thr ) THEN         ! Check for threshold
-
-                  !PRINT *, '  [WL_COARE3P6] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4)
-
-                  !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! momentum integral
-                  ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
-
-                  IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN         !check threshold for warm layer existence
-                     !******************************************
-                     ! Compute the absorption profile
-                     !******************************************
-                     DO jl = 1, 5                           !loop 5 times for zfr
-                        zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) &
-                           &        + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz
-                        zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
-                        IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
-                     END DO
-
-                  ELSE
-                     !***********************
-                     ! Warm layer wiped out
-                     !***********************
-                     zfr  = 0.75
-                     zdz  = H_wl_max
-                     zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
-
-                  END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp )
-
-                  IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 !
-                  
-                  ! Warm layer correction
-                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz)
-                  pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
-                  
-               END IF ! IF ( zQabs >= Qabs_thr )
-
-            END IF ! IF ( isd_sol >= 21600 ) THEN  ! (21600 == 6am)
-
-            IF ( iwait == 0 ) THEN
-               IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN
-                  !PRINT *, '  [WL_COARE3P6] WE UPDATE ACCUMULATED FLUXES !!!'
-                  Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
-                  Tau_ac(ji,jj) = ztac !
-                  IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1
-               END IF
-            END IF
-
-            H_wl(ji,jj) = zdz
-
-            IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj)
-
-         END DO
-      END DO
-      
-   END SUBROUTINE WL_COARE3P6
-
-
-   !!======================================================================
-END MODULE sbcblk_skin_coare3p6
