Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbc_phy.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbc_phy.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbc_phy.F90	(revision 15551)
@@ -257,6 +257,8 @@
         zmask = 1._wp
         ztpot = ptpot
+        zta   = 0._wp
       ELSE
         zmask = 0._wp 
+        ztpot = 0._wp
         zta   = pta
       ENDIF
@@ -585,5 +587,5 @@
       REAL(wp) :: ztptv
       !!-------------------------------------------------------------------
-      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
+      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
       !
       zsstv = virt_temp_sclr( psst, pssq )   ! virtual potential SST
@@ -610,13 +612,15 @@
       INTEGER  ::   ji, jj
 
-      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         IF( l_ptqa_l_prvd ) THEN  !!GS: "IF" inside loop needs to be removed
+      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
+      IF( l_ptqa_l_prvd ) THEN
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
                &                                pta_layer=pta_layer(ji,jj ),  pqa_layer=pqa_layer(ji,jj ) )
-         ELSE
+         END_2D
+      ELSE
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
-         END IF
-      END_2D
+         END_2D
+      END IF
 
    END FUNCTION Ri_bulk_vctr
Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk.F90	(revision 15551)
@@ -275,6 +275,6 @@
          IF( ln_ANDREAS )      &
             & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' )
-         IF( nn_fsbc /= 1 ) &
-            & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
+         !IF( nn_fsbc /= 1 ) &
+         !   & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
       END IF
 
@@ -555,5 +555,6 @@
          ! Potential temperature of air at z=rn_zqt (most reanalysis products provide absolute temp., not potential temp.)
          IF( ln_tair_pot ) THEN
-            theta_air_zt(:,:)        = sf(jp_tair )%fnow(:,:,1)
+            ! 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...)
@@ -584,6 +585,5 @@
             qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
          ENDIF
-         tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)    !#LB: should it be POTENTIAL temperature instead ????
-         !tatm_ice(:,:) = theta_air_zt(:,:)         !#LB: THIS! ?
+         tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)    !#LB: should it be POTENTIAL temperature (theta_air_zt) instead ????
          qatm_ice(:,:) = q_air_zt(:,:)
          tprecip(:,:)  = sf(jp_prec)%fnow(:,:,1) * rn_pfac
@@ -649,5 +649,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean
       REAL(wp), DIMENSION(jpi,jpj) ::   zsspt             ! potential sea-surface temperature [K]
-      REAL(wp), DIMENSION(jpi,jpj) ::   zpre, ztabs
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpre, ztabs       ! air pressure [Pa] & absolute temperature [K]
       REAL(wp), DIMENSION(jpi,jpj) ::   zztmp1, zztmp2
       !!---------------------------------------------------------------------
@@ -762,5 +762,5 @@
 
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
-         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
+         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
          WHERE ( fr_i(:,:) > 0.001_wp )
             ! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
@@ -768,5 +768,5 @@
             pssq(:,:)  = zztmp2(:,:)
          END WHERE
-         !! ptsk and pssq have been updated!!!
+         ! apply potential temperature increment to abolute SST
          ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
       END IF
@@ -789,10 +789,8 @@
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
 
-         zpre(:,:) = pres_temp( q_zu(:,:), pslp(:,:), rn_zu, ptpot=theta_zu(:,:), pta=ztabs(:,:) )
-         rhoa(:,:) = rho_air( ztabs(:,:), q_zu(:,:), zpre(:,:) )
-         !!GS: for debug, to be removed 
-         CALL iom_put( "pres_zu", zpre(:,:)*tmask(:,:,1) )
-         CALL iom_put( "slp", pslp(:,:)*tmask(:,:,1) )
-         CALL iom_put( "tabs_zu", (ztabs(:,:)-rt0)*tmask(:,:,1) )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            zpre(ji,jj) = pres_temp( q_zu(ji,jj), pslp(ji,jj), rn_zu, ptpot=theta_zu(ji,jj), pta=ztabs(ji,jj) )
+            rhoa(ji,jj) = rho_air( ztabs(ji,jj), q_zu(ji,jj), zpre(ji,jj) )
+         END_2D
 
          CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
@@ -808,17 +806,17 @@
 
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         IF( wndm(ji,jj) > 0._wp ) THEN
-            zztmp = taum(ji,jj) / wndm(ji,jj)
+            IF( wndm(ji,jj) > 0._wp ) THEN
+              zztmp = taum(ji,jj) / wndm(ji,jj)
 #if defined key_cyclone
-            ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-            ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
 #else
-            ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
-            ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
+               ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
+               ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
 #endif
-         ELSE
-            ztau_i(ji,jj) = 0._wp
-            ztau_j(ji,jj) = 0._wp
-         ENDIF
+            ELSE
+               ztau_i(ji,jj) = 0._wp
+               ztau_j(ji,jj) = 0._wp
+            ENDIF
          END_2D
 
@@ -867,8 +865,7 @@
       ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius
 
-      CALL iom_put( "sspt", zsspt(:,:)-rt0 )
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius
-         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference...
+         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference
       ENDIF
       !
@@ -984,4 +981,5 @@
       !! ** Method  :   compute momentum using bulk formulation
       !!                formulea, ice variables and read atmospheric fields.
+      !!                NB: ice drag coefficient is assumed to be a constant
       !!---------------------------------------------------------------------
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa]
@@ -1016,5 +1014,5 @@
       ! potential sea-ice surface temperature [K]
       zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
-      
+
       ! sea-ice <-> atmosphere bulk transfer coefficients
       SELECT CASE( nblk_ice )
@@ -1046,6 +1044,7 @@
       END SELECT
 
-      ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
-      IF( iom_use('spt_ice')) CALL iom_put("spt_ice", (zsipt-rt0)*ztmp)
+      IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) &
+         & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
+
       IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
       IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
@@ -1131,4 +1130,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       REAL(wp), DIMENSION(jpi,jpj)     ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
@@ -1136,7 +1136,7 @@
       !
       zcoef_dqlw = 4._wp * emiss_i * stefan             ! local scalars
-      !
       zztmp = 1. / ( 1. - albo )
       dqla_ice(:,:,:) = 0._wp
+
       ! Heat content per unit mass (J/kg)
       zcptrain(:,:) = (      ptair        - rt0 ) * rcp  * tmask(:,:,1)
Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 15551)
@@ -189,5 +189,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
-      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta
+      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
       !
       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt
@@ -323,5 +323,5 @@
 
          IF(( l_use_cs ).OR.( l_use_wl )) THEN
-            zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+            zpre(:,:)  = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
             zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
          ENDIF
Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 15551)
@@ -179,5 +179,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
-      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta
+      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
       !
       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt
Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 15551)
@@ -184,5 +184,5 @@
       REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length...
       REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
-      REAL(wp), DIMENSION(jpi,jpj) :: zrhoa, zpre, zta
+      REAL(wp), DIMENSION(jpi,jpj) :: zrhoa, zpre, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
       !
       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst  ! to back up the initial bulk SST
@@ -356,5 +356,4 @@
             !! Cool-skin contribution
 
-
             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
                &                   ztmp1, ztmp0, Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
Index: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbccpl.F90	(revision 15550)
+++ NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbccpl.F90	(revision 15551)
@@ -54,5 +54,5 @@
 #endif
 
-   USE sbc_phy, ONLY : pp_cldf
+   USE sbc_phy, ONLY : pp_cldf, rpref
 
    IMPLICIT NONE
@@ -220,7 +220,4 @@
 #endif
 
-   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rho0)
-
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
 
@@ -1189,4 +1186,5 @@
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
+      REAL(wp) ::   r1_grau                ! = 1.e0 / (grav * rho0)
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
       !!----------------------------------------------------------------------
