Index: /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 6726)
+++ /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 6727)
@@ -9,5 +9,5 @@
    !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier) additions:
    !!                           -  new bulk routine for efficiency
-   !!                           -  WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !!!!
+   !!                           -  WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
    !!                           -  file names and file characteristics in namelist
    !!                           -  Implement reading of 6-hourly fields
@@ -18,17 +18,16 @@
    !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE
    !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk
-   !!
-   !!            4.0  !  2016-06  (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to
-   !!                             the CORE algorithm anymore
-   !!                             => based on AeroBulk (http://aerobulk.sourceforge.net/)
+   !!            4.0  !  2016-06  (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
+   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   sbc_blk         : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
-   !!   blk_oce         : computes momentum, heat and freshwater fluxes over ocean
-   !!   blk_ice         : computes momentum, heat and freshwater fluxes over ice
-   !!   rho_air         : density of (moist) air (depends on T_air, q_air and SLP
-   !!   cp_air          : specific heat of (moist) air (depends spec. hum. q_air)
-   !!    q_sat          : saturation humidity as a unction of SLP and temperature
+   !!   sbc_blk       : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
+   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
+   !!   blk_ice       : computes momentum, heat and freshwater fluxes over ice
+   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
+   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
+   !!   q_sat         : saturation humidity as a function of SLP and temperature
+   !!   L_vap         : latent heat of vaporization of water as a function of temperature
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -49,10 +48,9 @@
    USE par_ice_2      ! LIM-2 parameters
 #endif
-
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
    USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003) 
    USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
    USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31) 
-
+   !
    USE iom            ! I/O manager library
    USE in_out_manager ! I/O manager
@@ -66,11 +64,12 @@
    PRIVATE
 
-   PUBLIC   sbc_blk         ! routine called in sbcmod module
+   PUBLIC   sbc_blk       ! routine called in sbcmod module
 #if defined key_lim2 || defined key_lim3
-   PUBLIC   blk_ice_tau     ! routine called in sbc_ice_lim module
-   PUBLIC   blk_ice_flx     ! routine called in sbc_ice_lim module
+   PUBLIC   blk_ice_tau   ! routine called in sbc_ice_lim module
+   PUBLIC   blk_ice_flx   ! routine called in sbc_ice_lim module
 #endif
 
-   !!Lolo: should ultimately be moved in the module with all physical constants ?
+!!Lolo: should ultimately be moved in the module with all physical constants ?
+!!gm  : In principle, yes.
    REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
    REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
@@ -381,45 +380,35 @@
 
       ! ... specific humidity at SST and IST tmask(
-      zsq(:,:) = 0.98 * q_sat(zst(:,:), sf(jp_slp)%fnow(:,:,1))
+      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
       !!
       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
       !!    (since reanalysis products provide T at z, not theta !)
-      ztpot =  sf(jp_tair)%fnow(:,:,1) + gamma_moist(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1)) * rn_zqt
-
-            SELECT CASE ( nblk )
-            CASE( np_NCAR      )   ;   WRITE(numout,*) '         "NCAR" algorithm        (Large and Yeager 2008)'
-            CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '         "COARE 3.0" algorithm   (Fairall et al. 2003)'
-            CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '         "COARE 3.5" algorithm   (Edson et al. 2013)'
-            CASE( np_ECMWF     )   ;   WRITE(numout,*) '         "ECMWF" algorithm       (IFS cycle 31)'
-            END SELECT
+      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
 
       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
-         !
-      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &
-            &                                            Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
-      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &
-            &                                            Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
-      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &
-            &                                            Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
-      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &
-            &                                            Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
+      !
+      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
+         &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
+      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
+         &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
+      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
+         &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
+      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
+         &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
       CASE DEFAULT
          CALL ctl_stop( 'STOP', 'sbc_blk: non-existing bulk formula selected' )
       END SELECT
 
-      ! Computing true air density :
-      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN
-         !! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
-         zrhoa = rho_air(zt_zu(:,:), zq_zu(:,:), sf(jp_slp)%fnow(:,:,1))
-      ELSE
-         !! At zt:
-         zrhoa = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
+      !                          ! Compute true air density :
+      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
+         zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) )
+      ELSE                                      ! At zt:
+         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       END IF
 
-      ! ... tau module, i and j component
-      DO jj = 1, jpj
+      DO jj = 1, jpj             ! tau module, i and j component
          DO ji = 1, jpi
-            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd(ji,jj) ! using bulk wind speed
+            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd(ji,jj)   ! using bulk wind speed
             taum  (ji,jj) = zztmp * wndm  (ji,jj)
             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
@@ -428,8 +417,7 @@
       END DO
 
-      ! ... add the HF tau contribution to the wind stress module?
-      IF( lhftau ) THEN
-         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
-      ENDIF
+      !                          ! add the HF tau contribution to the wind stress module
+      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
+
       CALL iom_put( "taum_oce", taum )   ! output wind stress module
 
@@ -466,5 +454,5 @@
       ENDIF
 
-      zqla(:,:) = Lvap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
+      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
 
 
@@ -641,17 +629,16 @@
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies)
-      !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      REAL(wp) ::   zst2, zst3
-      REAL(wp) ::   zcoef_dqlw, zcoef_dqla
-      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable
-      !!
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3)
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
+      !!
+      INTEGER  ::   ji, jj, jl               ! dummy loop indices
+      REAL(wp) ::   zst2, zst3               ! local variable
+      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
+      REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw        ! long wave heat sensitivity over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb        ! sensible  heat sensitivity over ice
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa
       !!---------------------------------------------------------------------
@@ -666,5 +653,5 @@
       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
       !
-      zrhoa(:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
+      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       !
       zztmp = 1. / ( 1. - albo )
@@ -782,20 +769,22 @@
 
    END SUBROUTINE blk_ice_flx
+   
 #endif
-
 
    FUNCTION rho_air( ptak, pqa, pslp )
       !!-------------------------------------------------------------------------------
-      !! ** Purpose : compute density of (moist) air with eq. of state
+      !!                           ***  FUNCTION rho_air  ***
+      !!
+      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
       !!
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
       !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak, &  !: air temperature   [K]
-         &                                        pqa, &  !: air spec. hum.    [kg/kg]
-         &                                        pslp    !: pressure in       [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   !:              [kg/m^3]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
       !!-------------------------------------------------------------------------------
       !
-      rho_air = pslp/(R_dry*ptak*(1._wp + rctv0*pqa))
+      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
       !
    END FUNCTION rho_air
@@ -804,13 +793,15 @@
    FUNCTION cp_air( pqa )
       !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
       !! ** Purpose : Compute specific heat (Cp) of moist air
       !!
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
       !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa     !: air spec. hum. [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             :: cp_air  !: [J/K/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
       !!-------------------------------------------------------------------------------
       !
-      Cp_air = Cp_dry + Cp_vap*pqa
+      Cp_air = Cp_dry + Cp_vap * pqa
       !
    END FUNCTION cp_air
@@ -819,4 +810,6 @@
    FUNCTION q_sat( ptak, pslp )
       !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION q_sat  ***
+      !!
       !! ** Purpose : Specific humidity at saturation in [kg/kg]
       !!              Based on accurate estimate of "e_sat"
@@ -825,11 +818,10 @@
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
       !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: &
-         &                            ptak,       &  !: air temperature [K]
-         &                            pslp           !: sea level atmospheric pressure  [Pa]
-      REAL(wp), DIMENSION(jpi,jpj) :: q_sat
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
       !
       INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: ze_sat, ztmp     ! local scalar
+      REAL(wp) ::   ze_sat, ztmp   ! local scalar
       !!----------------------------------------------------------------------------------
       !
@@ -837,13 +829,12 @@
          DO ji = 1, jpi
             !
-            ztmp = rt0/ptak(ji,jj)
+            ztmp = rt0 / ptak(ji,jj)
             !
-            ! Vapour pressure at saturation [hPa] :
-            !   WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)       &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) &
+            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
+            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
+               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
-            !
-            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )     ! 0.01 because SLP is in [Pa]
+               !
+            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
             !
          END DO
@@ -855,4 +846,6 @@
    FUNCTION gamma_moist( ptak, pqa )
       !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION gamma_moist  ***
+      !!
       !! ** Purpose : Compute the moist adiabatic lapse-rate.
       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
@@ -861,6 +854,7 @@
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
       !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             :: gamma_moist
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
       !
       INTEGER  ::   ji, jj         ! dummy loop indices
@@ -871,5 +865,5 @@
          DO ji = 1, jpi
             zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
-            ziRT = 1./(R_dry*ptak(ji,jj))    ! 1/RT
+            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
             gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
          END DO
@@ -879,17 +873,19 @@
 
 
-   FUNCTION Lvap( psst )
+   FUNCTION L_vap( psst )
       !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap  ***
+      !!
       !! ** Purpose : Compute the latent heat of vaporization of water from temperature
       !!
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
       !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             :: Lvap   !: [J/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst   !: water temperature [K]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
       !!----------------------------------------------------------------------------------
       !
-      Lvap = (2.501 - 0.00237*(psst(:,:) - rt0))*1.E6
-      !
-   END FUNCTION Lvap
+      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
+      !
+   END FUNCTION L_vap
 
    !!======================================================================
Index: /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90
===================================================================
--- /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90	(revision 6726)
+++ /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90	(revision 6727)
@@ -18,7 +18,7 @@
    !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com
    !!
-   !!====================================================================================
-
-   !! History :  3.6  !  2016-02  (L.Brodeau) Original code
+   !!======================================================================
+   !! History :  3.6  !  2016-02  (L.Brodeau)   Original code
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -27,26 +27,26 @@
    !!                   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 wrk_nemo        ! work arrays
-   USE timing          ! Timing
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE sbcwave, ONLY   :  cdn_wave ! wave module
+   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 sbcwave, ONLY  :  cdn_wave ! wave module
 #if defined key_lim3 || defined key_cice
-   USE sbc_ice         ! Surface boundary condition: ice fields
+   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 in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! distribued memory computing library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE prtctl         ! Print control
+   USE lib_fortran    ! to use key_nosignedzero
+
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC :: TURB_COARE
+   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
 
    !! COARE own values for given constants:
@@ -71,5 +71,5 @@
       !!
       !! ** Method : Monin Obukhov Similarity Theory
-      !!======================================================================================
+      !!----------------------------------------------------------------------
       !!
       !! INPUT :
@@ -92,23 +92,19 @@
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
       !!    *  U_blk  : bulk wind at 10m                                      [m/s]
-      !!
-      !!
-      !!============================================================================
-
-
-      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(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(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu            [m/s]
+      !!----------------------------------------------------------------------
+      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(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(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]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
-      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) ::   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 at 10m                          [m/s]
-
+      !
       INTEGER :: j_itt
       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
@@ -123,9 +119,10 @@
       REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt
       REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2
-
+      !!----------------------------------------------------------------------
+      !
       IF( nn_timing == 1 )  CALL timing_start('turb_coare')
 
-      CALL wrk_alloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu)
-      CALL wrk_alloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
+      CALL wrk_alloc( jpi,jpj,   u_star, t_star, q_star, zeta_u, dt_zu, dq_zu)
+      CALL wrk_alloc( jpi,jpj,   znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
 
       l_zt_equal_zu = .FALSE.
@@ -136,5 +133,5 @@
       !! First guess of temperature and humidity at height zu:
       t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
-      q_zu = MAX(q_zt , 1.E-6)  !               "
+      q_zu = MAX(q_zt , 1.e-6)  !               "
 
       !! Pot. temp. difference (and we don't want it to be 0!)
@@ -258,6 +255,4 @@
 
 
-
-
    FUNCTION alfa_charn( pwnd )
       !!-------------------------------------------------------------------
@@ -270,6 +265,5 @@
       !! Wind greater than 18 m/s :  alfa = 0.018
       !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://sourceforge.net/p/aerobulk)
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
@@ -278,4 +272,5 @@
       INTEGER  ::   ji, jj         ! dummy loop indices
       REAL(wp) :: zw, zgt10, zgt18
+      !!-------------------------------------------------------------------
       !
       DO jj = 1, jpj
@@ -296,5 +291,4 @@
       !
    END FUNCTION alfa_charn
-
 
 
@@ -315,4 +309,5 @@
       INTEGER  ::   ji, jj         ! dummy loop indices
       REAL(wp) ::     zqa          ! local scalar
+      !!-------------------------------------------------------------------
       !
       DO jj = 1, jpj
Index: /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90
===================================================================
--- /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90	(revision 6726)
+++ /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90	(revision 6727)
@@ -47,5 +47,5 @@
    PRIVATE
 
-   PUBLIC ::   TURB_COARE3P5
+   PUBLIC ::   TURB_COARE3P5   ! called by sbcblk.F90
 
    !                                   ! COARE own values for given constants:
Index: /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90	(revision 6726)
+++ /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90	(revision 6727)
@@ -19,7 +19,7 @@
    !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
    !!----------------------------------------------------------------------
-   !!
-   !! History :  4.0  !  2016-06  (L.Brodeau) Original code
-   !!
+   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   turb_ecmwf  : computes the bulk turbulent transfer coefficients
@@ -47,5 +47,5 @@
    PRIVATE
 
-   PUBLIC :: TURB_ECMWF
+   PUBLIC ::   TURB_ECMWF   ! called by sbcblk.F90
 
    !                   !! ECMWF own values for given constants, taken form IFS documentation...
Index: /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90
===================================================================
--- /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 6726)
+++ /branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 6727)
@@ -15,6 +15,5 @@
    !!                         L. Brodeau, 2015
    !!=====================================================================
-   !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar
-   !!                                         of old sbcblk_core.F90
+   !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90
    !!----------------------------------------------------------------------
 
@@ -45,7 +44,7 @@
    PRIVATE
 
-   PUBLIC :: TURB_NCAR
-
-   !! NCAR own values for given constants:
+   PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90
+
+   !                              ! NCAR own values for given constants:
    REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature...
    
