Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11209)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11221)
@@ -26,17 +26,16 @@
    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 iom             ! I/O manager library
-   USE lib_mpp         ! distribued memory computing library
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
    USE lib_fortran     ! to use key_nosignedzero
 
-   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
 
    IMPLICIT NONE
@@ -45,4 +44,6 @@
    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90
    
+   INTEGER , PARAMETER ::   nb_itt = 4        ! number of itterations
+
    !!----------------------------------------------------------------------
 CONTAINS
@@ -59,20 +60,4 @@
       !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
       !!
-      !! ** Method : Monin Obukhov Similarity Theory
-      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
-      !!
-      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
-      !!
-      !! ** Last update: Laurent Brodeau, June 2014:
-      !!    - handles both cases zt=zu and zt/=zu
-      !!    - optimized: less 2D arrays allocated and less operations
-      !!    - better first guess of stability by checking air-sea difference of virtual temperature
-      !!       rather than temperature difference only...
-      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
-      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
-      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
-      !!      => 'vkarmn' and 'grav'
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!
       !! INPUT :
@@ -95,4 +80,7 @@
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
       !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
@@ -111,7 +99,6 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
-      INTEGER ::   j_itt
-      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
-      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
+      INTEGER :: j_itt
+      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient
@@ -124,5 +111,5 @@
       !
       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 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
       U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
@@ -137,5 +124,5 @@
          ztmp0   (:,:) = cdn_wave(:,:)
       ELSE
-         ztmp0 = cd_neutral_10m( U_blk )
+      ztmp0 = cd_neutral_10m( U_blk )
       ENDIF
 
@@ -144,9 +131,12 @@
       !! Initializing transf. coeff. with their first guess neutral equivalents :
       Cd = ztmp0
-      Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 )
-      Ch = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
+      Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 )
+      Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))
       stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd)
  
-      IF( ln_cdgw )   Cen = Ce  ; Chn = Ch
+      IF( ln_cdgw ) THEN
+	Cen = Ce
+	Chn = Ch
+      ENDIF
 
       !! Initializing values at z_u with z_t values:
@@ -169,5 +159,5 @@
          !! Stability parameters :
          zeta_u   = zu*ztmp0
-         zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u )
+         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
          zpsi_h_u = psi_h( zeta_u )
 
@@ -176,5 +166,5 @@
             !! Array 'stab' is free for the moment so using it to store 'zeta_t'
             stab = zt*ztmp0
-            stab = SIGN( MIN(ABS(stab),10.0_wp), stab )  ! Temporaty array stab == zeta_t !!!
+            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!!
             stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
             t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
@@ -187,9 +177,9 @@
             stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
             Cd   = stab * stab
-            ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
+            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
-            ztmp1 = 1. + Chn * ztmp0     
+            ztmp1 = 1._wp + Chn * ztmp0     
             Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
-            ztmp1 = 1. + Cen * ztmp0
+            ztmp1 = 1._wp + Cen * ztmp0
             Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
 
@@ -205,20 +195,20 @@
 
          stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
-         Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
+         Cx_n10  = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
          Chn(:,:) = Cx_n10
 
          !! Update of transfer coefficients:
-         ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
+         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
          Cd      = ztmp0 / ( ztmp1*ztmp1 )
          stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
 
-         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
+         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
          ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
-         ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
+         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
          Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
 
-         Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
+         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
          Cen(:,:) = Cx_n10
-         ztmp1 = 1. + Cx_n10*ztmp0
+         ztmp1 = 1._wp + Cx_n10*ztmp0
          Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
          ENDIF
@@ -255,7 +245,7 @@
             zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
             !
-            cd_neutral_10m(ji,jj) = 1.e-3 * ( &
-               &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s
-               &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s
+            cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &
+               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
+               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
             !
             cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
