Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8892)
@@ -480,5 +480,8 @@
       ii = ii + 1
       ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
-      
+at_i_b(:,:) = 0.0      
+u_ice_b(:,:) = 0.0      
+v_ice_b(:,:) = 0.0      
+
       ! * Ice thickness distribution variables
       ii = ii + 1
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8892)
@@ -103,5 +103,5 @@
       !                             !-----------------------------------------------------
       IF ( jpl > 1 )   CALL ice_itd_reb( kt )
-
+write(numout,*) "RSRH out of ice_itd_reb in icecor", kt ; flush(numout)
       !                             !-----------------------------------------------------
       CALL ice_var_zapsmall         !  Zap small values                                  !
@@ -133,4 +133,5 @@
       CASE( 1 )                        !--- dyn trend diagnostics
          !
+write(numout,*) "RSRH GG icecor CASE 1", kt ; flush(numout)
 !!gm   here I think the number of ice cat is too small to use a SUM instruction...
          DO jj = 1, jpj
@@ -152,4 +153,5 @@
       CASE( 2 )                        !--- thermo trend diagnostics & ice aging
          !
+write(numout,*) "RSRH GG icecor CASE 2", kt ; flush(numout)
          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
          !
@@ -173,4 +175,5 @@
          !
       END SELECT
+write(numout,*) "RSRH HH icecor", kt ; flush(numout)
       !
       ! controls
@@ -179,4 +182,5 @@
       IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints
       IF( nn_timing == 1 )   CALL timing_stop ('icecor')                                                             ! timing
+write(numout,*) "RSRH leaving icecor", kt ; flush(numout)
       !
    END SUBROUTINE ice_cor
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8892)
@@ -146,5 +146,10 @@
          !
       CASE ( jp_purecpl )           !--- coupled formulation
+write(numout,*) "RSRH call cpl_ice_flx at_i_b",kt,glob_sum(at_i_b(:,:));flush(numout)
+write(numout,*) "RSRH call cpl_ice_flx alb_ice",kt,glob_sum(alb_ice(:,:));flush(numout)
+write(numout,*) "RSRH call cpl_ice_flx sst_m",kt,glob_sum(sst_m(:,:));flush(numout)
+
                                 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
+write(numout,*) "RSRH call cpl_ice_flx",kt,glob_sum(at_i_b(:,:));flush(numout)
          IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8892)
@@ -365,7 +365,17 @@
          DO jl = 1, jpl
             a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac
+WRITE(numout,*) "RSRH iceistate AA: a_ip_frac", glob_sum(a_ip_frac(:,:,:)) ; flush(numout)
+WRITE(numout,*) "RSRH iceistate AA: zswitch", glob_sum(zswitch(:,:)) ; flush(numout)
+WRITE(numout,*) "RSRH iceistate BB: h_ip pre calc", glob_sum(h_ip(:,:,:)) ; flush(numout)
+
             h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac
+WRITE(numout,*) "RSRH iceistate BB: h_ip post calc", glob_sum(h_ip(:,:,:)) ; flush(numout)
          END DO
+WRITE(numout,*) "RSRH iceistate CC: a_ip_frac pre a_ip calc", glob_sum(a_ip_frac(:,:,:)) ; flush(numout)
+WRITE(numout,*) "RSRH iceistate CC: a_i pre a_ip calc", glob_sum(a_i(:,:,:)) ; flush(numout)
+WRITE(numout,*) "RSRH iceistate CC: a_ip pre a_ip calc", glob_sum(a_ip(:,:,:)) ; flush(numout)
          a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
+WRITE(numout,*) "RSRH iceistate DD: h_ip pre v_ip calc", glob_sum(h_ip(:,:,:)) ; flush(numout)
+WRITE(numout,*) "RSRH iceistate DD: a_ip pre v_ip calc", glob_sum(a_ip(:,:,:)) ; flush(numout)
          v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
 
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8892)
@@ -26,5 +26,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lib_fortran    ! fortran utilities (sum + no signed zero)
    USE prtctl         ! Print control
 
@@ -108,5 +108,6 @@
          zdhice(:,:) = 0._wp
          zhbnew(:,:) = 0._wp
-
+write(numout,*) "RSRH PE iceitd AA h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd AA a_i", sum(a_i(:,:,:)); flush(numout)
          CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
          CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
@@ -197,4 +198,6 @@
          DO jl = 1, jpl
             !
+write(numout,*) "RSRH PE iceitd BB h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd BB a_i", sum(a_i(:,:,:)); flush(numout)
             CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
             CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
@@ -241,4 +244,6 @@
                END DO
                !
+write(numout,*) "RSRH PE iceitd CC h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd CC a_i", sum(a_i(:,:,:)); flush(numout)
                CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
                CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
@@ -294,7 +299,13 @@
          ! 7) Make sure h_i >= minimum ice thickness hi_min
          !----------------------------------------------------------------------------------------------
+write(numout,*) "RSRH PE iceitd DD h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd DD a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd DD a_ip", sum(a_ip(:,:,:)); flush(numout)
          CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
+write(numout,*) "RSRH PE iceitd EE h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd EE a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd EE a_ip", sum(a_ip(:,:,:)); flush(numout)
          
          DO ji = 1, npti
@@ -306,7 +317,13 @@
          END DO
          !
+write(numout,*) "RSRH PE iceitd FF h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd FF a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd FF a_ip", sum(a_ip(:,:,:)); flush(numout)
          CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
+write(numout,*) "RSRH PE iceitd GG h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd GG a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd GG a_ip", sum(a_ip(:,:,:)); flush(numout)
          !
       ENDIF
@@ -527,4 +544,7 @@
       END WHERE
       !
+write(numout,*) "RSRH PE iceitd HH h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd HH a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd HH a_ip", sum(a_ip(:,:,:)); flush(numout)
       CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
       CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
@@ -536,4 +556,7 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
+write(numout,*) "RSRH PE iceitd II h_i", sum(h_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd II a_i", sum(a_i(:,:,:)); flush(numout)
+write(numout,*) "RSRH PE iceitd II a_ip", sum(a_ip(:,:,:)); flush(numout)
       !
    END SUBROUTINE ice_itd_shiftice
@@ -576,7 +599,14 @@
          END DO
          !
+write(numout,*) "RSRH PE iceitd KK h_i",jl, sum(h_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd KK a_i",jl, sum(a_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd KK a_ip",jl, sum(a_ip(:,:,jl)); flush(numout)
+        !
 !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
+write(numout,*) "RSRH PE iceitd LL h_i",jl, sum(h_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd LL a_i",jl, sum(a_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd LL a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
          !
          DO ji = 1, npti
@@ -606,4 +636,5 @@
       !                       !-----------------------------------------
       DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
+write(numout,*) "RSRH PE iceitd LL for jl=",jl ;flush(numout)
          !                    !-----------------------------------------
          npti = 0 ; nptidx(:) = 0
@@ -616,7 +647,15 @@
             END DO
          END DO
-         !
+write(numout,*) "RSRH PE iceitd LL after 1st loop for jl=",jl ;flush(numout)
+write(numout,*) "RSRH PE iceitd get summing LL after 1st loop for jl=",jl ;flush(numout)
+         !
+write(numout,*) "RSRH PE iceitd MM pre h_i",jl, sum(h_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd MM pre a_i",jl, sum(a_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd MM pre a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl+1)    ) ! jl+1 is ok
          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl+1)    ) ! jl+1 is ok
+write(numout,*) "RSRH PE iceitd MM h_i",jl, sum(h_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd MM a_i",jl, sum(a_i(:,:,jl)); flush(numout)
+write(numout,*) "RSRH PE iceitd MM a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
          DO ji = 1, npti
             jdonor(ji,jl) = jl + 1
@@ -624,4 +663,5 @@
             zdvice(ji,jl) = v_i_1d(ji)
          END DO
+write(numout,*) "RSRH PE iceitd MM after 2nd loop for jl=",jl ;flush(numout)
          !
          IF( npti > 0 ) THEN
@@ -632,6 +672,8 @@
             zdvice(1:npti,jl) = 0._wp
          ENDIF
+write(numout,*) "RSRH PE iceitd NN after 3rd loop for jl=",jl ;flush(numout)
          !
       END DO
+write(numout,*) "RSRH PE iceitd NN out of NN loop for jl=",jl ;flush(numout)
       !
    END SUBROUTINE ice_itd_reb
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8892)
@@ -107,5 +107,5 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- !
          !                                   !-----------------------!
-         !
+         !meaning
          kt_ice = kt                              ! -- Ice model time step
          !
@@ -116,4 +116,6 @@
          v_oce(:,:) = ssv_m(:,:)
          !
+write(numout,*) "RSRH ice_stp 00 before eos_stp ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp before eos_stp ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
          CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
          t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
@@ -134,6 +136,10 @@
          CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary) 
          !
+write(numout,*) "RSRH ice_stp 00 before ice_dyn ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp before ice_dyn ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
          IF( ln_icedyn .AND. .NOT.lk_c1d )   &
             &                           CALL ice_dyn( kt )            ! -- Ice dynamics
+write(numout,*) "RSRH ice_stp 00 after ice_dyn ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp aft ice_dyn ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
          !
          !                          !==  lateral boundary conditions  ==!
@@ -145,6 +151,11 @@
          !
          !                          !==  previous lead fraction and ice volume for flux calculations
+write(numout,*) "RSRH ice_stp 00 be4 glo2eqv ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp b4 glo2eqv ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
          CALL ice_var_glo2eqv            ! h_i and h_s for ice albedo calculation
+write(numout,*) "RSRH ice_stp 00 be4 ice_var_agg ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          CALL ice_var_agg(1)             ! at_i for coupling 
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp b4 store_filds ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 be4 store_filds ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          CALL store_fields               ! Store now ice values
 
@@ -162,17 +173,29 @@
          !    qprec_ice, qevap_ice
          !------------------------------------------------------!
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp before forc ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 be4 forc ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
                                         CALL ice_forcing_flx( kt, ksbc )
          !----------------------------!
          ! --- ice thermodynamics --- !
          !----------------------------!
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp before thd ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 be4 thd", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
+write(numout,*) "RSRH ice_stp 00 qns,sfx,emp before icecor ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp 00 before icecor ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          !
          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
+write(numout,*) "RSRH ice_stp AA qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp AA ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
 # if defined key_agrif
          IF( .NOT. Agrif_Root() )       CALL agrif_update_lim3( kt )
 # endif
          CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
+write(numout,*) "RSRH ice_stp BB qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp BB ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
+write(numout,*) "RSRH ice_stp CC qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp CC ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
                                      !
 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
@@ -182,14 +205,24 @@
 !!# endif
                                         CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
+write(numout,*) "RSRH ice_stp DD qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp DD ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
 !!# if defined key_agrif
 !!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
 !!# endif
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
+write(numout,*) "RSRH ice_stp EE qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp EE ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          !
                                         CALL ice_wri( kt )            ! -- Ice outputs 
+write(numout,*) "RSRH ice_stp FF qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp FF ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          !
          IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file 
+write(numout,*) "RSRH ice_stp HH qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp HH ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          !
          IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
+write(numout,*) "RSRH ice_stp II qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp II ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
          !
       ENDIF   ! End sea-ice time step only
@@ -199,4 +232,6 @@
       !-------------------------!
       IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses
+write(numout,*) "RSRH ice_stp JJ qns,sfx,emp  ", kt , glob_sum(qns(:,:)), glob_sum(sfx(:,:)),glob_sum(emp(:,:)) ; flush(numout)
+write(numout,*) "RSRH ice_stp JJ ", kt , glob_sum(ssv_m(:,:)), glob_sum(v_ice(:,:)),glob_sum(ssu_m(:,:)), glob_sum(u_ice(:,:)) ; flush(numout)
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
@@ -360,4 +395,8 @@
       
       ! ice velocities & total concentration
+
+
+write(numout,*) "RSRH store_fields a_i_b sum", glob_sum(a_i_b(:,:,:));flush(numout)
+
       at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
       u_ice_b(:,:) = u_ice(:,:)
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8892)
@@ -99,6 +99,13 @@
       ! controls
       IF( nn_timing == 1 )   CALL timing_start('icethd')                                                             ! timing
+
+write(numout,*) "RSRH icethd AA", kt; flush(numout)
       IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-
+write(numout,*) "RSRH icethd BB", kt; flush(numout)
+
+write(numout,*) "RSRH icethd BB u_ice", kt, glob_sum(u_ice(:,:)); flush(numout)
+write(numout,*) "RSRH icethd BB v_ice", kt, glob_sum(v_ice(:,:)); flush(numout)
+write(numout,*) "RSRH icethd BB ssu_m", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH icethd BB ssv_m", kt, glob_sum(ssv_m(:,:)); flush(numout)
       IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
@@ -106,6 +113,11 @@
          WRITE(numout,*) '~~~~~~~'
       ENDIF
+write(numout,*) "RSRH icethd CC", kt; flush(numout)
       
       CALL ice_var_glo2eqv
+write(numout,*) "RSRH icethd DD", kt; flush(numout)
+      
+      ! RSRH. zfric needs to be properly initialised....Duuuuhh!
+      zfric(:,:) = 0.0 
 
       !---------------------------------------------!
@@ -113,6 +125,23 @@
       !---------------------------------------------!
       IF( ln_icedyn ) THEN
+write(numout,*) "RSRH icethd EE in", kt; flush(numout)
          zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
+
+write(numout,*) "RSRH icethd EE u_ice", kt, glob_sum(u_ice(:,:)); flush(numout)
+write(numout,*) "RSRH icethd EE v_ice", kt, glob_sum(v_ice(:,:)); flush(numout)
+write(numout,*) "RSRH icethd EE ssu_m", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH icethd EE ssv_m", kt, glob_sum(ssv_m(:,:)); flush(numout)
+write(numout,*) "RSRH icethd EE zu_io", kt, glob_sum(zu_io(:,:)); flush(numout)
+write(numout,*) "RSRH icethd EE zv_io", kt, glob_sum(zv_io(:,:)); flush(numout)
+
+write(numout,*) "RSRH PE icethd EE u_ice", kt, sum(u_ice(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd EE v_ice", kt, sum(v_ice(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd EE ssu_m", kt, sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd EE ssv_m", kt, sum(ssv_m(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd EE zu_io", kt, sum(zu_io(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd EE zv_io", kt, sum(zv_io(:,:)); flush(numout)
+    
+
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1
@@ -122,5 +151,7 @@
             END DO
          END DO
+write(numout,*) "RSRH icethd EE", kt; flush(numout)
       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
+write(numout,*) "RSRH icethd FF in", kt; flush(numout)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -131,7 +162,25 @@
          END DO
       ENDIF
+write(numout,*) "RSRH icethd GG", kt; flush(numout)
+write(numout,*) "RSRH icethd GG pre zfric link", kt; flush(numout)
+write(numout,*) "RSRH icethd GG sum pre zfric link", kt, glob_sum(zfric(:,:)); flush(numout)
+write(numout,*) "RSRH PE icethd GG sum pre zfric link SUM", kt, sum(zfric(:,:)); flush(numout)
+
       CALL lbc_lnk( zfric, 'T',  1. )
-      !
+ write(numout,*) "RSRH icethd GG after zfric link", kt; flush(numout)
+     !
       ftr_ice(:,:,:) = 0._wp  ! initialization (part of solar radiation transmitted through the ice)
+
+write(numout,*) "RSRH GG at_i", glob_sum(at_i(:,:)); flush(numout)
+write(numout,*) "RSRH GG at_i_b", glob_sum(at_i_b(:,:)); flush(numout)
+write(numout,*) "RSRH GG qsr_oce", glob_sum(qsr_oce(:,:)); flush(numout)
+write(numout,*) "RSRH GG qns_oce", glob_sum(qns_oce(:,:)); flush(numout)
+write(numout,*) "RSRH GG frq_m", glob_sum(frq_m(:,:)); flush(numout)
+write(numout,*) "RSRH GG qemp_oce", glob_sum(qemp_oce(:,:)); flush(numout)
+write(numout,*) "RSRH GG e3t_m", glob_sum(e3t_m(:,:)); flush(numout)
+write(numout,*) "RSRH GG t_bo", glob_sum(t_bo(:,:)); flush(numout)
+write(numout,*) "RSRH GG sst_m", glob_sum(sst_m(:,:)); flush(numout)
+
+write(numout,*) "RSRH GG zfric", glob_sum(zfric(:,:)); flush(numout)
 
       !--------------------------------------------------------------------!
@@ -185,4 +234,5 @@
          END DO
       END DO
+write(numout,*) "RSRH icethd HH", kt; flush(numout)
       
       ! In case we bypass open-water ice formation
@@ -206,4 +256,5 @@
          &           - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt 
                                                                                 !    (fhld should be 0 while bott growth)
+write(numout,*) "RSRH icethd II", kt; flush(numout)
       !-------------------------------------------------------------------------------------------!
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
@@ -211,5 +262,6 @@
       DO jl = 1, jpl
 
-         ! select ice covered grid points
+ write(numout,*) "RSRH icethd JJ for Jl",jl, kt; flush(numout)
+        ! select ice covered grid points
          npti = 0 ; nptidx(:) = 0
          DO jj = 1, jpj
@@ -222,10 +274,13 @@
          END DO
 
+write(numout,*) "RSRH icethd KK for Jl",jl, kt; flush(numout)
          IF( lk_mpp )         CALL mpp_ini_ice( npti , numout )
 
+write(numout,*) "RSRH icethd LL for Jl",jl, kt; flush(numout)
          IF( npti > 0 ) THEN  ! If there is no ice, do nothing.
             !                                                                
                               CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
             !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
+write(numout,*) "RSRH icethd MM for Jl",jl, kt; flush(numout)
             !
             s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here) 
@@ -239,8 +294,11 @@
                               CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
             ENDIF
+write(numout,*) "RSRH icethd NN for Jl",jl, kt; flush(numout)
             !
                               CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
+write(numout,*) "RSRH icethd OO for Jl",jl, kt; flush(numout)
             !
                               CALL ice_thd_temp                     ! --- temperature update --- !
+write(numout,*) "RSRH icethd PP for Jl",jl, kt; flush(numout)
             !
 !!gm please create a new logical (l_thd_lam or a better explicit name) set one for all in icestp.F90 module
@@ -254,8 +312,11 @@
             END IF
             !
+write(numout,*) "RSRH icethd QQ for Jl",jl, kt; flush(numout)
             IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- !
             !
+write(numout,*) "RSRH icethd RR for Jl",jl, kt; flush(numout)
                               CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
             !                                                       ! --- & Move to 2D arrays --- !
+write(numout,*) "RSRH icethd SS for Jl",jl, kt; flush(numout)
             !
             IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
@@ -267,6 +328,8 @@
 
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+write(numout,*) "RSRH icethd UU for Jl",jl, kt; flush(numout)
       !
                            CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
+write(numout,*) "RSRH icethd VV for Jl",jl, kt; flush(numout)
       !                                                    !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
       !                   
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8892)
@@ -44,4 +44,5 @@
    !
    USE in_out_manager ! I/O manager
+   use lib_fortran
    USE iom            ! NetCDF library
    USE lib_mpp        ! distribued memory computing library
@@ -1133,4 +1134,8 @@
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
+
+write(numout,*) "RSRH in sbc_cpl_rcv",kt
+write(numout,*) "RSRH in sbc_cpl_rcv qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
+
       !
       !                                                      ! ======================================================= !
@@ -1222,4 +1227,5 @@
       ENDIF
       !                                                      ! ========================= !
+write(numout,*) "RSRH in sbc_cpl_rcv BB qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !                                                      !    wind stress module     !   (taum)
       !                                                      ! ========================= !
@@ -1284,4 +1290,5 @@
 
       !                                                      ! ================== !
+write(numout,*) "RSRH in sbc_cpl_rcv CC qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !                                                      ! atmosph. CO2 (ppm) !
       !                                                      ! ================== !
@@ -1421,4 +1428,5 @@
       ENDIF
       
+write(numout,*) "RSRH in sbc_cpl_rcv QQ qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !                                                      ! ========================= !
       IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
@@ -1452,4 +1460,5 @@
          ELSE                   ;   emp(:,:) =                              zemp(:,:)
          ENDIF
+write(numout,*) "RSRH in sbc_cpl_rcv VV qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
          !
          !                                                       ! non solar heat flux over the ocean (qns)
@@ -1488,4 +1497,5 @@
          !
       ENDIF
+write(numout,*) "RSRH in sbc_cpl_rcv ZZ qns, emp, sfx",kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !
       CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
@@ -1778,4 +1788,8 @@
       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i )
 
+write(numout,*) "RSRH in sbc_cpl_ice_flx AA qns, emp, sfx", glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx AA emp_tot, emp_ice, emp_oce", glob_sum(emp_tot(:,:)), glob_sum(emp_ice(:,:)), glob_sum(emp_oce(:,:));flush(numout)
+
+
       IF (sn_rcv_emp%clcat == 'yes') THEN
          CALL wrk_alloc( jpi, jpj, jpl, zevap_ice)
@@ -1798,9 +1812,16 @@
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con case"
          zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con casejpr_snow ",glob_sum(frcv(jpr_snow)%z3(:,:,1)); flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con casejpr_rain ",glob_sum(frcv(jpr_rain)%z3(:,:,1)); flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con casejpr_tevp ",glob_sum(frcv(jpr_tevp)%z3(:,:,1)); flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con casejpr_ievp ",glob_sum(frcv(jpr_ievp)%z3(:,:,1)); flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB con case picefr ",glob_sum(picefr(:,:)); flush(numout)
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB oce and ice case"
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
          zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
@@ -1808,4 +1829,8 @@
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
       END SELECT
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB what case?",TRIM( sn_rcv_emp%cldes )
+
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB emp_tot, emp_ice, emp_oce", glob_sum(emp_tot(:,:)), glob_sum(emp_ice(:,:)), glob_sum(emp_oce(:,:));flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx BB zemp_tot, zemp_ice", glob_sum(zemp_tot(:,:)), glob_sum(zemp_ice(:,:));flush(numout)
 
 #if defined key_lim3
@@ -1906,4 +1931,5 @@
 
 #endif
+write(numout,*) "RSRH in sbc_cpl_ice_flx LL emp_tot, emp_ice, emp_oce", glob_sum(emp_tot(:,:)), glob_sum(emp_ice(:,:)), glob_sum(emp_oce(:,:));flush(numout)
       ! outputs
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
@@ -1975,4 +2001,5 @@
       ! Heat content per unit mass of rain (J/kg)
       zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) ) 
+write(numout,*) "RSRH in sbc_cpl_ice_flx NN zcptsnw,zcptn", glob_sum(zcptsnw(:,:)), glob_sum(zcptn(:,:));flush(numout)
 
       ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
@@ -1991,4 +2018,5 @@
 !!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
 !!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice
+write(numout,*) "RSRH in sbc_cpl_ice_flx NN zsprecip,zcptrain", glob_sum(zsprecip(:,:)), glob_sum(zcptrain(:,:));flush(numout)
       
       ! --- total non solar flux (including evap/precip) --- !
@@ -1997,14 +2025,21 @@
       ! --- in case both coupled/forced are active, we must mix values --- ! 
       IF( ln_mixcpl ) THEN
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.1 zqns_tot,zqns_ice, zqns_oce", glob_sum(zqns_tot(:,:)), glob_sum(zqns_ice(:,:,:)), glob_sum(zqns_oce(:,:));flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.1 xcplmask", glob_sum(xcplmask(:,:,0));flush(numout)
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.1 xcplmask", glob_sum(Zmsk(:,:));flush(numout)
          qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
          qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.2 zqns_tot,zqns_ice, zqns_oce", glob_sum(zqns_tot(:,:)), glob_sum(zqns_ice(:,:,:)), glob_sum(zqns_oce(:,:));flush(numout)
          DO jl=1,jpl
             qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
             qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
          ENDDO
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.3 zqns_tot,zqns_ice, zqns_oce", glob_sum(zqns_tot(:,:)), glob_sum(zqns_ice(:,:,:)), glob_sum(zqns_oce(:,:));flush(numout)
          qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
          qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
          qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
+write(numout,*) "RSRH in sbc_cpl_ice_flx OO.4 zqns_tot,zqns_ice, zqns_oce", glob_sum(zqns_tot(:,:)), glob_sum(zqns_ice(:,:,:)), glob_sum(zqns_oce(:,:));flush(numout)
       ELSE
+write(numout,*) "RSRH in sbc_cpl_ice_flx PP zqns_tot,zqns_ice, zqns_oce", glob_sum(zqns_tot(:,:)), glob_sum(zqns_ice(:,:,:)), glob_sum(zqns_oce(:,:));flush(numout)
          qns_tot  (:,:  ) = zqns_tot  (:,:  )
          qns_oce  (:,:  ) = zqns_oce  (:,:  )
@@ -2015,4 +2050,5 @@
          qemp_ice (:,:  ) = zqemp_ice (:,:  )
       ENDIF
+write(numout,*) "RSRH in sbc_cpl_ice_flx PP qns_tot, qns_ice, qns_oce", glob_sum(qns_tot(:,:)), glob_sum(qns_ice(:,:)), glob_sum(qns_oce(:,:));flush(numout)
 
 #else
@@ -2049,4 +2085,5 @@
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
+write(numout,*) "RSRH in sbc_cpl_ice_flx MM qns, emp, sfx", glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
@@ -2202,4 +2239,5 @@
          CALL wrk_dealloc( jpi,jpj,1,zevap_ice)
       ENDIF
+write(numout,*) "RSRH in sbc_cpl_ice_flx ZZ qns, emp, sfx", glob_sum(qns(:,:)), glob_sum(emp(:,:)),glob_sum( sfx(:,:));flush(numout)
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 8892)
@@ -20,4 +20,5 @@
    !
    USE in_out_manager ! I/O manager
+   use lib_fortran
    USE prtctl         ! Print control
    USE iom            ! IOM library
@@ -59,8 +60,14 @@
       !
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
+write(numout,*) "RSRH sbcssm tsb", kt, glob_sum(tsb(:,:,:,jp_tem)); flush(numout)
+write(numout,*) "RSRH sbcssm tsn", kt, glob_sum(tsn(:,:,:,jp_tem)); flush(numout)
+write(numout,*) "RSRH sbcssm tsa", kt, glob_sum(tsa(:,:,:,jp_tem)); flush(numout)
+write(numout,*) "RSRH PE sbcssm tsb", kt, sum(tsb(:,:,:,jp_tem)); flush(numout)
+write(numout,*) "RSRH PE sbcssm tsn", kt, sum(tsn(:,:,:,jp_tem)); flush(numout)
+write(numout,*) "RSRH PE sbcssm tsa", kt, sum(tsa(:,:,:,jp_tem)); flush(numout)
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zts(ji,jj,jp_tem) = tsn(ji,jj,mikt(ji,jj),jp_tem)
-            zts(ji,jj,jp_sal) = tsn(ji,jj,mikt(ji,jj),jp_sal)
+            zts(ji,jj,jp_tem) = tsb(ji,jj,mikt(ji,jj),jp_tem)
+            zts(ji,jj,jp_sal) = tsb(ji,jj,mikt(ji,jj),jp_sal)
          END DO
       END DO
@@ -68,13 +75,30 @@
       IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
+write(numout,*) "RSRH sbcssm nnfsbc=1", kt; flush(numout)
+write(numout,*) "RSRH sbcssm ub 00", kt, glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm vb 00", kt, glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm ssu_m 00", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH sbcssm ssv_m 00", kt, glob_sum(ssv_m(:,:)); flush(numout)
+write(numout,*) "RSRH sbcssm sst_m 00", kt, glob_sum(sst_m(:,:)); flush(numout)
+
          ssu_m(:,:) = ub(:,:,1)
          ssv_m(:,:) = vb(:,:,1)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem)
+write(numout,*) "RSRH sbcssm ub 01", kt, glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm vb 01", kt, glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm ssu_m 01", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH sbcssm ssv_m 01", kt, glob_sum(ssv_m(:,:)); flush(numout)
+         IF( l_useCT )  THEN    
+   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+write(numout,*) "RSRH sbcssm sst_m AA", kt, glob_sum(sst_m(:,:)); flush(numout)
+         ELSE                    
+   sst_m(:,:) = zts(:,:,jp_tem)
+write(numout,*) "RSRH sbcssm sst_m BB", kt, glob_sum(sst_m(:,:)); flush(numout)
          ENDIF
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
+         IF( ln_apr_dyn ) THEN   
+   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    
+   ssh_m(:,:) = sshn(:,:)
          ENDIF
          !
@@ -84,7 +108,9 @@
          !
       ELSE
+write(numout,*) "RSRH sbcssm nnfsbc NOT 1", kt; flush(numout)
          !                                                ! ----------------------------------------------- !
          IF( kt == nit000 .AND. .NOT. l_ssm_mean ) THEN   !   Initialisation: 1st time-step, no input means !
-            !                                             ! ----------------------------------------------- !
+ write(numout,*) "RSRH sbcssm nnfsbc NOT 1 BB", kt; flush(numout)
+           !                                             ! ----------------------------------------------- !
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) 'sbc_ssm : mean fields initialised to instantaneous values'
@@ -107,4 +133,5 @@
             !                                             ! ---------------------------------------- !
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
+ write(numout,*) "RSRH sbcssm nnfsbc NOT 1 CC", kt; flush(numout)
             !                                             ! ---------------------------------------- !
             ssu_m(:,:) = 0._wp     ! reset to zero ocean mean sbc fields
@@ -119,8 +146,22 @@
          !                                                !        Cumulate at each time step        !
          !                                                ! ---------------------------------------- !
+write(numout,*) "RSRH sbcssm DD add ub vb", kt; flush(numout)
+write(numout,*) "RSRH sbcssm  ub DD", kt, glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm  vb DD", kt, glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH sbcssm  ssu_m DD", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH sbcssm  ssv_m DD", kt, glob_sum(ssv_m(:,:)); flush(numout)
+write(numout,*) "RSRH sbcssm  sst_m DD", kt, glob_sum(sst_m(:,:)); flush(numout)
          ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1)
          ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem)
+write(numout,*) "RSRH ub DD aft", kt, glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH vb DD aft", kt, glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH ssu_m DD aft", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH ssv_m DD aft", kt, glob_sum(ssv_m(:,:)); flush(numout)
+         IF( l_useCT )  THEN    
+ write(numout,*) "RSRH sbcssm EE update sst_m", kt; flush(numout)
+              sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+         ELSE                    
+ write(numout,*) "RSRH sbcssm FF update sst_m", kt; flush(numout)
+              sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem)
          ENDIF
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
@@ -145,5 +186,10 @@
             e3t_m(:,:) = e3t_m(:,:) * zcoef     ! mean vertical scale factor [m]
             frq_m(:,:) = frq_m(:,:) * zcoef     ! mean fraction of solar net radiation absorbed in the 1st T level [-]
-            !
+write(numout,*) "RSRH ub EE aft", kt, glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH vb EE aft", kt, glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH ssu_m EE aft", kt, glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH ssv_m EE aft", kt, glob_sum(ssv_m(:,:)); flush(numout)
+write(numout,*) "RSRH sst_m EE aft", kt, glob_sum(sst_m(:,:)); flush(numout)
+           !
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -235,5 +281,12 @@
          ENDIF
       ENDIF
-      !
+
+write(numout,*) "RSRH ub FF aft", glob_sum(ub(:,:,1)); flush(numout)
+write(numout,*) "RSRH vb FF aft",  glob_sum(vb(:,:,1)); flush(numout)
+write(numout,*) "RSRH ssu_m FF aft",  glob_sum(ssu_m(:,:)); flush(numout)
+write(numout,*) "RSRH ssv_m FF aft",  glob_sum(ssv_m(:,:)); flush(numout)
+
+
+     !
       IF( .NOT. l_ssm_mean ) THEN   ! default initialisation. needed by iceistate
          !
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8892)
@@ -27,4 +27,5 @@
    USE trd_oce        ! trends: ocean variables
    USE trdtra         ! trends manager: tracers 
+use lib_fortran
    !
    USE in_out_manager ! I/O manager
@@ -91,4 +92,7 @@
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
+
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn AA", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
+
       !
 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
@@ -119,4 +123,5 @@
          sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
       ENDIF
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn BB", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !                             !==  Now sbc tracer content fields  ==!
       DO jj = 2, jpj
@@ -136,4 +141,9 @@
          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
       ENDIF
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn CC", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
+write(numout,*) "RSRH tra_sbc qns, emp, sfx FF", kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)), glob_sum(sfx(:,:)); flush(numout)
+write(numout,*) "RSRH tra_sbc tsn FF", kt, glob_sum(tsn(:,:,1,1)); flush(numout)
+write(numout,*) "RSRH tra_sbc e3t_n FF", kt, glob_sum(e3t_n(:,:,1)); flush(numout)
+write(numout,*) "RSRH tra_sbc sbc_tsc_b, sbc_tsc FF", kt, glob_sum(sbc_tsc_b(:,:,1)), glob_sum(sbc_tsc(:,:,1)); flush(numout)
       !
       DO jn = 1, jpts               !==  update tracer trend  ==!
@@ -144,4 +154,5 @@
          END DO
       END DO
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn DD", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !                  
       IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
@@ -158,4 +169,5 @@
 !!gm         probably taken into account in r1_hisf_tbl : to be verified
       IF( ln_isf ) THEN
+write(numout,*) "RSRH tra_sbc risf_tsc_b, risf_tsc FF", kt, glob_sum(risf_tsc_b(:,:,1)), glob_sum(risf_tsc(:,:,1)); flush(numout)
          zfact = 0.5_wp
          DO jj = 2, jpj
@@ -183,4 +195,5 @@
          END DO
       END IF
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn EE", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !
       !----------------------------------------
@@ -189,4 +202,5 @@
       !
       IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
+write(numout,*) "RSRH tra_sbc rnf_tsc_b, rnf_tsc FF", kt, glob_sum(rnf_tsc_b(:,:,1)), glob_sum(rnf_tsc(:,:,1)); flush(numout)
          zfact = 0.5_wp
          DO jj = 2, jpj 
@@ -204,4 +218,5 @@
          END DO  
       ENDIF
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn FF", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
       IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
@@ -226,4 +241,5 @@
          END DO
       ENDIF
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn GG", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
@@ -239,4 +255,5 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
+write(numout,*) "RSRH tra_sbc tsb, tsa, tsn HH", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !
    END SUBROUTINE tra_sbc
Index: /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8891)
+++ /branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8892)
@@ -38,4 +38,8 @@
    !
    USE iom              ! xIOs server
+use in_out_manager
+use lib_fortran
+use oce
+
 
    IMPLICIT NONE
@@ -116,4 +120,6 @@
       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
 
+write(numout,*) "RSRH step tsb, tsa, tsn AA", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Ocean physics update                (ua, va, tsa used as workspace)
@@ -124,4 +130,5 @@
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+write(numout,*) "RSRH step tsb, tsa, tsn BB", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
       !
@@ -139,4 +146,5 @@
          avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
       ENDIF
+write(numout,*) "RSRH step tsb, tsa, tsn CC", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
@@ -151,4 +159,5 @@
                          CALL zdf_mxl( kstp )         ! mixed layer depth
 
+write(numout,*) "RSRH step tsb, tsa, tsn DD", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
                                                       ! write TKE or GLS information in the restart file
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
@@ -173,4 +182,5 @@
          ENDIF
       ENDIF
+write(numout,*) "RSRH step tsb, tsa, tsn EE", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !                                                                   ! eddy diffusivity coeff.
       IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )       !       and/or eiv coeff.
@@ -197,4 +207,5 @@
                &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 !!jc: fs simplification
+write(numout,*) "RSRH step tsb, tsa, tsn FF", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
                             
                          ua(:,:,:) = 0._wp            ! set dynamics trends to zero
@@ -220,4 +231,5 @@
                             CALL wzv        ( kstp )              ! now cross-level velocity 
       ENDIF
+write(numout,*) "RSRH step tsb, tsa, tsn GG", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
                          CALL dyn_bfr       ( kstp )  ! bottom friction
@@ -237,9 +249,13 @@
       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
                          CALL dia_ar5( kstp )         ! ar5 diag
+write(numout,*) "RSRH step tsb, tsa, tsn HH", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+write(numout,*) "RSRH step tsb, tsa, tsn HH.1", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
                          CALL dia_wri( kstp )         ! ocean model: outputs
+write(numout,*) "RSRH step tsb, tsa, tsn HH.2", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       !
       IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output
       
+write(numout,*) "RSRH step tsb, tsa, tsn HH.3", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 #if defined key_top
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -247,4 +263,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL trc_stp       ( kstp )  ! time-stepping
+write(numout,*) "RSRH step tsb, tsa, tsn HH.4", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 #endif
 
@@ -257,14 +274,21 @@
          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
                          CALL tra_sbc       ( kstp )  ! surface boundary condition
+write(numout,*) "RSRH step tsb, tsa, tsn HH.5", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
+write(numout,*) "RSRH step tsb, tsa, tsn HH.6", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
+write(numout,*) "RSRH step tsb, tsa, tsn HH.7", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( lk_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
+write(numout,*) "RSRH step tsb, tsa, tsn HH.8", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends
+write(numout,*) "RSRH step tsb, tsa, tsn HH.9", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
       IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
+write(numout,*) "RSRH step tsb, tsa, tsn HH.10", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 #if defined key_agrif
       IF(.NOT. Agrif_Root())  & 
                &         CALL Agrif_Sponge_tra        ! tracers sponge
 #endif
-                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection
+ write(numout,*) "RSRH step tsb, tsa, tsn II", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
+                        CALL tra_adv       ( kstp )  ! horizontal & vertical advection
                          CALL tra_ldf       ( kstp )  ! lateral mixing
 
@@ -298,4 +322,5 @@
       !
       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
+ write(numout,*) "RSRH step tsb, tsa, tsn JJ", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
 !!gm : This does not only concern the dynamics ==>>> add a new title
@@ -328,5 +353,4 @@
                          CALL ctl_stop( 'step: indic < 0' )
                          CALL dia_wri_state( 'output.abort', kstp )
-                               CALL ctl_stop('STOP','NEMO failure in stp')
       ENDIF
       IF( kstp == nit000 ) THEN
@@ -336,9 +360,11 @@
             &    CALL FLUSH    ( numoni )     ! flush output namelist ice (if exist)
       ENDIF
+ write(numout,*) "RSRH step tsb, tsa, tsn KK", kstp, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+!!gm why lk_oasis and not lk_cpl ????
+      !IF( lk_oasis   )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
