Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/cfgs/SHARED/field_def_nemo-oce.xml	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/cfgs/SHARED/field_def_nemo-oce.xml	(revision 15197)
@@ -856,4 +856,8 @@
      <field id="utrd_tfr"       long_name="i-trend: top friction (3D) in top layer"         unit="m/s^2"                        />   
      <field id="utrd_tfr2d"     long_name="i-trend: top friction: 2D component"             unit="m/s^2"    grid_ref="grid_U_2D" />
+     <field id="utrd_iceoc_iceimp"  long_name="i-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_U_2D" />
+     <field id="utrd_tfre_bt_iceimp"  long_name="i-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_U_2D" />
+     <field id="utrd_tfri_iceimp"  long_name="i-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_U_2D" />
+     <field id="utrd_tfr_iceimp"  long_name="i-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_U_2D" />
      <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        />   
      <field id="utrd_tot2d"     long_name="i-trend: total momentum trend before atf: 2D component"    unit="m/s^2"  grid_ref="grid_U_2D" />   
@@ -903,4 +907,8 @@
      <field id="vtrd_tfr"       long_name="j-trend: top friction (3D) in top layer"         unit="m/s^2"                        />   
      <field id="vtrd_tfr2d"     long_name="j-trend: top friction: 2D component"             unit="m/s^2"    grid_ref="grid_V_2D" />
+     <field id="vtrd_iceoc_iceimp"  long_name="j-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_V_2D" />
+     <field id="vtrd_tfre_bt_iceimp"  long_name="j-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_V_2D" />
+     <field id="vtrd_tfri_iceimp"  long_name="j-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_V_2D" />
+     <field id="vtrd_tfr_iceimp"  long_name="j-trend: top friction: as if implicit ice-ocean drag"    unit="m/s^2" grid_ref="grid_V_2D" />
      <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />   
      <field id="vtrd_tot2d"     long_name="j-trend: total momentum trend before atf: 2D component"    unit="m/s^2"  grid_ref="grid_V_2D" />   
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/ICE/iceupdate.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/ICE/iceupdate.F90	(revision 15197)
@@ -365,4 +365,8 @@
       ELSE
          zflagi = 1._wp
+         IF( l_trddyn ) THEN
+            ! Save drag with right sign to update alternative top drag trend diagnostic
+            rCdU_ice(:,:) = -r1_rau0 * tmod_io(:,:) * at_i(:,:) * tmask(:,:,1) 
+         ENDIF
       ENDIF
       !
@@ -384,4 +388,7 @@
                uiceoc(ji,jj) = zat_u * zutau_ice
                viceoc(ji,jj) = zat_v * zvtau_ice
+               ! save part of the ice-ocean stress as would do for the implicit ice-ocean drag calculation
+               uiceoc_iceimp(ji,jj) = zat_u * 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * u_ice(ji,jj)
+               viceoc_iceimp(ji,jj) = zat_v * 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * v_ice(ji,jj)
             ENDIF
             !
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/DYN/dynzdf.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/DYN/dynzdf.F90	(revision 15197)
@@ -16,4 +16,5 @@
    USE dom_oce        ! ocean space and time domain variables 
    USE sbc_oce        ! surface boundary condition: ocean
+   USE sbc_ice, ONLY: rCdU_ice  ! ice-ocean drag coefficient
    USE zdf_oce        ! ocean vertical physics variables
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
@@ -174,4 +175,16 @@
             ENDIF
          END IF
+         IF( l_trddyn .AND. (.NOT.ln_drgice_imp) ) THEN
+            ! Output the contribution to the top friction trend as if using implicit ice-ocean drag
+            DO jj = 2, jpjm1        
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
+                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
+                  ztrdu_fr(ji,jj,1) = 0.5*( rCdU_ice(ji+1,jj)+rCdU_ice(ji,jj) ) * ua_b(ji,jj) / ze3ua
+                  ztrdv_fr(ji,jj,1) = 0.5*( rCdU_ice(ji,jj+1)+rCdU_ice(ji,jj) ) * va_b(ji,jj) / ze3va
+               END DO
+            END DO
+            CALL trd_dyn( ztrdu_fr(:,:,1), ztrdv_fr(:,:,1), jpdyn_tfre_bt_iceimp, kt )
+         END IF
      ENDIF
       !
@@ -521,4 +534,12 @@
             ztrdv_fr(:,:,1) = ( viceoc_b(:,:) + viceoc(:,:) ) / ( e3v_n(:,:,1) * rau0 )
             CALL trd_dyn( ztrdu_fr(:,:,1), ztrdv_fr(:,:,1), jpdyn_iceoc, kt )
+            IF( .NOT. ln_drgice_imp ) THEN
+               !                           ! trends due to ice-ocean drag in surface layer 
+               !                           ! as if using implicit ice ocean drag
+               ztrdu_fr(:,:,:) = 0._wp    ;   ztrdv_fr(:,:,:) = 0._wp
+               ztrdu_fr(:,:,1) = ( uiceoc_iceimp_b(:,:) + uiceoc_iceimp(:,:) ) / ( e3u_n(:,:,1) * rau0 )
+               ztrdv_fr(:,:,1) = ( viceoc_iceimp_b(:,:) + viceoc_iceimp(:,:) ) / ( e3v_n(:,:,1) * rau0 )
+               CALL trd_dyn( ztrdu_fr(:,:,1), ztrdv_fr(:,:,1), jpdyn_iceoc_iceimp, kt )
+            ENDIF
          ENDIF
          !
@@ -558,4 +579,16 @@
                CALL trd_dyn( ztrdu_fr, ztrdv_fr, jpdyn_tfri, kt )
             ENDIF
+            IF( .NOT. ln_drgice_imp ) THEN  ! trends as if using implicit ice-ocean drag
+               ztrdu_fr(:,:,:) = 0._wp    ;   ztrdv_fr(:,:,:) = 0._wp
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     ztrdu_fr(ji,jj,1) = 0.5 * ( rCdU_ice(ji+1,jj) + rCdU_ice(ji,jj) )  & 
+     &                                      * un(ji,jj,1) / e3u_n(ji,jj,1)
+                     ztrdv_fr(ji,jj,1) = 0.5 * ( rCdU_ice(ji,jj+1) + rCdU_ice(ji,jj) )  &
+     &                                      * vn(ji,jj,1) / e3v_n(ji,jj,1)
+                  END DO
+               END DO
+               CALL trd_dyn( ztrdu_fr(:,:,1), ztrdv_fr(:,:,1), jpdyn_tfri_iceimp, kt )
+            ENDIF
          ENDIF
          !
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbc_oce.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbc_oce.F90	(revision 15197)
@@ -113,4 +113,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   uiceoc, uiceoc_b !: sea surface i-stress (ice-ocean)             [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   viceoc, viceoc_b !: sea surface j-stress (ice-ocean)             [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   uiceoc_iceimp, uiceoc_iceimp_b !: sea surface i-stress (ice-ocean)             [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   viceoc_iceimp, viceoc_iceimp_b !: sea surface j-stress (ice-ocean)             [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum             !: module of sea surface stress (at T-point)    [N/m2] 
    !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbcmod.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/SBC/sbcmod.F90	(revision 15197)
@@ -422,4 +422,6 @@
             uiceoc_b(:,:) = uiceoc(:,:)      
             viceoc_b(:,:) = viceoc(:,:)      
+            uiceoc_iceimp_b(:,:) = uiceoc_iceimp(:,:)      
+            viceoc_iceimp_b(:,:) = viceoc_iceimp(:,:)      
          ENDIF        
          !
@@ -541,4 +543,6 @@
             uiceoc_b(:,:) = uiceoc(:,:)      
             viceoc_b(:,:) = viceoc(:,:)      
+            uiceoc_iceimp_b(:,:) = uiceoc_iceimp(:,:)      
+            viceoc_iceimp_b(:,:) = viceoc_iceimp(:,:)      
          ENDIF        
          !
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trd_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trd_oce.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trd_oce.F90	(revision 15197)
@@ -90,4 +90,9 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfre_bt  = 25  !: bottom friction due to barotropic currents for baroclinic trend (ln_dynspg_ts=.TRUE.)
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri     = 26  !: implicit bottom friction for baroclinic trend (ln_drgimp=.TRUE.)
+   !                               ================     !: EXTRA DIAGNOSTIC FLAGS - AS IF USING IMPLICIT ICE-OCEAN DRAG
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_iceoc_iceimp = 27  !: ice-ocean drag trend as if using implicit ice-ocean drag
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tfri_iceimp  = 28  !: implicit bottom friction for baroclinic trend (ln_drgimp=.FALSE.)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tfre_bt_iceimp=29  !: top friction due to barotropic currents for baroclinic trend (ln_dynspg_ts=.FALSE.)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tfr_iceimp   = 30  !: top friction trend as if using implicit ice-ocean drag
    !
    !!----------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trddyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trddyn.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trddyn.F90	(revision 15197)
@@ -49,4 +49,6 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_tfr2d, zvtrd_tfr2d
    REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_bfr2d, zvtrd_bfr2d
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_tfr_iceimp, zvtrd_tfr_iceimp
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_iceoc_iceimp, zvtrd_iceoc_iceimp
 
    !! * Substitutions
@@ -237,4 +239,18 @@
             DEALLOCATE( zue, zve, zutrd_tfr, zvtrd_tfr )
          ENDIF
+         IF( ALLOCATED( zutrd_tfr_iceimp ) ) THEN
+            ! Finalise top friction trend as if using implicit ice-ocean drag
+            IF( ALLOCATED( zutrd_iceoc_iceimp ) ) THEN
+               ! Add trend due to ice-ocean stress at the surface
+               zutrd_tfr_iceimp(:,:) = zutrd_tfr_iceimp(:,:) + zutrd_iceoc_iceimp(:,:)
+               zvtrd_tfr_iceimp(:,:) = zvtrd_tfr_iceimp(:,:) + zvtrd_iceoc_iceimp(:,:)
+               DEALLOCATE( zutrd_iceoc_iceimp, zvtrd_iceoc_iceimp )
+            ENDIF
+            ! Simulate "removal of depth-mean component" as if it was a 3D field:
+            zutrd_tfr_iceimp(:,:) = ( zutrd_tfr_iceimp(:,:) * ( 1._wp - r1_hu_a(:,:) ) ) * umask(:,:,1)
+            zvtrd_tfr_iceimp(:,:) = ( zvtrd_tfr_iceimp(:,:) * ( 1._wp - r1_hv_a(:,:) ) ) * vmask(:,:,1)
+            CALL trd_dyn_iom_2d( zutrd_tfr_iceimp, zvtrd_tfr_iceimp, jpdyn_tfr_iceimp, kt )
+            DEALLOCATE( zutrd_tfr_iceimp, zvtrd_tfr_iceimp )
+         ENDIF
          IF( ALLOCATED( zutrd_bfr ) ) THEN
             ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) )
@@ -353,4 +369,27 @@
           zutrd_iceoc(:,:) = putrd(:,:)
           zvtrd_iceoc(:,:) = pvtrd(:,:)
+
+      CASE( jpdyn_iceoc_iceimp )
+          !
+          ! Save surface ice-ocean stress trend as if using implicit ice-ocean drag.
+          ! Will be add to TFR_ICEIMP trend later.
+          IF( .NOT. ALLOCATED(zutrd_iceoc_iceimp) ) ALLOCATE( zutrd_iceoc_iceimp(jpi,jpj), zvtrd_iceoc_iceimp(jpi,jpj) )
+          zutrd_iceoc_iceimp(:,:) = putrd(:,:)
+          zvtrd_iceoc_iceimp(:,:) = pvtrd(:,:)
+          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_iceoc_iceimp, kt )
+
+      CASE( jpdyn_tfre_bt_iceimp, jpdyn_tfri_iceimp )
+          !
+          ! Add various top friction (as if using implicit ice ocean drag)
+          ! terms for baroclinic trend to saved quantity.
+          ! Any depth-mean component removed later when TFR_ICEIMP trend written out. 
+          IF( .NOT. ALLOCATED(zutrd_tfr_iceimp) ) THEN 
+             ALLOCATE( zutrd_tfr_iceimp(jpi,jpj), zvtrd_tfr_iceimp(jpi,jpj) )
+             zutrd_tfr_iceimp(:,:) = 0.0
+             zvtrd_tfr_iceimp(:,:) = 0.0
+          ENDIF
+          zutrd_tfr_iceimp(:,:) = zutrd_tfr_iceimp(:,:) + putrd(:,:) 
+          zvtrd_tfr_iceimp(:,:) = zvtrd_tfr_iceimp(:,:) + pvtrd(:,:)
+          CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
 
       CASE( jpdyn_tau )
@@ -530,4 +569,12 @@
       CASE( jpdyn_tfr )      ;   CALL iom_put( "utrd_tfr2d", putrd )      ! top friction depth-mean trend
                                  CALL iom_put( "vtrd_tfr2d", pvtrd )
+      CASE( jpdyn_iceoc_iceimp ) ;   CALL iom_put( "utrd_iceoc_iceimp", putrd )  ! total top friction for top layer                    
+                                     CALL iom_put( "vtrd_iceoc_iceimp", pvtrd )  ! as if with implicit ice-ocean drag
+      CASE( jpdyn_tfre_bt_iceimp );  CALL iom_put( "utrd_tfre_bt_iceimp", putrd )  ! total top friction for top layer                    
+                                     CALL iom_put( "vtrd_tfre_bt_iceimp", pvtrd )  ! as if with implicit ice-ocean drag
+      CASE( jpdyn_tfri_iceimp )   ;  CALL iom_put( "utrd_tfri_iceimp", putrd )  ! total top friction for top layer                    
+                                     CALL iom_put( "vtrd_tfri_iceimp", pvtrd )  ! as if with implicit ice-ocean drag
+      CASE( jpdyn_tfr_iceimp )   ;   CALL iom_put( "utrd_tfr_iceimp", putrd )  ! total top friction for top layer                    
+                                     CALL iom_put( "vtrd_tfr_iceimp", pvtrd )  ! as if with implicit ice-ocean drag
       CASE( jpdyn_tot )      ;   CALL iom_put( "utrd_tot2d", putrd )      ! total 2D trend, excluding time filter
                                  CALL iom_put( "vtrd_tot2d", pvtrd )
Index: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trdini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trdini.F90	(revision 15196)
+++ NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends_iceoc_drag/src/OCE/TRD/trdini.F90	(revision 15197)
@@ -85,5 +85,7 @@
       ! Allocate (partial) ice-ocean stresses (only used for dynamics trends diagnostics). 
       IF( l_trddyn .and. nn_ice == 2 ) ALLOCATE( uiceoc(jpi,jpj), uiceoc_b(jpi,jpj), &
-                                                 viceoc(jpi,jpj), viceoc_b(jpi,jpj) )
+                                                 viceoc(jpi,jpj), viceoc_b(jpi,jpj), &
+                                                 uiceoc_iceimp(jpi,jpj), uiceoc_iceimp_b(jpi,jpj), &
+                                                 viceoc_iceimp(jpi,jpj), viceoc_iceimp_b(jpi,jpj) )
 
 !!gm check the stop below      
