Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/cfgs/SHARED/namelist_ref	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/cfgs/SHARED/namelist_ref	(revision 12953)
@@ -677,4 +677,5 @@
    !
    ln_drgimp   = .true.    !  implicit top/bottom friction flag
+      ln_drgice_imp = .false. ! implicit ice-ocean drag
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/ICE/iceupdate.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/ICE/iceupdate.F90	(revision 12953)
@@ -26,4 +26,5 @@
    USE icectl         ! sea-ice: control prints
    USE bdy_oce , ONLY : ln_bdy
+   USE zdfdrg
    !
    USE in_out_manager ! I/O manager
@@ -323,4 +324,5 @@
       REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
       REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
+      REAL(wp) ::   zflagi                          !   -      -
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('ice_update_tau')
@@ -355,4 +357,12 @@
       !
       !                                      !==  every ocean time-step  ==!
+      IF ( ln_drgice_imp ) THEN
+         ! Save drag with right sign to update top drag in the ocean implicit
+         ! friction 
+         rCdU_ice(:,:) = -r1_rau0 * tmod_io(:,:) * at_i(:,:) * tmask(:,:,1) 
+         zflagi = 0._wp
+      ELSE
+         zflagi = 1._wp
+      ENDIF
       !
       DO jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle
@@ -364,6 +374,6 @@
                &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji  ,jj+1,1) )
             !                                                   ! linearized quadratic drag formulation
-            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
-            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - zflagi * pu_oce(ji,jj) )
+            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - zflagi * pv_oce(ji,jj) )
             !                                                   ! stresses at the ocean surface
             utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynnxt.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynnxt.F90	(revision 12953)
@@ -48,4 +48,5 @@
    USE prtctl         ! Print control
    USE timing         ! Timing
+   USE zdfdrg, ONLY: ln_drgice_imp, rCdU_top
 #if defined key_agrif
    USE agrif_oce_interp
@@ -99,4 +100,5 @@
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
@@ -354,4 +356,40 @@
       ENDIF
       !
+      IF ( iom_use("utau") ) THEN
+         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
+            ALLOCATE(zutau(jpi,jpj)) 
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  jk = miku(ji,jj) 
+                  zutau(ji,jj) = utau(ji,jj) & 
+                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk) 
+               END DO
+            END DO
+            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
+            CALL iom_put(  "utau", zutau(:,:) )
+            DEALLOCATE(zutau)
+         ELSE
+            CALL iom_put(  "utau", utau(:,:) )
+         ENDIF
+      ENDIF
+      !
+      IF ( iom_use("vtau") ) THEN
+         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
+            ALLOCATE(zvtau(jpi,jpj))
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  jk = mikv(ji,jj)
+                  zvtau(ji,jj) = vtau(ji,jj) &
+                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
+               END DO
+            END DO
+            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
+            CALL iom_put(  "vtau", zvtau(:,:) )
+            DEALLOCATE(zvtau)
+         ELSE
+            CALL iom_put(  "vtau", vtau(:,:) )
+         ENDIF
+      ENDIF
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynspg_ts.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynspg_ts.F90	(revision 12953)
@@ -1465,5 +1465,5 @@
       !                    !==  Set the barotropic drag coef.  ==!
       !
-      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
          
          DO jj = 2, jpjm1
@@ -1528,5 +1528,5 @@
       !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
       !
-      IF( ln_isfcav ) THEN
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynzdf.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynzdf.F90	(revision 12953)
@@ -141,5 +141,5 @@
             END DO
          END DO
-         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
+         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
             DO jj = 2, jpjm1        
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -258,5 +258,5 @@
             END DO
          END DO
-         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
+         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
@@ -423,5 +423,5 @@
             END DO
          END DO
-         IF ( ln_isfcav ) THEN
+         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/SBC/sbcmod.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/SBC/sbcmod.F90	(revision 12953)
@@ -564,7 +564,4 @@
       ENDIF
       !
-      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
-      CALL iom_put( "vtau", vtau )   ! j-wind stress
-      !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
          CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask )
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfdrg.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfdrg.F90	(revision 12953)
@@ -46,4 +46,5 @@
    LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0)
    LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag
+   LOGICAL , PUBLIC ::   ln_drgice_imp ! implicit ice-ocean drag
 
    !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist *
@@ -71,4 +72,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCd0_top, rCd0_bot   !: precomputed top/bottom drag coeff. at t-point (>0)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCdU_top, rCdU_bot   !: top/bottom drag coeff. at t-point (<0)  [m/s]
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCdU_ice             !: ice-ocean drag (used if ln_drgice_imp=T)
 
    !! * Substitutions
@@ -231,5 +233,5 @@
       INTEGER   ::   ios, ioptio   ! local integers
       !!
-      NAMELIST/namdrg/ ln_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp
+      NAMELIST/namdrg/ ln_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp
       !!----------------------------------------------------------------------
       !
@@ -254,4 +256,5 @@
          WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer
          WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp
+         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp
       ENDIF
       !
@@ -264,4 +267,6 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
       !
+      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) & 
+         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' )
       !
       !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
@@ -274,6 +279,13 @@
       !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities)
       !
-      IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting
-         ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) )
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting
+         ALLOCATE( rCdU_top(jpi,jpj) )
+         rCdU_top(:,:) = 0._wp
+      ENDIF
+      !
+      IF (ln_drgice_imp) ALLOCATE( rCdU_ice(jpi,jpj) )
+      !
+      IF( ln_isfcav ) THEN
+         ALLOCATE( rCd0_top(jpi,jpj))
          CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in
             &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out
Index: /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfphy.F90	(revision 12952)
+++ /NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/ZDF/zdfphy.F90	(revision 12953)
@@ -252,4 +252,12 @@
       ENDIF
       !
+      IF ( ln_drgice_imp) THEN
+         IF ( ln_isfcav ) THEN
+            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:)
+         ELSE
+            rCdU_top(:,:) = rCdU_ice(:,:)
+         ENDIF
+      ENDIF
+      ! 
       !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
       !
