Index: NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/OCE/SBC/sbcblk.F90	(revision 13185)
+++ NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/OCE/SBC/sbcblk.F90	(revision 13258)
@@ -74,17 +74,22 @@
 #endif
 
-   INTEGER , PUBLIC            ::   jpfld         ! maximum number of files to read
-   INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   ! index of 10m wind velocity (i-component) (m/s)    at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   ! index of 10m wind velocity (j-component) (m/s)    at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   ! index of 10m air temperature             (Kelvin)
-   INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   ! index of specific humidity               ( % )
-   INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   ! index of solar heat                      (W/m2)
-   INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   ! index of Long wave                       (W/m2)
-   INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   ! index of total precipitation (rain+snow) (Kg/m2/s)
-   INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   ! index of snow (solid prcipitation)       (kg/m2/s)
-   INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   ! index of sea level pressure              (Pa)
-   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
-
+   INTEGER , PUBLIC, PARAMETER ::   jp_wndi  =  1   ! index of 10m wind velocity (i-component) (m/s)    at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin)
+   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % )
+   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2)
+   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2)
+   INTEGER , PUBLIC, PARAMETER ::   jp_prec  =  7   ! index of total precipitation (rain+snow) (Kg/m2/s)
+   INTEGER , PUBLIC, PARAMETER ::   jp_snow  =  8   ! index of snow (solid prcipitation)       (kg/m2/s)
+   INTEGER , PUBLIC, PARAMETER ::   jp_slp   =  9   ! index of sea level pressure              (Pa)
+   INTEGER , PUBLIC, PARAMETER ::   jp_uoatm = 10   ! index of surface current (i-component)
+   !                                                !          seen by the atmospheric forcing (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component)
+   !                                                !          seen by the atmospheric forcing (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read
+
+   ! Warning: keep this structure allocatable for Agrif...
    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read)
 
@@ -98,13 +103,16 @@
    LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
    !
+   LOGICAL  ::   ln_crt_fbk     ! Add surface current feedback to the wind stress computation  (Renault et al. 2020)
+   REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta
+   REAL(wp) ::   rn_stau_b      ! 
+   !
    REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation
    REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation
-   REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress
    REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements
    REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements
    !
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL)
+   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
+   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
 
    LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
@@ -113,4 +121,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
    !
    INTEGER  ::   nhumi          ! choice of the bulk algorithm
@@ -162,13 +171,17 @@
       !!
       CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
-      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i        ! array of namelist informations on the fields to read
-      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
-      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
-      TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        "
+      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
+      TYPE(FLD_N) ::   sn_wndi, sn_wndj , sn_humi, sn_qsr      ! informations about the fields to be read
+      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        "
+      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        "
+      TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        "
+      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
-         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       &
+         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     &
+         &                 sn_hpgi, sn_hpgj,                                          &
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , rn_zqt, rn_zu,                                    &
-         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
+         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot,           &
+         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback
          &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB
       !!---------------------------------------------------------------------
@@ -242,16 +255,17 @@
       !                                   !* set the bulk structure
       !                                      !- store namelist information in an array
-      IF( ln_blk ) jpfld = 9
-      IF( ln_abl ) jpfld = 11
-      ALLOCATE( slf_i(jpfld) )
-      !
-      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
-      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
-      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
-      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
-      slf_i(jp_slp ) = sn_slp
-      IF( ln_abl ) THEN
-         slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj
-      END IF
+      !
+      slf_i(jp_wndi ) = sn_wndi    ;   slf_i(jp_wndj ) = sn_wndj
+      slf_i(jp_qsr  ) = sn_qsr     ;   slf_i(jp_qlw  ) = sn_qlw
+      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi
+      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow
+      slf_i(jp_slp  ) = sn_slp
+      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm
+      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj
+      !
+      IF( .NOT. ln_abl ) THEN   ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know...
+         slf_i(jp_hpgi)%clname = 'NOT USED'
+         slf_i(jp_hpgj)%clname = 'NOT USED'
+      ENDIF
       !
       !                                      !- allocate the bulk structure
@@ -264,17 +278,27 @@
       DO jfpr= 1, jpfld
          !
-         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero)
-            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
-            sf(jfpr)%fnow(:,:,1) = 0._wp
+         IF(   ln_abl    .AND.                                                      &
+            &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   &
+            &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN
+            ipka = jpka   ! ABL: some fields are 3D input
+         ELSE
+            ipka = 1
+         ENDIF
+         !
+         ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) )
+         !
+         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default)
+            IF(     jfpr == jp_slp  ) THEN
+               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa
+            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN
+               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents
+            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN
+               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case
+            ELSE
+               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
+               CALL ctl_stop( ctmp1 )
+            ENDIF
          ELSE                                                  !-- used field  --!
-            IF(   ln_abl    .AND.                                                      &
-               &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   &
-               &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input
-               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) )
-               IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) )
-            ELSE                                                                                ! others or Bulk fields are 2D fiels
-               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
-               IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) )
-            ENDIF
+            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation
             !
             IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   &
@@ -327,8 +351,13 @@
          WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
          WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
-         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
          WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
          WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
          WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
+         WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk
+         IF(ln_crt_fbk) THEN
+         WRITE(numout,*) '         Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta'
+         WRITE(numout,*) '            Alpha                                         rn_stau_a    = ', rn_stau_a
+         WRITE(numout,*) '            Beta                                          rn_stau_b    = ', rn_stau_b
+         ENDIF
          !
          WRITE(numout,*)
@@ -429,14 +458,15 @@
       !                                            ! compute the surface ocean fluxes using bulk formulea
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
-         CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in
-            &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in
-            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in
-            &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs)
-            &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out
-
-         CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in
-            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in
-            &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in
-            &                zsen, zevp )                                            !   <=> in out
+         CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m,        &   !   <<= in
+            &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in (wl/cs)
+            &                tsk_m, zssq, zcd_du, zsen, zevp )                         !   =>> out
+
+         CALL blk_oce_2(     sf(jp_tair )%fnow(:,:,1), sf(jp_qsr  )%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in
+            &                zsen, zevp )                                              !   <=> in out
       ENDIF
       !
@@ -470,8 +500,8 @@
 
 
-   SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp
-      &                  pslp , pst   , pu   , pv,        &  ! inp
-      &                  pqsr , pqlw  ,                   &  ! inp
-      &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out
+   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi,         &  ! inp
+      &                      pslp , pst  , pu   , pv,            &  ! inp
+      &                      puatm, pvatm, pqsr , pqlw ,         &  ! inp
+      &                      ptsk , pssq , pcd_du, psen, pevp   )   ! out
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_oce_1  ***
@@ -498,4 +528,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s]
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   puatm  ! surface current seen by the atm at T-point (i-component) [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pvatm  ! surface current seen by the atm at T-point (j-component) [m/s]
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   !
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   !
@@ -508,5 +540,9 @@
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zztmp                ! local variable
+      REAL(wp) ::   zstmax, zstau
+#if defined key_cyclone
       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
+#endif
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
@@ -532,18 +568,16 @@
       zwnd_j(:,:) = 0._wp
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
-      DO_2D_00_00
-         pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
-         pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
+      DO_2D_11_11
+         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_11_11
+         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
       END_2D
 #endif
-      DO_2D_00_00
-         zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-         zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
-      END_2D
-      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
-      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
-      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
-         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
-
       ! ----------------------------------------------------------------------------- !
       !      I   Solar FLUX                                                           !
@@ -593,10 +627,12 @@
          !#LB: because AGRIF hates functions that return something else than a scalar, need to
          !     use scalar version of gamma_moist() ...
-         DO_2D_11_11
-            ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
-         END_2D
-      ENDIF
-
-
+         IF( ln_tpot ) THEN
+            DO_2D_11_11
+               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
+            END_2D
+         ELSE
+            ztpot = ptair(:,:)
+         ENDIF
+      ENDIF
 
       !! Time to call the user-selected bulk parameterization for
@@ -674,22 +710,47 @@
          pevp(:,:) = pevp(:,:) * tmask(:,:,1)
 
-         ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array...
-         zcd_oce = 0._wp
-         WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm
-         zwnd_i = zcd_oce * zwnd_i
-         zwnd_j = zcd_oce * zwnd_j
-
-         CALL iom_put( "taum_oce", taum )   ! output wind stress module
+         DO_2D_11_11
+            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)
+#else
+               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
+         END_2D
+
+         IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
+            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_01_01   ! 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) )
+            END_2D
+         ENDIF
 
          ! ... 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 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))
-            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
-               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
+         DO_2D_00_00              ! 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))
          END_2D
-         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+
+         IF( ln_crt_fbk ) THEN
+            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
+         ELSE
+            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+         ENDIF
+
+         CALL iom_put( "taum_oce", taum )   ! output wind stress module
 
          IF(sn_cfctl%l_prtctl) THEN
@@ -862,5 +923,4 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
       REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature
       REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays
@@ -873,10 +933,7 @@
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO_2D_00_00
-         zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  )
-         zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  )
-         wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+      DO_2D_11_11
+         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
       END_2D
-      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
       !
       ! Make ice-atm. drag dependent on ice concentration
@@ -898,11 +955,11 @@
 
       IF( ln_blk ) THEN
-         ! ------------------------------------------------------------- !
-         !    Wind stress relative to the moving ice ( U10m - U_ice )    !
-         ! ------------------------------------------------------------- !
-         zztmp1 = rn_vfac * 0.5_wp
+         ! ---------------------------------------------------- !
+         !    Wind stress relative to nonmoving ice ( U10m )    !
+         ! ---------------------------------------------------- !
+         ! supress moving ice in wind stress computation as we don't know how to do it properly...
          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) ) )
+            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj)
+            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj)
          END_2D
          !
@@ -914,9 +971,9 @@
             pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) )
          END_2D
-         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
+         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
          !
          IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
             &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
-      ELSE
+      ELSE ! ln_abl
          zztmp1 = 11637800.0_wp
          zztmp2 =    -5897.8_wp
@@ -1382,5 +1439,5 @@
          !
       END_2D
-      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1., pch, 'T', 1. )
+      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1.0_wp, pch, 'T', 1.0_wp )
       !
    END SUBROUTINE Cdn10_Lupkes2015
