Index: NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 12582)
+++ NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 12623)
@@ -98,5 +98,5 @@
       &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
       !!
@@ -184,7 +184,7 @@
       LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::  u_star, t_star, q_star
-      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu     
-      REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air
+      REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
+      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
+      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air
       REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length...
       REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
@@ -196,9 +196,7 @@
       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
       !!----------------------------------------------------------------------------------
-
       IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
 
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
 
       !! Initializations for cool skin and warm layer:
@@ -413,26 +411,26 @@
       !!----------------------------------------------------------------------------------
       DO_2D_11_11
-         !
-         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
-         !
+      !
+      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
+      !
       END_2D
    END FUNCTION psi_m_ecmwf
@@ -458,23 +456,23 @@
       !
       DO_2D_11_11
-         !
-         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
-         !
+      !
+      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
+      !
       END_2D
    END FUNCTION psi_h_ecmwf
