diff --git a/src/ICE/ice.F90 b/src/ICE/ice.F90
index 5b0b82ecac7b1221f976e7cdad4d6fbeb722bad9..00b824385f83b2141eec714d42b63a73946afe40 100644
--- a/src/ICE/ice.F90
+++ b/src/ICE/ice.F90
@@ -453,6 +453,11 @@ MODULE ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top     !: Surface conduction flux (W/m2)
    !
    !!----------------------------------------------------------------------
+   !! * Only for atmospheric coupling
+   !!----------------------------------------------------------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
+   !
+   !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
    !! Software governed by the CeCILL license (see ./LICENSE)
@@ -551,6 +556,10 @@ CONTAINS
       ii = ii + 1
       ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
 
+      ! * For atmospheric coupling
+      ii = ii + 1
+      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) )
+
       ice_alloc = MAXVAL( ierr(:) )
       IF( ice_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
       !
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index a15f18d01c5af53b905a14daef8c631be8cf0aba..d2dbc9698cb3f29d76b55313130795c186105856 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -996,24 +996,22 @@ CONTAINS
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE ice_dyn_1d2d ***
       !!
-      !! ** Purpose :   move arrays from 1d to 2d and the reverse
+      !! ** Purpose :   move arrays between 1d <=> 2d forms and
+      !!                            between 2d <=> 3d forms
       !!-----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
+      INTEGER, INTENT(in) ::   kn   ! 1= from 2D/3D to 1D/2D   
+      !                               2= from 1D/2D to 2D/3D
       !
       INTEGER ::   jl, jk   ! dummy loop indices
       !!-----------------------------------------------------------------------
       !
       SELECT CASE( kn )
-      !                    !---------------------!
-      CASE( 1 )            !==  from 2D to 1D  ==!
-         !                 !---------------------!
+      !                    !---------------------------!
+      CASE( 1 )            !==  from 2D/3D to 1D/2D  ==!
+         !                 !---------------------------!
          ! fields used but not modified
          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
-         ! the following fields are modified in this routine
-         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
-         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
-         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
@@ -1035,9 +1033,9 @@ CONTAINS
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
          !
-         !                 !---------------------!
-      CASE( 2 )            !==  from 1D to 2D  ==!
-         !                 !---------------------!
+         !                 !---------------------------!
+      CASE( 2 )            !==  from 1D/2D to 2D/3D  ==!
+         !                 !---------------------------!
          CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
diff --git a/src/ICE/icerst.F90 b/src/ICE/icerst.F90
index ebb65f886360ec07175031de86ee3ffadc2553a9..f7917a5c9b81b711bd4a96453f9c95a50565bc76 100644
--- a/src/ICE/icerst.F90
+++ b/src/ICE/icerst.F90
@@ -324,6 +324,11 @@ CONTAINS
             ENDIF
          ENDIF
 
+         ! If this is a coupled model we need to pick up a_i for use as a_i_last_couple
+         IF (ln_cpl) then
+            a_i_last_couple = a_i
+         ENDIF
+         
          IF(.NOT.lrxios) CALL iom_delay_rst( 'READ', 'ICE', numrir )   ! read only ice delayed global communication variables
          !                 ! ---------------------------------- !
       ELSE                 ! == case of a simplified restart == !
diff --git a/src/ICE/icetab.F90 b/src/ICE/icetab.F90
index 2b1b880f55bc92ee24f8c057ae6655fd4267ef43..f16393877dd801afc41d980482d7fe78ab544b8e 100644
--- a/src/ICE/icetab.F90
+++ b/src/ICE/icetab.F90
@@ -32,14 +32,15 @@ MODULE icetab
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d )
+   SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab2d, tab3d )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tab_2d_1d  ***
+      !!                  ***  ROUTINE tab_3d_2d  ***
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   ndim1d   ! 1d size
       INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field
-      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab1d    ! output 1D field
+
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab3d    ! input 3D field
+      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab2d    ! output 2D field
       !
       INTEGER ::   jl, jn, jid, jjd
       !!----------------------------------------------------------------------
@@ -47,7 +48,7 @@ CONTAINS
          DO jn = 1, ndim1d
             jid          = MOD( tab_ind(jn) - 1 , jpi ) + 1
             jjd          =    ( tab_ind(jn) - 1 ) / jpi + 1
-            tab1d(jn,jl) = tab2d(jid,jjd,jl)
+            tab2d(jn,jl) = tab3d(jid,jjd,jl)
          END DO
       END DO
    END SUBROUTINE tab_3d_2d
@@ -72,14 +73,14 @@ CONTAINS
    END SUBROUTINE tab_2d_1d
 
 
-   SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d )
+   SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab2d, tab3d )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tab_2d_1d  ***
+      !!                  ***  ROUTINE tab_2d_3d  ***
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   ndim1d    ! 1D size
       INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index
-      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab2d     ! output 2D field
+      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab2d     ! input 2D field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab3d     ! output 3D field
       !
       INTEGER ::   jl, jn, jid, jjd
       !!----------------------------------------------------------------------
@@ -87,7 +88,7 @@ CONTAINS
          DO jn = 1, ndim1d
             jid               = MOD( tab_ind(jn) - 1 ,  jpi ) + 1
             jjd               =    ( tab_ind(jn) - 1 ) / jpi  + 1
-            tab2d(jid,jjd,jl) = tab1d(jn,jl)
+            tab3d(jid,jjd,jl) = tab2d(jn,jl)
          END DO
       END DO
    END SUBROUTINE tab_2d_3d
diff --git a/src/OCE/DOM/domvvl.F90 b/src/OCE/DOM/domvvl.F90
index 70fc4651dfafadd7037c8baacbdf5496c237e5a8..8b6c643bcddd66bc7cb38b0ad161f7ef1a612872 100644
--- a/src/OCE/DOM/domvvl.F90
+++ b/src/OCE/DOM/domvvl.F90
@@ -692,7 +692,7 @@ CONTAINS
       !!                - vertical interpolation: simple averaging
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)   ::  pe3_out   ! output interpolated e3
       CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
       !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
       !
diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index 44cc31e74db222f362470b6ad075c5fc90285632..ddb41232e94fc983254c8dce2c41cf156b91a9b0 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -167,7 +167,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       REAL(wp), DIMENSION(jpi,jpj)  :: zhV! fluxes
-      REAL(dp), DIMENSION(jpi,jpj)  :: zhU! fluxes
+      REAL(wp), DIMENSION(jpi,jpj)  :: zhU! fluxes
 !!st#if defined key_qco 
 !!st      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
 !!st#endif
@@ -270,6 +270,11 @@ CONTAINS
             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
          END_2D
+      ELSE
+         ! Ensure zhU and zhV are initialised to SOMETHING at all points to avoid referencing
+         ! uninitialsed values in halos later on!
+         zhU(:,:) = 0._wp
+         zhV(:,:) = 0._wp
       ENDIF
       !
       !                                   !=  Add bottom stress contribution from baroclinic velocities  =!
@@ -502,13 +507,13 @@ CONTAINS
          IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
          !      
          !                             ! resulting flux at mid-step (not over the full domain)
+
          DO_2D( 1, 0, 1, 1 )   ! not jpi-column
             zhU(ji,jj) = e2u(ji,jj) * ua_e(ji,jj) * zhup2_e(ji,jj)
          END_2D
          DO_2D( 1, 1, 1, 0 )   ! not jpj-row
             zhV(ji,jj) = e1v(ji,jj) * va_e(ji,jj) * zhvp2_e(ji,jj)
          END_2D
-         !
 #if defined key_agrif
          ! Set fluxes during predictor step to ensure volume conservation
          IF( ln_bt_fw )   CALL agrif_dyn_ts_flux( jn, zhU, zhV )
@@ -519,6 +524,12 @@ CONTAINS
             CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V
             !
          ENDIF    
+
+         ! It seems safest to do this here since zhU and zhV are not initially calculated in halos
+         ! by this code or by wad_Umsk, but halo values (ji-1 and jj-1) ARE required in the zhdiv 
+         ! sea level calculation. The current trunk (Feb 2024) has resolved all these issues by rewriting.
+         CALL lbc_lnk( 'dynspg_ts',  zhU, 'U', -1._wp)
+         CALL lbc_lnk( 'dynspg_ts',  zhV, 'V', -1._wp)
          !
          !
          !     Compute Sea Level at step jit+1
@@ -529,15 +540,16 @@ CONTAINS
             zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
             ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( ssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
          END_2D
-         !
-         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp,  zhU, 'U', -1._dp)
-         CALL lbc_lnk( 'dynspg_ts',  zhV, 'V', -1._wp )
+
+         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp)
+
          !
          ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
          IF( ln_bdy )   CALL bdy_ssh( ssha_e )
 #if defined key_agrif
          CALL agrif_ssh_ts( jn )
 #endif
+
          !
          !                             ! Sum over sub-time-steps to compute advective velocities
          za2 = wgtbtp2(jn)             ! zhU, zhV hold fluxes extrapolated at jn+0.5
@@ -556,6 +568,7 @@ CONTAINS
          !  
          ! Sea Surface Height at u-,v-points (vvl case only)
          IF( .NOT.ln_linssh ) THEN
+
 #if defined key_qcoTest_FluxForm
             !                                ! 'key_qcoTest_FluxForm' : simple ssh average
             DO_2D( 1, 0, 1, 1 )
@@ -623,6 +636,7 @@ CONTAINS
          !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
          !--         h     \                                                                                                 /  --!
          !------------------------------------------------------------------------------------------------------------------------!
+
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
             DO_2D( 0, 0, 0, 0 )
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
@@ -676,7 +690,7 @@ CONTAINS
          ENDIF
        
          IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
-            DO_2D( 0, 0, 0, 0 )
+           DO_2D( 0, 0, 0, 0 )
                hu_e (ji,jj) =    hu_0(ji,jj) + zsshu_a(ji,jj)
                hur_e(ji,jj) = ssumask(ji,jj) / (  hu_e(ji,jj) + 1._wp - ssumask(ji,jj)  )
                hv_e (ji,jj) =    hv_0(ji,jj) + zsshv_a(ji,jj)
@@ -686,10 +700,10 @@ CONTAINS
          !
          IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
             CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
-                 &                   , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
-                 &                   , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
+                                     , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
+                                     , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
          ELSE
-            CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
+            CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
          ENDIF
          !                                                 ! open boundaries
          IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
@@ -728,7 +742,6 @@ CONTAINS
          !                                          ! Sum sea level
          pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
 
-         !                                                 ! ==================== !
       END DO                                               !        end loop      !
       !                                                    ! ==================== !
       ! -----------------------------------------------------------------------------
@@ -1085,7 +1098,13 @@ CONTAINS
       !                             ! Allocate time-splitting arrays
       IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
       !
-      !                             ! read restart when needed
+
+      ! RSRH. I've just copied the following from the current NEMO main. It shouldnt affect evolution (!) but 
+      ! does help with debugging and testing!
+      ! init some arrays for debug sette
+      ssha_e(:,:) = 0._wp
+      !
+      !                      !: restart/initialise
       CALL ts_rst( nit000, 'READ' )
       !
    END SUBROUTINE dyn_spg_ts_init
@@ -1165,8 +1184,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       INTEGER  ::   ji ,jj                             ! dummy loop indices
       REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  :: pht, phu, phv, punb, pvnb, zhV
-      REAL(dp), DIMENSION(jpi,jpj), INTENT(in   )  :: zhU
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  :: pht, phu, phv, punb, pvnb, zhU, zhV
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
       !!----------------------------------------------------------------------
       SELECT CASE( nvor_scheme )
@@ -1272,8 +1290,7 @@ CONTAINS
       !! ** Action  :  ptmsk : wetting & drying t-mask
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  :: phV, pu, pv! ocean velocities and transports
-      REAL(dp), DIMENSION(jpi,jpj), INTENT(inout)  :: phU! ocean velocities and transports
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
       !
       INTEGER  ::   ji, jj   ! dummy loop indices
diff --git a/src/OCE/LBC/mppini.F90 b/src/OCE/LBC/mppini.F90
index 1ba6029bad619a4cafcc79f50dea386a3ae62733..6bc449de86f9a16f9cb6c57a2520bcc54a390972 100644
--- a/src/OCE/LBC/mppini.F90
+++ b/src/OCE/LBC/mppini.F90
@@ -1436,7 +1436,7 @@ ENDIF
      ! This is rather confused by the fact that jpjglo has a value BIGGER
      ! than it did at pre 4.2... e.g. for ORCA1 it's set to 333 instead of 332
      ! which is rather baffling and confuses some dimensioning calculations
-     ! if we're not very very careful 
+     ! which previously worked fine pre 4.2. 
 
 
       ! Set up dimensions for old style coupling exchanges on extended grid
diff --git a/src/OCE/SBC/cpl_oasis3.F90 b/src/OCE/SBC/cpl_oasis3.F90
index 56c3954b6718c8dfe527aee9a80ca05423c314e0..2769435c9591b489d89d3da235826938deda68d9 100644
--- a/src/OCE/SBC/cpl_oasis3.F90
+++ b/src/OCE/SBC/cpl_oasis3.F90
@@ -45,7 +45,6 @@ MODULE cpl_oasis3
    PRIVATE
 #if ! defined key_oasis3
    ! Dummy interface to oasis_get if not using oasis
-   ! RSRH Is this really needed
    INTERFACE oasis_get
       MODULE PROCEDURE oasis_get_1d, oasis_get_2d
    END INTERFACE
@@ -253,8 +252,8 @@ CONTAINS
       CALL oasis_def_partition ( id_part_2d, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos
 
       ! RSRH Set up 2D box partition for compatibility with existing weights files
-      ! so we DONT HAVE TO GENERATE AND MANAGE MULTIPLE SETS OF WEIGHTS PURLEY BECAUSE OF AN
-      ! ARBITRARY CHANGE IN THE SOURCE CODE!
+      ! so we don't have to generate and manage multiple sets of weights purely because of 
+      ! the changes to nemo 4.2+ code!
 
       ! This is just a hack for global cyclic models for the time being
       Ni0glo_ext = jpiglo
@@ -569,13 +568,14 @@ CONTAINS
       INTEGER                   , INTENT(  out) ::   kinfo     ! OASIS3 info argument
       !!
       INTEGER                                   ::   jc,jm     ! local loop index
-      LOGICAL                                   ::   llaction, ll_1st
+      LOGICAL                                   ::   llaction, ll_1st, lrcv
 
       !!--------------------------------------------------------------------
       !
       ! receive local data from OASIS3 on every process
       !
       kinfo = OASIS_idle
+      lrcv=.FALSE.
       !
 
       DO jc = 1, srcv(kid)%nct
@@ -594,9 +594,10 @@ CONTAINS
                   &  WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
 
                IF( llaction ) THEN   ! data received from oasis do not include halos
-                                     ! RSRH, but DO cater for wrap columns when using pre vn 4.2 format remapping weights. 
+                                     ! but DO still cater for wrap columns when using pre vn4.2 compatible remapping weights. 
 
                   kinfo = OASIS_Rcv
+                  lrcv=.TRUE. 
                   IF( ll_1st ) THEN
                      pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) =   exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
 
@@ -626,14 +627,17 @@ CONTAINS
 
          ENDDO
 
-         !--- Call lbc_lnk to populate halos of received fields.
-         IF( .NOT. ll_1st ) THEN
-
-            CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
-
-         ENDIF
-
       ENDDO
+
+      ! RSRH I've changed this since:
+      ! 1) it seems multi cat fields may not be properly updated in halos when called on a per 
+      !    category basis(?) 
+      ! 2) it's more efficient to have a single call (and simpler to understand) to update ALL 
+      !    categories at the same time!
+      !--- Call lbc_lnk to populate halos of received fields.
+      IF (lrcv) then
+         CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,:), srcv(kid)%clgrid, srcv(kid)%nsgn )
+      endif
       !
    END SUBROUTINE cpl_rcv
 
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index 98dcb8b78c206b192ffab1cfa3f46b1026e6648f..b0b669fadf596c196f46fd56e6c6ec331a40441f 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -33,7 +33,7 @@ MODULE sbccpl
 #endif
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      !
-   USE oce     , ONLY : ts, uu, vv, ssh, fraqsr_1lev
+   USE oce   
 #if defined key_medusa
    USE oce , ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl,  &
                         PCO2a_in_cpl, Dust_in_cpl
@@ -254,9 +254,7 @@ MODULE sbccpl
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                ! all fields recieved from the atmosphere
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
-#if defined key_si3 || defined key_cice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
-#endif
+
 
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
 
@@ -1460,6 +1458,14 @@ CONTAINS
                IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
                   ! Temporary code for HadGEM3 - will be removed eventually.
                   ! Only applies when we have only taux on U grid and tauy on V grid
+
+                  !RSRH these MUST be initialised because the halos are not explicitly set 
+                  ! but they're passed to repcmo and used directly in calculations, so if 
+                  ! they point at junk in memory then bad things will happen!
+                  ! (You can prove this by running with preset NaNs). 
+                  ztx(:,:)=0.0
+                  zty(:,:)=0.0
+
                   DO_2D( 0, 0, 0, 0 )                                       
                           ztx(ji,jj)=0.25*vmask(ji,jj,1)                &
                              *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &
@@ -1574,6 +1580,12 @@ CONTAINS
          !
       ENDIF
 
+#if defined key_medusa
+      IF (ln_medusa) THEN
+        IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
+        IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
+      ENDIF
+#endif
       !                                                      ! ================== !
       !                                                      ! atmosph. CO2 (ppm) !
       !                                                      ! ================== !
@@ -1871,6 +1883,7 @@ CONTAINS
          antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
       ENDIF
       !
+      IF (ln_timing) CALL timing_stop('sbc_cpl_rcv')
    END SUBROUTINE sbc_cpl_rcv
 
 
@@ -2070,13 +2083,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !
 #if defined key_si3 || defined key_cice
-      !
-      IF( kt == nit000 ) THEN
-         ! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
-         IF( .NOT.ALLOCATED(a_i_last_couple) )   ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
-         ! initialize to a_i for the 1st time step
-         a_i_last_couple(:,:,:) = a_i(:,:,:)
-      ENDIF
+
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
       ziceld(:,:) = 1._wp - picefr(:,:)
@@ -2276,23 +2283,7 @@ CONTAINS
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
 !!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
       !
-      !                                                      ! ========================= !
-      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
-      !                                                      ! ========================= !
-      CASE ('coupled')
-         IF (ln_scale_ice_flux) THEN
-            WHERE( a_i(:,:,:) > 1.e-10_wp )
-               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-            ELSEWHERE
-               qml_ice(:,:,:) = 0.0_wp
-               qcn_ice(:,:,:) = 0.0_wp
-            END WHERE
-         ELSE
-            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
-            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
-         ENDIF
-      END SELECT
+
       !
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
@@ -2702,8 +2693,12 @@ CONTAINS
       REAL(wp) ::   zumax, zvmax
       REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp5, ztmp6 ! RSRH temporary work arrays 
+                                                         ! to avoid intent conflicts in repcmo calls
       !!----------------------------------------------------------------------
       !
+      IF (ln_timing) CALL timing_start('sbc_cpl_snd')
+
       isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
       info = OASIS_idle
 
@@ -3069,13 +3064,20 @@ CONTAINS
                           *(zoty1(ji,jj)+zoty1(ji+1,jj)    &
                           +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
                END_2D
-              
-               ! Ensure any N fold and wrap columns are updated
-               CALL lbc_lnk('zotx1', ztmp1, 'V', -1.0)
-               CALL lbc_lnk('zoty1', ztmp2, 'U', -1.0)
-	               
+              	               
                ikchoix = -1
-               CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)
+               ! zotx1 and zoty1 are input only to repcmo while ztmp5 and ztmp6
+               ! are the newly calculated (output) values.
+               ! Don't make the mistake of using zotx1 and zoty1 twice in this
+               ! call for both input and output fields since it creates INTENT
+               ! conflicts. 
+               CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,ztmp5,ztmp6,ikchoix)
+               zotx1(:,:)=ztmp5(:,:)
+               zoty1(:,:)=ztmp6(:,:)
+
+               ! Ensure any N fold and wrap columns are updated. 
+               CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
+
             ENDIF
          ENDIF
          !
diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90
index 878841a8eb28f6f811b3d685d18b756228bd1053..ec1ba7773fe44ab3e7e968b18018ec579fe14ca7 100644
--- a/src/OCE/TRA/traadv_fct.F90
+++ b/src/OCE/TRA/traadv_fct.F90
@@ -87,6 +87,7 @@ CONTAINS
       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
       REAL(wp), DIMENSION(A2D(nn_hls),jpk)         :: zwi, zwz, ztu, ztv, zltu, zltv, ztw
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)         :: zwi_in, ztw_in !temp arrays to avoid intent conflicts
       REAL(dp), DIMENSION(A2D(nn_hls),jpk)         :: zwx, zwy
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
@@ -194,7 +195,9 @@ CONTAINS
          END_3D
 
          IF ( ll_zAimp ) THEN
-            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+            ! We need to use separate copies of zwi to avoid intent conflicts!
+            zwi_in(:,:,:) = zwi(:,:,:)
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
@@ -314,7 +317,10 @@ CONTAINS
                ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
             END_3D
             !
-            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
+
+            ! We need to use separate copies of ztw to avoid intent conflicts!
+            ztw_in(:,:,:) = ztw(:,:,:)
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
             !
             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
@@ -724,6 +730,8 @@ CONTAINS
       REAL(wp) ::   zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v   !   -      -
       REAL(wp) ::   ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi_in, ztw_in ! Read-only copies to avoid INTENT conflicts
+                                                                  ! in calls to tridia_solver
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
       LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
@@ -739,6 +747,8 @@ CONTAINS
       zwy_3d(:,:,:) = 0._wp
       zwz(:,:,:) = 0._wp
       zwi(:,:,:) = 0._wp
+
+      zwt(:,:,:) = 0._wp
       !
       l_trd = .FALSE.            ! set local switches
       l_hst = .FALSE.
@@ -820,7 +830,10 @@ CONTAINS
          END DO
 
          IF ( ll_zAimp ) THEN
-            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+
+            ! We need to use separate copies of zwi to avoid intent conflicts!
+            zwi_in(:,:,:) = zwi(:,:,:)
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
             DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
@@ -934,7 +947,9 @@ CONTAINS
                ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
             END_3D
             !
-            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
+            ! We need to use separate copies of ztw to avoid intent conflicts!
+            ztw_in(:,:,:) = ztw(:,:,:)
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
             !
             DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
diff --git a/src/OCE/nemogcm.F90 b/src/OCE/nemogcm.F90
index e8461cf9a1511e338c8af507cbbd0e54b459c994..80d3653c37cd0ce5bd84856cfc5bea8160964c57 100644
--- a/src/OCE/nemogcm.F90
+++ b/src/OCE/nemogcm.F90
@@ -174,12 +174,12 @@ CONTAINS
                IF ( istp ==         nitend ) elapsed_time = zstptiming - elapsed_time
             ENDIF
             !
+            IF (lk_oasis) THEN
+               CALL sbc_cpl_snd( istp, Nbb, Nnn )  ! Coupling to atmos
+            ENDIF
 #  if defined key_qco   ||   defined key_linssh
             CALL stp_MLF( istp )
 #  else
-	    IF (lk_oasis) THEN
-               CALL sbc_cpl_snd( istp, Nbb, Nnn )  ! Coupling to atmos
-            ENDIF
             CALL stp    ( istp )
 #  endif
             istp = istp + 1
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 17e44e46e363fcce36bcb1fc77c5c85fdb6efa51..176f04b9c36fd91ef694af10dbfa0b1f755a9899 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -443,7 +443,7 @@ CONTAINS
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
+      !!IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
       !
 #if defined key_xios
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -561,6 +561,7 @@ CONTAINS
          CALL lbc_lnk( 'finalize_lbc', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, &
             &                          r3u_f(:,:),   'U', 1._wp, r3v_f(:,:),   'V', 1._wp )
       ENDIF
+
       !                                        !* BDY open boundaries
       IF( ln_bdy )   THEN
                                CALL bdy_tra( kt, Kbb, pts,      Kaa )