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 11216)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11217)
@@ -17,5 +17,5 @@
    !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
-   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle) 
+   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)
    !!----------------------------------------------------------------------
 
@@ -24,10 +24,10 @@
    !!   sbc_blk       : bulk formulation as ocean surface boundary condition
    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
-   !!             sea-ice case only : 
+   !!             sea-ice case only :
    !!   blk_ice_tau   : provide the air-ice stress
    !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
    !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
    !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
-   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -45,8 +45,8 @@
    USE icethd_dh      ! for CALL ice_thd_snwblow
 #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_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 sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
    !
    USE iom            ! I/O manager library
@@ -58,5 +58,5 @@
    USE sbcblk_phy     !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc...
 
-   
+
    IMPLICIT NONE
    PRIVATE
@@ -68,5 +68,5 @@
    PUBLIC   blk_ice_flx   ! routine called in icesbc
    PUBLIC   blk_ice_qcn   ! routine called in icesbc
-#endif 
+#endif
 
    INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
@@ -96,5 +96,5 @@
    REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
-!!gm ref namelist initialize it so remove the setting to false below
+   !!gm ref namelist initialize it so remove the setting to false below
    LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
    LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
@@ -111,5 +111,5 @@
    LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB
    !LB.
-   
+
    INTEGER  ::   nblk           ! choice of the bulk algorithm
    !                            ! associated indices:
@@ -159,5 +159,5 @@
       !! ** Purpose :   choose and initialize a bulk formulae formulation
       !!
-      !! ** Method  : 
+      !! ** Method  :
       !!
       !!----------------------------------------------------------------------
@@ -173,5 +173,5 @@
          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
-         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
+         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
          &                 ln_skin, ln_humi_dpt                                                        ! cool-skin / warm-layer param. !LB
@@ -181,5 +181,5 @@
       IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
       !
-      !                             !** read bulk namelist  
+      !                             !** read bulk namelist
       REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
       READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
@@ -194,9 +194,17 @@
       !                             !** initialization of the chosen bulk formulae (+ check)
       !                                   !* select the bulk chosen in the namelist and check the choice
-                                                               ioptio = 0
-      IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
-      IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
-      IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
-      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
+      ioptio = 0
+      IF( ln_NCAR      ) THEN
+         nblk =  np_NCAR        ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ln_COARE_3p0 ) THEN
+         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ln_COARE_3p5 ) THEN
+         nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ln_ECMWF     ) THEN
+         nblk =  np_ECMWF       ;   ioptio = ioptio + 1
+      ENDIF
       !
       !LB:
@@ -208,13 +216,13 @@
       END IF
       !LB.
-         
+
       IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
       !
       IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
          IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
-         IF( sn_qsr%ln_tint ) THEN 
+         IF( sn_qsr%ln_tint ) THEN
             CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
                &           '              ==> We force time interpolation = .false. for qsr' )
-            sn_qsr%ln_tint = .false.
+            sn_qsr%ln_tint = .FALSE.
          ENDIF
       ENDIF
@@ -245,26 +253,26 @@
       !
       IF ( ln_wave ) THEN
-      !Activated wave module but neither drag nor stokes drift activated
+         !Activated wave module but neither drag nor stokes drift activated
          IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
             CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
-      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
+            !drag coefficient read from wave model definable only with mfs bulk formulae and core
+         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN
+            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
          ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
-             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
          ENDIF
       ELSE
-      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
-         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
-         &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
-         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
-         &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &  
-         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
-      ENDIF 
-      !
-      !           
+         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
+            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
+            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
+            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
+            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &
+            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
+      ENDIF
+      !
+      !
       IF(lwp) THEN                     !** Control print
          !
-         WRITE(numout,*)                  !* namelist 
+         WRITE(numout,*)                  !* namelist
          WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
          WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
@@ -344,7 +352,9 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
          qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
-         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
-         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
-         ENDIF 
+         IF( ln_dm2dc ) THEN
+            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
+         ELSE
+            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
+         ENDIF
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          !LB:
@@ -411,5 +421,5 @@
 
       ! ... components ( U10m - U_oce ) at T-point (unmasked)
-!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
+      !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
       zwnd_i(:,:) = 0._wp
       zwnd_j(:,:) = 0._wp
@@ -440,6 +450,8 @@
       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
       zztmp = 1. - albo
-      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
-      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
+      IF( ln_dm2dc ) THEN
+         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
+      ELSE
+         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
       ENDIF
 
@@ -451,4 +463,9 @@
       ! ... specific humidity at SST and IST tmask(
       zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
+
+      IF ( ln_skin ) THEN
+         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
+         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
+      END IF
 
       !LB:
@@ -458,5 +475,5 @@
          IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-      ELSE         
+      ELSE
          zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity!
       END IF
@@ -467,49 +484,66 @@
       !!    (since reanalysis products provide T at z, not theta !)
       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * 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, zqair, wndm,   &  ! NCAR-COREv2
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-
-      !LB: Skin!!!         
-      CASE( np_COARE_3p0 )
-         IF ( ln_skin ) THEN            
+      
+
+      IF ( ln_skin ) THEN
+         
+         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
+
+         CASE( np_COARE_3p0 )
             CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
                &                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     )
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF
+               &                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 DEFAULT
+            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' )
+         END SELECT
+
+         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq:
+         WHERE ( fr_i < 0.001_wp )
+            ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep it!!!
             zst(:,:) = zst(:,:)*tmask(:,:,1)
             zsq(:,:) = zsq(:,:)*tmask(:,:,1)
-         ELSE
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'            
+         ELSEWHERE
+            ! we forget about the update...
+            zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo
+            zsq(:,:) = zqla(:,:) !LB:  "   "   "
+         END WHERE
+
+         !LB: Update of tsk, the official array for skin temperature
+         tsk(:,:) = zst(:,:)
+
+      ELSE !IF ( ln_skin )
+
+         
+         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
+            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+
+         CASE( np_COARE_3p0 )
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'
             CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-         END IF
-            
-      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-
-      !LB: Skin!!!
-      CASE( np_ECMWF     )
-         IF ( ln_skin ) THEN            
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
-               &                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))
-            zst(:,:) = zst(:,:)*tmask(:,:,1)
-            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
-         ELSE
+
+         CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5
+            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+
+         CASE( np_ECMWF     )
             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-         END IF
-         !LB.
-         
-      CASE DEFAULT
-         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
-      END SELECT
-
-
-      !LB: cool-skin/warm-layer:
-      IF ( ln_skin ) tsk(:,:) = zst(:,:)
+
+         CASE DEFAULT
+            CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
+         END SELECT
+
+      END IF ! IF ( ln_skin )
+
+
 
 
@@ -521,6 +555,6 @@
       END IF
 
-!!      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.
+      !!      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.
 
       DO jj = 1, jpj             ! tau module, i and j component
@@ -612,20 +646,20 @@
       !
       !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
-         CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB
-         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
-         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
-         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
-         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
-         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
-         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
-         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
-         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
-         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
-         CALL iom_put( 'snowpre', sprecip )                 ! Snow
-         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
-         IF ( ln_skin ) THEN
-            CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
-            CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
-         END IF
+      CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB
+      CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
+      CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
+      CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
+      CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
+      CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
+      CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
+      CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
+      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
+      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
+      CALL iom_put( 'snowpre', sprecip )                 ! Snow
+      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
+      IF ( ln_skin ) THEN
+         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
+         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
+      END IF
       !!#LB. ENDIF
       !
@@ -650,5 +684,5 @@
    !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
    !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
-   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
    !!----------------------------------------------------------------------
 
@@ -693,9 +727,9 @@
          Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
       ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
-         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
-      ENDIF
-
-!!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
-!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
+         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
+      ENDIF
+
+      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
+      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
 
       ! local scalars ( place there for vector optimisation purposes)
@@ -776,6 +810,8 @@
 
       zztmp = 1. / ( 1. - albo )
-      WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
-      ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp
+      WHERE( ptsu(:,:,:) /= 0._wp )
+         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
+      ELSEWHERE
+         z1_st(:,:,:) = 0._wp
       END WHERE
       !                                     ! ========================== !
@@ -866,5 +902,5 @@
       DO jl = 1, jpl
          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
-         !                         ! But we do not have Tice => consider it at 0degC => evap=0 
+         !                         ! But we do not have Tice => consider it at 0degC => evap=0
       END DO
 
@@ -873,10 +909,10 @@
       zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
       !
-      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
       ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
       ELSEWHERE                                                         ! zero when hs>0
-         qtr_ice_top(:,:,:) = 0._wp 
+         qtr_ice_top(:,:,:) = 0._wp
       END WHERE
       !
@@ -891,5 +927,5 @@
       !
    END SUBROUTINE blk_ice_flx
-   
+
 
    SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
@@ -900,5 +936,5 @@
       !!                to force sea ice / snow thermodynamics
       !!                in the case conduction flux is emulated
-      !!                
+      !!
       !! ** Method  :   compute surface energy balance assuming neglecting heat storage
       !!                following the 0-layer Semtner (1976) approach
@@ -925,5 +961,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
       !!---------------------------------------------------------------------
-      
+
       ! -------------------------------------!
       !      I   Enhanced conduction factor  !
@@ -933,5 +969,5 @@
       !
       zgfac(:,:,:) = 1._wp
-      
+
       IF( ld_virtual_itd ) THEN
          !
@@ -939,6 +975,6 @@
          zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
          zfac3 = 2._wp / zepsilon
-         !   
-         DO jl = 1, jpl                
+         !
+         DO jl = 1, jpl
             DO jj = 1 , jpj
                DO ji = 1, jpi
@@ -948,7 +984,7 @@
             END DO
          END DO
-         !      
-      ENDIF
-      
+         !
+      ENDIF
+
       ! -------------------------------------------------------------!
       !      II   Surface temperature and conduction flux            !
@@ -960,5 +996,5 @@
          DO jj = 1 , jpj
             DO ji = 1, jpi
-               !                    
+               !
                zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
                   &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
@@ -977,16 +1013,16 @@
                qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
                qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  &
-                             &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
+                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
 
                ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
-               hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
+               hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)
 
             END DO
          END DO
          !
-      END DO 
-      !      
+      END DO
+      !
    END SUBROUTINE blk_ice_qcn
-   
+
 
    SUBROUTINE Cdn10_Lupkes2012( Cd )
@@ -994,9 +1030,9 @@
       !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
       !!
-      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
+      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
       !!                 to make it dependent on edges at leads, melt ponds and flows.
       !!                 After some approximations, this can be resumed to a dependency
       !!                 on ice concentration.
-      !!                
+      !!
       !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
       !!                 with the highest level of approximation: level4, eq.(59)
@@ -1010,5 +1046,5 @@
       !!
       !!                 This new drag has a parabolic shape (as a function of A) starting at
-      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
+      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
       !!                 and going down to Cdi(say 1.4e-3) for A=1
       !!
@@ -1037,5 +1073,5 @@
       Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
          &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
-      
+
    END SUBROUTINE Cdn10_Lupkes2012
 
@@ -1046,12 +1082,12 @@
       !!
       !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
-      !!                 between sea-ice and atmosphere with distinct momentum 
-      !!                 and heat coefficients depending on sea-ice concentration 
+      !!                 between sea-ice and atmosphere with distinct momentum
+      !!                 and heat coefficients depending on sea-ice concentration
       !!                 and atmospheric stability (no meltponds effect for now).
-      !!                
+      !!
       !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
       !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
       !!                 it considers specific skin and form drags (Andreas et al. 2010)
-      !!                 to compute neutral transfert coefficients for both heat and 
+      !!                 to compute neutral transfert coefficients for both heat and
       !!                 momemtum fluxes. Atmospheric stability effect on transfert
       !!                 coefficient is also taken into account following Louis (1979).
@@ -1093,8 +1129,10 @@
 
       ! mean temperature
-      WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
-      ELSEWHERE                       ;   ztm_su(:,:) = rt0
+      WHERE( at_i_b(:,:) > 1.e-20 )
+         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
+      ELSEWHERE
+         ztm_su(:,:) = rt0
       ENDWHERE
-      
+
       ! Momentum Neutral Transfert Coefficients (should be a constant)
       zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
@@ -1105,5 +1143,5 @@
       ! Heat Neutral Transfert Coefficients
       zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details)
-     
+
       ! Atmospheric and Surface Variables
       zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
@@ -1117,13 +1155,13 @@
             zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
             zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
-            
+
             ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
             zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
             zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
-            
+
             ! Momentum and Heat Neutral Transfert Coefficients
             zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
-            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
-                       
+            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
+
             ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
             z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
@@ -1137,5 +1175,5 @@
                zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
             ENDIF
-            
+
             IF( zrib_i <= 0._wp ) THEN
                zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
@@ -1145,9 +1183,9 @@
                zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
             ENDIF
-            
+
             ! Momentum Transfert Coefficients (Eq. 38)
             Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
                &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
-            
+
             ! Heat Transfert Coefficients (Eq. 49)
             Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
