Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90	(revision 12770)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90	(revision 12799)
@@ -18,5 +18,5 @@
    USE dom_oce, ONLY  : tmask
    USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa
-   USE sbcblk         ! use rn_?fac
+   USE sbcblk         ! use rn_?fac, cdn_oce
    USE sbcblk_phy     ! use some physical constants for flux computation
    !
@@ -102,10 +102,10 @@
 #endif
       !
-      REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zwnd_i, zwnd_j
-      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF
-      !
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c
+      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j
+      REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF
+      !
+      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_a
+      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_b
+      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_c
       !
       INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
@@ -351,5 +351,5 @@
 #if defined key_si3
             zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
-            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) * rn_vfac
+            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) ) * rn_vfac
             zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
 #endif
@@ -425,5 +425,5 @@
 #if defined key_si3
             zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
-            zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) * rn_vfac
+            zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac
             zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
 #endif
@@ -624,7 +624,7 @@
 
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl( tab2d_1=pwndm  , clinfo1=' abl_stp: wndm   : ' )
-         CALL prt_ctl( tab2d_1=ptaui  , clinfo1=' abl_stp: utau   : ' )
-         CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' )
+         CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=umask,   &
+            &          tab2d_2=ptauj , clinfo2='          vtau   : ', mask2=vmask )
+         CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm   : ' )
       ENDIF
 
@@ -647,6 +647,8 @@
       CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
       !
-      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
-         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
+      IF(sn_cfctl%l_prtctl) THEN
+         CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask,   &
+            &          tab2d_2=ptauj_ice , clinfo2='          vtau_ice : ', mask2=vmask )
+      END IF
 #endif
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -722,5 +724,5 @@
          ! Compute vertical shear
          DO jk = 2, jpkam1
-            DO ji = 1,jpi
+            DO ji = 1, jpi
                zcff        = 1.0_wp / e3w_abl( jk )**2
                zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
@@ -889,7 +891,7 @@
             DO jk = 2, jpkam1
                DO ji = 1, jpi
-                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )
+                  zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
                   mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  &
-                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   )
+                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
                END DO
             END DO
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90	(revision 12770)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90	(revision 12799)
@@ -209,4 +209,5 @@
       ! ABL timestep
       rDt_abl = nn_fsbc * rn_Dt
+      IF(lwp) WRITE(numout,*) ' ABL timestep = ', rDt_abl,' s'
 
       ! Check parameters for dynamics
@@ -294,6 +295,4 @@
       ENDIF
 
-      rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI)
-
    END SUBROUTINE sbc_abl_init
 
@@ -335,55 +334,58 @@
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
 
-      !!-------------------------------------------------------------------------------------------
-      !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
-      !!-------------------------------------------------------------------------------------------
-
-      CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in
-         &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in
-         &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in
-         &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in
-         &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out
-
-#if defined key_si3
-      CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
-         &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),    &   !   <<= in
-         &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
-         &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
-#endif
-
-      !!-------------------------------------------------------------------------------------------
-      !! 3 - Advance ABL variables from now (n) to after (n+1)
-      !!-------------------------------------------------------------------------------------------
-
-      CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
-         &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
-         &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
-         &              sf(jp_slp )%fnow(:,:,1),                            &   !   <<= in
-         &              sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),   &   !   <<= in
-         &              zcd_du, zsen, zevp,                                 &   !   <=> in/out
-         &              wndm, utau, vtau, taum                              &   !   =>> out
-#if defined key_si3
-         &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
-         &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
-         &            , utau_ice, vtau_ice                                  &   !   =>> out
-#endif
-         &                                                                  )
-      !!-------------------------------------------------------------------------------------------
-      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
-      !!                                                                time swap is done in abl_stp
-      !!-------------------------------------------------------------------------------------------
-
-      CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            &
-         &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
-         &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
-         &            tsk_m, zsen, zevp                                )
-
-      CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)
-      IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file
-
-#if defined key_si3
-      ! Avoid a USE abl in icesbc module
-      sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
-#endif
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+
+        !!-------------------------------------------------------------------------------------------
+        !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
+        !!-------------------------------------------------------------------------------------------
+  
+        CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in
+           &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in
+           &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in
+           &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in
+           &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out
+
+#if defined key_si3
+         CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
+            &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),    &   !   <<= in
+            &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
+            &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
+#endif
+
+         !!-------------------------------------------------------------------------------------------
+         !! 3 - Advance ABL variables from now (n) to after (n+1)
+         !!-------------------------------------------------------------------------------------------
+   
+         CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
+            &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
+            &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
+            &              sf(jp_slp )%fnow(:,:,1),                            &   !   <<= in
+            &              sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),   &   !   <<= in
+            &              zcd_du, zsen, zevp,                                 &   !   <=> in/out
+            &              wndm, utau, vtau, taum                              &   !   =>> out
+#if defined key_si3
+            &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
+            &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
+            &            , utau_ice, vtau_ice                                  &   !   =>> out
+#endif
+            &                                                                  )
+         !!-------------------------------------------------------------------------------------------
+         !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
+         !!                                                                time swap is done in abl_stp
+         !!-------------------------------------------------------------------------------------------
+
+         CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            &
+            &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
+            &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
+            &            tsk_m, zsen, zevp                                )
+   
+         CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)
+         IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file
+
+#if defined key_si3
+         ! Avoid a USE abl in icesbc module
+         sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
+#endif
+      END IF
 
    END SUBROUTINE sbc_abl
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom.F90	(revision 12770)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/IOM/iom.F90	(revision 12799)
@@ -1193,8 +1193,9 @@
                            &         'we accept this case, even if there is a possible mix-up between z and time dimension' )   
                      idmspc = idmspc - 1
-                  ELSE
-                     CALL ctl_stop( TRIM(clinfo), 'To keep iom lisibility, when reading a '//clrankpv//'D array,'         ,   &
-                        &                         'we do not accept data with '//cldmspc//' spatial dimensions',   &
-                        &                         'Use ncwa -a to suppress the unnecessary dimensions' )
+                  !!GS: possibility to read 3D ABL atmopsheric forcing and use 1st level to force BULK simulation
+                  !ELSE
+                  !   CALL ctl_stop( TRIM(clinfo), 'To keep iom lisibility, when reading a '//clrankpv//'D array,',   &
+                  !      &                         'we do not accept data with '//cldmspc//' spatial dimensions'  ,   &
+                  !      &                         'Use ncwa -a to suppress the unnecessary dimensions' )
                   ENDIF
             ENDIF
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 12770)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90	(revision 12799)
@@ -658,11 +658,11 @@
 
       IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
-         !! FL do we need this multiplication by tmask ... ???
          DO_2D_11_11
-            zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1)
+            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), phumi(ji,jj), pslp(ji,jj) )
          END_2D
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
