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 11209)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11215)
@@ -270,5 +270,5 @@
          WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
          WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
-         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
+         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p5
          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
          WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
@@ -472,6 +472,19 @@
       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+
+      !LB: Skin!!!         
+      CASE( np_COARE_3p0 )
+         IF ( ln_skin ) THEN            
+            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
+               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1))
+            zst(:,:) = zst(:,:)*tmask(:,:,1)
+            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
+         ELSE
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'            
+            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         END IF
+            
       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
@@ -480,17 +493,13 @@
       CASE( np_ECMWF     )
          IF ( ln_skin ) THEN            
-            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!'
-            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30)
             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1))
-            !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)...
             zst(:,:) = zst(:,:)*tmask(:,:,1)
             zsq(:,:) = zsq(:,:)*tmask(:,:,1)
-            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => AFTER ZST(40:50,30) =', zst(40:50,30)
          ELSE
             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
-               &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
          END IF
          !LB.
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11209)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11215)
@@ -13,5 +13,5 @@
    !!                     (https://github.com/brodeau/aerobulk)
    !!
-   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
@@ -26,15 +26,17 @@
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
-   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE iom             ! I/O manager library
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
    USE sbcwave, ONLY   :  cdn_wave ! wave module
 #if defined key_si3 || defined key_cice
    USE sbc_ice         ! Surface boundary condition: ice fields
 #endif
-   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-   USE lib_mpp        ! distribued memory computing library
-   USE prtctl         ! Print control
-   USE lib_fortran    ! to use key_nosignedzero
+   USE lib_fortran     ! to use key_nosignedzero
+
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
+   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
@@ -52,7 +54,8 @@
 CONTAINS
 
-   SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
+   SUBROUTINE turb_coare( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
-      &                   Cdn, Chn, Cen                       )
+      &                   Cdn, Chn, Cen,                      &
+      &                   Qsw, rad_lw, slp                   )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare  ***
@@ -63,4 +66,8 @@
       !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
       !!
+      !!                Applies the cool-skin warm-layer correction of the SST to T_s
+      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
+      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
+      !!                are provided as (optional) arguments!
       !!
       !! INPUT :
@@ -69,9 +76,19 @@
       !!    *  zu   : height for wind speed (generally 10m)                   [m]
       !!    *  U_zu : scalar wind speed at 10m                                [m/s]
-      !!    *  sst  : SST                                                     [K]
       !!    *  t_zt : potential air temperature at zt                         [K]
-      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
       !!
+      !! INPUT/OUTPUT:
+      !! -------------
+      !!    *  T_s  : SST or skin temperature                                 [K]
+      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
+      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
+      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
+      !!
+      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
+      !! ---------------
+      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
+      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
+      !!    *  slp    : sea-level pressure                                    [Pa]
       !!
       !! OUTPUT :
@@ -85,11 +102,11 @@
       !!
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
+      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
@@ -102,4 +119,8 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
+      !
       INTEGER :: j_itt
       LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
@@ -113,34 +134,50 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
-      !!----------------------------------------------------------------------
-      !
+      !
+      ! Cool skin:
+      LOGICAL :: l_use_skin = .FALSE.
+      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      !!----------------------------------------------------------------------------------
+      !
+      ! Cool skin ?
+      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
+         l_use_skin = .TRUE.
+      END IF
+      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
+
       l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
       IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
 
+      !! Initialization for cool skin:
+      IF( l_use_skin ) THEN
+         zsst   = T_s    ! save the bulk SST
+         T_s    = T_s - 0.25                      ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      END IF
+
       !! First guess of temperature and humidity at height zu:
-      t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions...
+      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
       q_zu = MAX( q_zt , 1.e-6_wp )   !               "
 
       !! Pot. temp. difference (and we don't want it to be 0!)
-      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
-      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
+      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
 
       znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
 
-      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
+      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution
       U_blk = SQRT(U_zu*U_zu + ztmp2)
 
-      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
+      ztmp2   = 10000._wp     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
       ztmp0   = LOG(zu*ztmp2)
       ztmp1   = LOG(10.*ztmp2)
-      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
-
-
-      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star
-      z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
-      z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
-      z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
+
+      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
+      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
       ztmp2  = vkarmn/ztmp0
@@ -149,13 +186,12 @@
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
 
-      ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
-      !ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
+      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
 
       !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
       ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       ztmp0 = ztmp0*ztmp2
-      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  BRN < 0
-         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  BRN > 0
-      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
+      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
+         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
+      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
       !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
@@ -174,8 +210,8 @@
          t_zu = t_zt - t_star/vkarmn*ztmp1
          q_zu = q_zt - q_star/vkarmn*ztmp1
-         q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
+         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
          !
-         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
-         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
+         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
       END IF
 
@@ -200,6 +236,6 @@
 
          !! Roughness lengthes z0, z0t (z0q = z0t) :
-         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
-         ztmp1 = z0*u_star/znu_a                           ! Re_r: roughness Reynolds number
+         z0    = ztmp2*ztmp1/grav + 0.11_wp*znu_a/u_star  ! Roughness length (eq.6)
+         ztmp1 = z0*u_star/znu_a                          ! Re_r: roughness Reynolds number
          z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
 
@@ -227,10 +263,30 @@
             t_zu = t_zt - t_star/vkarmn*ztmp1
             q_zu = q_zt - q_star/vkarmn*ztmp1
-            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
-            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
          END IF
 
-      END DO
-      !
+         !! SKIN related part
+         !! -----------------
+         IF( l_use_skin ) THEN
+            !! compute transfer coefficients at zu : lolo: verifier...
+            ztmp0 = u_star/U_blk
+            Ch   = ztmp0*t_star/dt_zu
+            Ce   = ztmp0*q_star/dq_zu
+            ! Non-Solar heat flux to the ocean:
+            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
+            ztmp2 = T_s*T_s
+            ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux
+               &       +    (rad_lw - emiss_w*stefan*ztmp2*ztmp2)                  ! Net longwave flux
+            !! Updating the values of the skin temperature T_s and q_s :
+            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
+            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
+         END IF
+
+         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
+            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
+         END IF
+
+      END DO !DO j_itt = 1, nb_itt
+
       ! compute transfer coefficients at zu :
       ztmp0 = u_star/U_blk
@@ -259,5 +315,5 @@
       !! Wind greater than 18 m/s :  alfa = 0.018
       !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
+      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
@@ -297,5 +353,5 @@
       !!       form from Beljaars and Holtslag (1991)
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
@@ -349,5 +405,5 @@
       !! form from Beljaars and Holtslag (1991)
       !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
+      !! Author: L. Brodeau, June 2016 / AeroBulk
       !!         (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11209)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11215)
@@ -14,5 +14,5 @@
    !!                     (https://github.com/brodeau/aerobulk)
    !!
-   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
@@ -38,6 +38,6 @@
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF)
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
+   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
@@ -57,19 +57,4 @@
    INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations
 
-   !! Cool-Skin / Warm-Layer related parameters:
-   REAL(wp),    PARAMETER :: &
-      &  rdt0    = 3600.*1.5, & !: time step
-      &  rd0     = 3. ,      &  !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005)
-      &  rNu0    = 0.5          !: Nu (exponent of temperature profile) Eq.11
-   !                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
-   !                            !: 0.3 to respect a warming of +3 K in calm
-   !                            !: condition for the insolation peak of +1000W/m^2
-   INTEGER,    PARAMETER :: &
-      &  nb_itt_wl = 10         !: number of sub-itterations for solving the differential equation in warm-layer part
-   !                            !:  => use "nb_itt_wl = 1" for No itteration! => way cheaper !!!
-   !                            !:    => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0)
-   !                            !:    => in that case no need for sub-itterations !
-   !                            !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions)
-   !                            !:  => Otherwize use "nb_itt_wl = 10"
    !!----------------------------------------------------------------------
 CONTAINS
@@ -123,5 +108,5 @@
       !!
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
@@ -157,10 +142,9 @@
       LOGICAL :: l_use_skin = .FALSE.
       REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
-
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
       !!----------------------------------------------------------------------------------
-      !
+
       ! Cool skin ?
       IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
@@ -169,5 +153,5 @@
       IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
 
-      ! Identical first gess as in COARE, with IFS parameter values though
+      ! Identical first gess as in COARE, with IFS parameter values though...
       !
       l_zt_equal_zu = .FALSE.
@@ -182,6 +166,6 @@
 
       !! First guess of temperature and humidity at height zu:
-      t_zu = MAX( t_zt , 0.0_wp  )   ! who knows what's given on masked-continental regions...
-      q_zu = MAX( q_zt , 1.e-6_wp)   !               "
+      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
+      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
 
       !! Pot. temp. difference (and we don't want it to be 0!)
@@ -189,5 +173,5 @@
       dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
 
-      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
+      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
 
       ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution
@@ -214,6 +198,6 @@
       ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       func_m = ztmp0*ztmp2 ! temporary array !!
-      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u
-         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))           !  BRN > 0
+      func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u
+         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0
       !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
@@ -259,5 +243,4 @@
          Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
          !! Note: it is slightly different that the L we would get with the usual
-         !! expression, as in coare algorithm or in 'mod_phy.f90' (One_on_L_MO())
          Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
@@ -315,5 +298,4 @@
          IF( l_use_skin ) THEN
             !! compute transfer coefficients at zu : lolo: verifier...
-            Cd = vkarmn*vkarmn/(func_m*func_m)
             Ch = vkarmn*vkarmn/(func_m*func_h)
             ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q
@@ -357,5 +339,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
@@ -405,5 +387,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
