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 11841)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11962)
@@ -142,22 +142,9 @@
       !!-------------------------------------------------------------------
       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
-         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
+         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), tsk(jpi,jpj), STAT=sbc_blk_alloc )
       !
       CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
       IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
    END FUNCTION sbc_blk_alloc
-
-   !LB:
-   INTEGER FUNCTION sbc_blk_cswl_alloc()
-      !!-------------------------------------------------------------------
-      !!             ***  ROUTINE sbc_blk_cswl_alloc ***
-      !!-------------------------------------------------------------------
-      !WRITE(numout,*) '*** LB: allocating tsk!'
-      ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc )
-      !WRITE(numout,*) '*** LB: done!'
-      CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc )
-      IF( sbc_blk_cswl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' )
-   END FUNCTION sbc_blk_cswl_alloc
-   !LB.
 
 
@@ -222,9 +209,6 @@
       !LB:
       !                             !** initialization of the cool-skin / warm-layer parametrization
-      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
-         IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
-         !                       ! allocate array(s) for cool-skin/warm-layer param.
-         IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
-      END IF
+      IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) &
+         & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
       !
       ioptio = 0
@@ -322,8 +306,7 @@
          !
          WRITE(numout,*)
-         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs !LB
-         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl !LB
+         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs
+         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl
          !
-         !LB:
          WRITE(numout,*)
          SELECT CASE( nhumi )              !* Print the choice of air humidity
@@ -332,5 +315,4 @@
          CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]'
          END SELECT
-         !LB.
          !
       ENDIF
@@ -390,5 +372,5 @@
          ENDIF
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
-         !LB:
+
          SELECT CASE( nhumi )
          CASE( np_humi_sph )
@@ -399,5 +381,5 @@
             qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file
          END SELECT
-         !LB.
+
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
          sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
@@ -445,5 +427,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] !LB
+      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg]
       !!---------------------------------------------------------------------
       !
@@ -510,8 +492,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
+         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
+         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
       END SELECT
@@ -529,19 +511,16 @@
 
          CASE( np_COARE_3p0 )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm
-            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
+            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
+               &                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) )
 
          CASE( np_COARE_3p6 )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm
-            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
+            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
+               &                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) )
 
          CASE( np_ECMWF     )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm
-            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
+            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  &
+               &                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) )
 
@@ -566,24 +545,20 @@
       ELSE !IF ( ln_skin_cs .OR. ln_skin_wl )
 
-
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
             !
          CASE( np_NCAR      )
-            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2
+            CALL turb_ncar    (      rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         &
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_COARE_3p0 )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!' !LBrm
-            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0
+            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_COARE_3p6 )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!' !LBrm
-            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
+            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_ECMWF     )
-            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' !LBrm
-            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
+            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  &
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
@@ -596,5 +571,5 @@
       !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
       !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
-      
+
       IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
          !! If zu == zt, then ensuring once for all that:
@@ -602,31 +577,26 @@
          q_zu(:,:) = zqair(:,:)
       END IF
-      
+
 
       !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
       ! -------------------------------------------------------------
-      
-      CALL BULK_FORMULA( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:), &
-         &                 wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1), &
-         &                 taum(:,:), zqsb(:,:), zqla(:,:),               &
-         &                 pEvap=zevap(:,:), prhoa=rhoa(:,:) )
-      
+
+      CALL BULK_FORMULA( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), &
+         &               Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:),           &
+         &               wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1),   &
+         &               taum(:,:), zqsb(:,:), zqla(:,:),                 &
+         &               pEvap=zevap(:,:), prhoa=rhoa(:,:) )
+
       zqla(:,:)  =  zqla(:,:) * tmask(:,:,1)
       zqsb(:,:)  =  zqsb(:,:) * tmask(:,:,1)
       taum(:,:)  =  taum(:,:) * tmask(:,:,1)
       zevap(:,:) = zevap(:,:) * tmask(:,:,1)
-      
+
       ! Tau i and j component on T-grid points, using array "Cd_atm" as a temporary array...
       Cd_atm = 0._wp
       WHERE ( wndm > 0._wp ) Cd_atm = taum / wndm
       zwnd_i = Cd_atm * zwnd_i
-      zwnd_j = Cd_atm * zwnd_j      
-      !DO jj = 1, jpj             ! tau i and j component on T-grid points
-      !   DO ji = 1, jpi
-      !      zztmp = taum(ji,jj) / MAX( wndm(ji,jj) , 0.01_wp )
-      !      zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-      !      zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
-      !   END DO
-      !END DO
+      zwnd_j = Cd_atm * zwnd_j
+
       !                          ! add the HF tau contribution to the wind stress module
       IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
@@ -652,5 +622,6 @@
       ! ----------------------------------------------------------------------------- !
 
-      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
+      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
+      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
       zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
 
@@ -848,8 +819,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
+         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
+         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
       END SELECT
