Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 13650)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 13655)
@@ -5,5 +5,5 @@
    !!   * bulk transfer coefficients C_D, C_E and C_H
    !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
-   !!   * the effective bulk wind speed at 10m U_blk
+   !!   * the effective bulk wind speed at 10m Ubzu
    !!   => all these are used in bulk formulas in sbcblk.F90
    !!
@@ -24,19 +24,9 @@
    !!                   returns the effective bulk wind speed at 10m
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
-   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 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 lib_mpp,        ONLY: ctl_stop         ! distribued memory computing library
+   USE in_out_manager, ONLY: nit000  ! I/O manager
+   USE sbc_phy         ! all thermodynamics functions, rho_air, q_sat, etc... !LB
    USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB
 
@@ -45,9 +35,7 @@
 
    PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
 
    !! ECMWF own values for given constants, taken form IFS documentation...
-   REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!!
+   REAL(wp), PARAMETER, PUBLIC :: charn0_ecmwf = 0.018_wp    ! Charnock constant (pretty high value here !!!
    !                                          !    =>  Usually 0.011 for moderate winds)
    REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1
@@ -57,5 +45,6 @@
    REAL(wp), PARAMETER ::   alpha_Q = 0.62    !
 
-   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
 
    !!----------------------------------------------------------------------
@@ -94,6 +83,6 @@
 
    SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
-      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           &
-      &                      Cdn, Chn, Cen,                                           &
+      &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                            &
+      &                      nb_iter, Cdn, Chn, Cen,                                           & ! optional output
       &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only
@@ -151,5 +140,5 @@
       !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
-      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
+      !!    *  Ubzu   : bulk wind speed at zu                                 [m/s]
       !!
       !!
@@ -171,7 +160,10 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
-      !
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
+      !
+      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
       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]
@@ -181,5 +173,5 @@
       REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
       !
-      INTEGER :: j_itt
+      INTEGER :: nbit, jit
       LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
       !
@@ -197,4 +189,7 @@
       !!----------------------------------------------------------------------------------
       IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
+
+      nbit = nb_iter0
+      IF( PRESENT(nb_iter) ) nbit = nb_iter
 
       l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
@@ -227,11 +222,11 @@
       znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
 
-      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
+      Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
 
       ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
       ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
-      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
-
-      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10)
+
+      z0     = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star
       z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
 
@@ -239,21 +234,19 @@
       z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
 
-      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
-
-      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
-
-      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
+      Cd     = MAX( (vkarmn/ztmp0)**2 , Cx_min )   ! first guess of Cd
+
+      ztmp0 = vkarmn2/LOG(zt/z0t)/Cd
+
+      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! 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 )
-      func_m = ztmp0*ztmp2 ! temporary array !!
-      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" !
-
+      func_h = (1._wp - ztmp1) *   ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & !  BRN < 0
+         &  +       ztmp1      * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 )                 !  BRN > 0
+      
       !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
       ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
 
-      u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on)
+      u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on)
       t_star = dt_zu*ztmp0
       q_star = dq_zu*ztmp0
@@ -276,8 +269,8 @@
 
 
-      !! First guess of inverse of Monin-Obukov length (1/L) :
+      !! First guess of inverse of Obukov length (1/L) :
       Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
 
-      !! Functions such as  u* = U_blk*vkarmn/func_m
+      !! Functions such as  u* = Ubzu*vkarmn/func_m
       ztmp0 = zu*Linv
       func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
@@ -285,10 +278,10 @@
 
       !! ITERATION BLOCK
-      DO j_itt = 1, nb_itt
+      DO jit = 1, nbit
 
          !! Bulk Richardson Number at z=zu (Eq. 3.25)
-         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
-
-         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
+         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN)
+
+         !! New estimate of the inverse of the Obukhon length (Linv == zeta/zu) :
          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
@@ -299,16 +292,16 @@
 
          !! Need to update roughness lengthes:
-         u_star = U_blk*vkarmn/func_m
+         u_star = Ubzu*vkarmn/func_m
          ztmp2  = u_star*u_star
          ztmp1  = znu_a/u_star
-         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
-         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
-         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
+         z0     = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp)
+         z0t    = MIN( ABS( alpha_H*ztmp1                           ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
+         z0q    = MIN( ABS( alpha_Q*ztmp1                           ) , 0.001_wp)
 
          !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8)
          ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
          !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
-         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
-         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
+         Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
+         ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0.
 
 
@@ -346,5 +339,5 @@
             !! Cool-skin contribution
 
-            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, &
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
                &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
 
@@ -359,5 +352,5 @@
          IF( l_use_wl ) THEN
             !! Warm-layer contribution
-            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, &
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
                &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
             CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
@@ -373,14 +366,14 @@
          ENDIF
 
-      END DO !DO j_itt = 1, nb_itt
-
-      Cd = vkarmn*vkarmn/(func_m*func_m)
-      Ch = vkarmn*vkarmn/(func_m*func_h)
-      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
-      Ce = vkarmn*vkarmn/(func_m*ztmp2)
-
-      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))
-      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))
-      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
+      END DO !DO jit = 1, nbit
+
+      Cd = MAX( vkarmn2/(func_m*func_m) , Cx_min )
+      Ch = MAX( vkarmn2/(func_m*func_h) , Cx_min )
+      ztmp2 = LOG(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
+      Ce = MAX( vkarmn2/(func_m*ztmp2)  , Cx_min )
+
+      IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min )
+      IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
+      IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0q)*LOG(zu/z0q)) , Cx_min )
 
       IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
@@ -408,29 +401,26 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
-      !!----------------------------------------------------------------------------------
+      REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc
+      !!----------------------------------------------------------------------------------
+      zc = 5._wp/0.35_wp
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-      !
-      zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
-      !
-      ! Unstable (Paulson 1970):
-      !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
-      zx = SQRT(ABS(1._wp - 16._wp*zzeta))
-      ztmp = 1._wp + SQRT(zx)
-      ztmp = ztmp*ztmp
-      psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
-         &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
-      !
-      ! Unstable:
-      ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
-         &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
-      !
-      ! Combining:
-      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
-      !
-      psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
-         &                +      stab  * psi_stab      ! (zzeta > 0) Stable
-      !
+            !
+            zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
+
+            ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
+            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
+            ztmp = 1._wp + zx
+            zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi
+
+            ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
+               &       - zta - 2._wp/3._wp*zc
+            !
+            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
+            !
+            psi_m_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
+               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
+            !
       END_2D
    END FUNCTION psi_m_ecmwf
@@ -452,27 +442,27 @@
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab
-      !!----------------------------------------------------------------------------------
+      REAL(wp) ::  zta, zx2, zpsi_unst, zpsi_stab, zstab, zc
+      !!----------------------------------------------------------------------------------
+      zc = 5._wp/0.35_wp
       !
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-      !
-      zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
-      !
-      zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
-      !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
-      ! Unstable (Paulson 1970) :
-      psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
-      !
-      ! Stable:
-      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-         &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
-      ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
-      !
-      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
-      !
-      !
-      psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
-         &                +    stab    * psi_stab        ! (zzeta > 0) Stable
-      !
+            !
+            zta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
+            !
+            ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
+            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
+            !
+            ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
+               &       - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp
+            !
+            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
+            !
+            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
+            !
+            psi_h_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
+               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable            
+            !
       END_2D
    END FUNCTION psi_h_ecmwf
