Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_ice.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_ice.F90	(revision 13719)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_ice.F90	(revision 13806)
@@ -20,5 +20,5 @@
 # endif
 # if defined key_cice
-   USE ice_domain_size, only: ncat 
+   USE ice_domain_size, only: ncat
 #endif
    USE lib_mpp          ! MPP library
@@ -32,5 +32,5 @@
 # if defined  key_si3
    LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .TRUE.   !: SI3 ice model
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE 
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE
 # endif
 # if defined  key_cice
@@ -47,6 +47,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_ice_top    !: solar flux transmitted below the ice surface      [W/m2]
 
@@ -87,9 +87,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point
-   
+
    ! variables used in the coupled interface
    INTEGER , PUBLIC, PARAMETER ::   jpl = ncat
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice 
-   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice
+
    ! already defined in ice.F90 for SI3
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
@@ -98,6 +98,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K]
 #endif
-
-   !#LB: REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] !#LB => moved to sbc_phy.F90 !!!
 
    !! arrays relating to embedding ice in the ocean
@@ -145,5 +143,5 @@
          &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , &
          &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , &
-         &                     STAT= ierr(3) )      
+         &                     STAT= ierr(3) )
       IF( ln_cpl )   ALLOCATE( h_i(jpi,jpj,jpl) , h_s(jpi,jpj,jpl) , STAT=ierr(4) )
 #endif
@@ -168,6 +166,6 @@
    LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model
    LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE ice model
-   !#LB: REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-] !#LB => moved to sbc_phy.F90 !!!  
-   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
+
+   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
@@ -178,6 +176,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt, botmelt
    !
-   !! arrays related to embedding ice in the ocean. 
-   !! These arrays need to be declared even if no ice model is required. 
+   !! arrays related to embedding ice in the ocean.
+   !! These arrays need to be declared even if no ice model is required.
    !! In the no ice model or traditional levitating ice cases they contain only zeros
    !! ---------------------
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_oce.F90	(revision 13719)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_oce.F90	(revision 13806)
@@ -159,5 +159,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
 
-   !#LB:
    !!----------------------------------------------------------------------
    !!                     Surface atmospheric fields
@@ -165,5 +164,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_air_zt       !: specific humidity of air at z=zt [kg/kg]ww
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: theta_air_zt   !: potential temperature of air at z=zt [K]
-   !#LB.
 
    
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90	(revision 13719)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90	(revision 13806)
@@ -107,8 +107,8 @@
    !
    !#LB:
-   LOGICAL  ::   ln_Cx_ice_cst      ! use constant ice-air bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i)
-   REAL(wp) ::   rn_Cd_i, rn_Ce_i, rn_Ch_i
-   LOGICAL  ::   ln_Cx_ice_LU12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012)
-   LOGICAL  ::   ln_Cx_ice_LG15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
+   LOGICAL  ::   ln_Cx_ice_cst             ! use constant air-ice bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i)
+   REAL(wp) ::   rn_Cd_i, rn_Ce_i, rn_Ch_i ! values for  "    "
+   LOGICAL  ::   ln_Cx_ice_LU12            ! air-ice bulk transfer coefficients based on Lupkes et al., 2012)
+   LOGICAL  ::   ln_Cx_ice_LG15            ! air-ice bulk transfer coefficients based on Lupkes & Gryanik, 2015)
    !#LB.
    !
@@ -137,5 +137,5 @@
    LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
    LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB
-   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature
+   LOGICAL  ::   ln_tair_pot    ! temperature read in files ("sn_tair") is already potential temperature (not absolute)
    !
    INTEGER  ::   nhumi          ! choice of the bulk algorithm
@@ -218,12 +218,12 @@
       NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, &   ! bulk algorithm
          &                 rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl,       &
-         &                 rn_pfac, rn_efac,                                &
-         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                &   ! current feedback
-         &                 ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tpot,  &
-         &                 ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i,        &
-         &                 ln_Cx_ice_LU12, ln_Cx_ice_LG15,                  &
-         &                 cn_dir,                                          &
-         &                 sn_wndi, sn_wndj, sn_qsr, sn_qlw ,               &   ! input fields
-         &                 sn_tair, sn_humi, sn_prec, sn_snow, sn_slp,      &
+         &                 rn_pfac, rn_efac,                                          &
+         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback
+         &                 ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot,        &
+         &                 ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i,                  &
+         &                 ln_Cx_ice_LU12, ln_Cx_ice_LG15,                            &
+         &                 cn_dir,                                                    &
+         &                 sn_wndi, sn_wndj, sn_qsr, sn_qlw ,                         &   ! input fields
+         &                 sn_tair, sn_humi, sn_prec, sn_snow, sn_slp,                &
          &                 sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj
 
@@ -567,11 +567,11 @@
          !      based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
          !      (most reanalysis products provide absolute temp., not potential temp.)
-         IF( ln_tpot ) THEN
+         IF( ln_tair_pot ) THEN
+            ! temperature read into file is already potential temperature, do nothing...
+            theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)            
+         ELSE
             ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
-            IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing air POTENTIAL temperature out of ABSOLUTE temperature!'
+            IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
             theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
-         ELSE
-            ! temperature read into file is already potential temperature
-            theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
          ENDIF
          !
@@ -683,13 +683,13 @@
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-      zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
-      zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
-      ! ... scalar wind at T-point (not masked)
-      wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) )
+         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
+         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
+         ! ... scalar wind at T-point (not masked)
+         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) )
       END_2D
 #else
       ! ... scalar wind module at T-point (not masked)
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-      wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
+         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
       END_2D
 #endif
@@ -787,10 +787,10 @@
       IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         zztmp = zU_zu(ji,jj)
-         wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod
-         pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
-         psen(ji,jj)   = zztmp * zch_oce(ji,jj)
-         pevp(ji,jj)   = zztmp * zce_oce(ji,jj)
-         rhoa(ji,jj)   = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
+            zztmp = zU_zu(ji,jj)
+            wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod
+            pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
+            psen(ji,jj)   = zztmp * zch_oce(ji,jj)
+            pevp(ji,jj)   = zztmp * zce_oce(ji,jj)
+            rhoa(ji,jj)   = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
          END_2D
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
@@ -825,8 +825,8 @@
             zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
             DO_2D( 0, 1, 0, 1 )   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
-            zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax
-            ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) )
-            ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
-            taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
+               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax
+               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) )
+               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
+               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
             END_2D
          ENDIF
@@ -836,10 +836,9 @@
          !     Note that coastal wind stress is not used in the code... so this extra care has no effect
          DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T
-         utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) &
-            &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
-         vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) &
-            &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
+            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) &
+               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
+            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) &
+               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
          END_2D
-
 
          IF( ln_crt_fbk ) THEN
@@ -870,5 +869,5 @@
    END SUBROUTINE blk_oce_1
 
-   
+
    SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, &   ! <<= in
       &                   ptsk, psen, plat, pevp     )   ! <<= in
@@ -1339,6 +1338,6 @@
          DO jl = 1, jpl
             DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-            zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
-            IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
+               zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
+               IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
             END_2D
          END DO
@@ -1354,25 +1353,25 @@
       DO jl = 1, jpl
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         !
-         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) ) )
-         ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
-         ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
-         zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
-         !
-         DO iter = 1, nit     ! --- Iterative loop
-            zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
-            zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
-            ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
-         END DO
-         !
-         ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
-         qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
-         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 ) )
-
-         ! --- 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)
+            !
+            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) ) )
+            ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
+            ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
+            zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
+            !
+            DO iter = 1, nit     ! --- Iterative loop
+               zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
+               zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
+               ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
+            END DO
+            !
+            ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
+            qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
+            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 ) )
+
+            ! --- 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)
 
          END_2D
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 13719)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 13806)
@@ -160,8 +160,8 @@
       REAL(wp) :: zalfa     !: thermal expansion coefficient of sea-water [1/K]
       REAL(wp) :: zdTwl_b, zdTwl_n  !: temp. diff. between "almost surface (right below viscous layer) and bottom of WL
-      REAL(wp) :: zfr, zeta 
-      REAL(wp) :: zusw, zusw2 
-      REAL(wp) :: zLa, zfLa 
-      REAL(wp) :: flg, zwf, zQabs 
+      REAL(wp) :: zfr, zeta
+      REAL(wp) :: zusw, zusw2
+      REAL(wp) :: zLa, zfLa
+      REAL(wp) :: flg, zwf, zQabs
       REAL(wp) :: ZA, ZB, zL1, zL2
       REAL(wp) :: zcst0, zcst1, zcst2, zcst3
Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcmod.F90	(revision 13719)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcmod.F90	(revision 13806)
@@ -46,5 +46,5 @@
    USE sbcssr         ! surface boundary condition: sea surface restoring
    USE sbcrnf         ! surface boundary condition: runoffs
-   USE sbcapr         ! surface boundary condition: atmo pressure 
+   USE sbcapr         ! surface boundary condition: atmo pressure
    USE sbcfwb         ! surface boundary condition: freshwater budget
    USE icbstp         ! Icebergs
@@ -137,5 +137,5 @@
          WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
          WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl
-!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist 
+!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist
          WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
          WRITE(numout,*) '         components of your executable              nn_components = ', nn_components
@@ -162,5 +162,5 @@
       IF( .NOT.ln_wave ) THEN
          ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false.
-      ENDIF 
+      ENDIF
       IF( ln_sdw ) THEN
          IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) &
@@ -187,5 +187,5 @@
       !                       !**  check option consistency
       !
-      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS) 
+      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS)
       SELECT CASE( nn_components )
       CASE( jp_iam_nemo )
@@ -219,5 +219,5 @@
       SELECT CASE( nn_ice )
       CASE( 0 )                        !- no ice in the domain
-      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)  
+      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)
       CASE( 2 )                        !- SI3  ice model
          IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl .OR. ln_usr ) )   &
@@ -227,8 +227,8 @@
             &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl or ln_abl or ln_usr = T' )
          IF( lk_agrif                                )   &
-            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
+            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' )
       CASE DEFAULT                     !- not supported
       END SELECT
-      IF( ln_diurnal .AND. .NOT. ln_blk  )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )
+      IF( ln_diurnal .AND. .NOT. (ln_blk.OR.ln_abl) )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )
       !
       !                       !**  allocate and set required variables
@@ -242,5 +242,5 @@
       !
       IF( sbc_ssr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_ssr arrays' )
-      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring 
+      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring
          qrp(:,:) = 0._wp
          erp(:,:) = 0._wp
@@ -331,5 +331,5 @@
          &   CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
       !
-   
+
       !                       !**  associated modules : initialization
       !
@@ -395,5 +395,5 @@
       !
       REAL(wp) ::     zthscl        ! wd  tanh scale
-      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt  
+      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt
 
       !!---------------------------------------------------------------------
@@ -427,8 +427,8 @@
       !
       !                                            !==  sbc formulation  ==!
-      !                                                   
+      !
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
-      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                        ! user defined formulation 
+      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                        ! user defined formulation
       CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                             ! flux formulation
       CASE( jp_blk     )
@@ -447,5 +447,5 @@
       IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
       !
-      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
+      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves
       !
       !                                            !==  Misc. Options  ==!
@@ -461,6 +461,6 @@
       IF( ln_icebergs    )   THEN
                                      CALL icb_stp( kt )           ! compute icebergs
-         ! Icebergs do not melt over the haloes. 
-         ! So emp values over the haloes are no more consistent with the inner domain values. 
+         ! Icebergs do not melt over the haloes.
+         ! So emp values over the haloes are no more consistent with the inner domain values.
          ! A lbc_lnk is therefore needed to ensure reproducibility and restartability.
          ! see ticket #2113 for discussion about this lbc_lnk.
@@ -476,5 +476,5 @@
       ! Special treatment of freshwater fluxes over closed seas in the model domain
       ! Should not be run if ln_diurnal_only
-      IF( l_sbc_clo      )   CALL sbc_clo( kt )   
+      IF( l_sbc_clo      )   CALL sbc_clo( kt )
 
 !!$!RBbug do not understand why see ticket 667
@@ -482,5 +482,5 @@
 !!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
       IF( ll_wd ) THEN     ! If near WAD point limit the flux for now
-         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
+         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999
          zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
                                                      ! depth above wd limit once
