Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90	(revision 12799)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90	(revision 13159)
@@ -632,5 +632,16 @@
 
       END SELECT
-
+      
+      IF( iom_use('Cd_oce') )   CALL iom_put("Cd_oce",   zcd_oce * tmask(:,:,1))
+      IF( iom_use('Ce_oce') )   CALL iom_put("Ce_oce",   zce_oce * tmask(:,:,1))
+      IF( iom_use('Ch_oce') )   CALL iom_put("Ch_oce",   zch_oce * tmask(:,:,1))
+      !! LB: mainly here for debugging purpose:
+      IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
+      IF( iom_use('q_zt') )     CALL iom_put("q_zt",     zqair       * tmask(:,:,1)) ! specific humidity       "
+      IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
+      IF( iom_use('q_zu') )     CALL iom_put("q_zu",     q_zu        * tmask(:,:,1)) ! specific humidity       "
+      IF( iom_use('ssq') )      CALL iom_put("ssq",      pssq        * tmask(:,:,1)) ! saturation specific humidity at z=0
+      IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu       * tmask(:,:,1)) ! bulk wind speed at z=zu
+      
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          !! ptsk and pssq have been updated!!!
@@ -643,14 +654,4 @@
          END WHERE
       END IF
-
-      !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef.
-      !!      CALL iom_put( "Ch_oce", zch_oce)  ! 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:
-         t_zu(:,:) = ztpot(:,:)
-         q_zu(:,:) = zqair(:,:)
-      ENDIF
-
 
       !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
@@ -671,5 +672,5 @@
             &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  &
             &               taum(:,:), psen(:,:), zqla(:,:),                   &
-            &               pEvap=pevp(:,:), prhoa=rhoa(:,:)                   )
+            &               pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
 
          zqla(:,:) = zqla(:,:) * tmask(:,:,1)
@@ -688,6 +689,6 @@
          ! ... utau, vtau at U- and V_points, resp.
          !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
-         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
-         DO_2D_10_10
+         !     Note that coastal wind stress is not used in the code... so this extra care has no effect
+         DO_2D_00_00
             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
@@ -893,23 +894,28 @@
          Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical
       ENDIF
-
-      !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef.
-      !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef.
-
+      
+      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice)
+      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice)
+      IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice)
+      
       ! local scalars ( place there for vector optimisation purposes)
       zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
 
       IF( ln_blk ) THEN
-         ! ------------------------------------------------------------ !
-         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
-         ! ------------------------------------------------------------ !
-         ! C-grid ice dynamics :   U & V-points (same as ocean)
-         DO_2D_00_00
-            putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             &
-               &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          &
-               &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) )
-            pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             &
-               &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          &
-               &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) )
+         ! ------------------------------------------------------------- !
+         !    Wind stress relative to the moving ice ( U10m - U_ice )    !
+         ! ------------------------------------------------------------- !
+         zztmp1 = rn_vfac * 0.5_wp
+         DO_2D_01_01    ! at T point 
+            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndi(ji,jj) - zztmp1 * ( puice(ji-1,jj  ) + puice(ji,jj) ) )
+            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndj(ji,jj) - zztmp1 * ( pvice(ji  ,jj-1) + pvice(ji,jj) ) )
+         END_2D
+         !
+         DO_2D_00_00    ! U & V-points (same as ocean).
+            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
+            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
+            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
+            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) )
+            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) )
          END_2D
          CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
@@ -1050,5 +1056,5 @@
       evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
       devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
-      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
+      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean  !LB: removed rn_efac here, correct???
 
       ! --- evaporation minus precipitation --- !
