Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icestp.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icestp.F90	(revision 11355)
@@ -145,4 +145,10 @@
 #endif
                                         CALL store_fields             ! Store now ice values
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_sbc_tau = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+         WRITE(9419,*) 'qlead(48,31) = ',qlead(48,31)
+      ENDIF
          !
          !------------------------------------------------!
@@ -159,4 +165,11 @@
                                         CALL ice_rst_opn( kt )        ! Open Ice restart file (if necessary) 
          !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_dyn = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+         WRITE(9419,*) 'qlead(48,31) = ',qlead(48,31)
+      ENDIF
+
          IF( ln_icedyn .AND. .NOT.lk_c1d )   &
             &                           CALL ice_dyn( kt )            ! -- Ice dynamics
@@ -165,4 +178,11 @@
          IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
          !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_var_glo2eqv = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+         WRITE(9419,*) 'qlead(48,31) = ',qlead(48,31)
+      ENDIF
+
          !                          !==  previous lead fraction and ice volume for flux calculations
                                         CALL ice_var_glo2eqv          ! h_i and h_s for ice albedo calculation
@@ -183,17 +203,54 @@
          !    qprec_ice, qevap_ice
          !------------------------------------------------------!
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_sbc_flx = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+         WRITE(9419,*) 'qlead(48,31) = ',qlead(48,31)
+      ENDIF
+
                                         CALL ice_sbc_flx( kt, ksbc )
          !----------------------------!
          ! --- ice thermodynamics --- !
          !----------------------------!
-         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_thd = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+         WRITE(9419,*) 'qlead(48,31) = ',qlead(48,31)
+      ENDIF
+
+         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics  
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_cor = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+      ENDIF
+    
          !
                                         CALL ice_cor( kt , 2 )        ! -- Corrections
          !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_var_glo2eqv = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+      ENDIF
+
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
          !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_update_flx = ',emp(48,31)
+         WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+      ENDIF
+
                                         CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
          !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp before ice_dia = ',emp(48,31)
+      ENDIF
+
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
          !
@@ -203,4 +260,9 @@
          !
          IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp ice_stp at end = ',emp(48,31)
+      ENDIF
+
          !
       ENDIF   ! End sea-ice time step only
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd.F90	(revision 11355)
@@ -329,4 +329,5 @@
       !
       INTEGER ::   jk   ! dummy loop indices
+      INTEGER ::   jn , jid, jjd
       !!-----------------------------------------------------------------------
       !
@@ -335,4 +336,17 @@
       CASE( 1 )            !==  from 2D to 1D  ==!
          !                 !---------------------!
+
+         DO jn = 1, npti
+           jid        = MOD( nptidx(jn) - 1 , jpi ) + 1
+           jjd        =    ( nptidx(jn) - 1 ) / jpi + 1
+           IF( jn == 54 ) THEN
+             IF( narea == 68) THEN
+               WRITE(968,*) 'in ice_thd_1d2d: ice_cat=',kl
+               WRITE(968,*) 'in ice_thd_1d2d: jid, jjd for point 54 is ',jid,jjd
+               WRITE(968,*) 'a_i for this point = ',a_i(jid,jjd)
+             END IF
+           END IF
+         END DO
+
          CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_dh.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_dh.F90	(revision 11355)
@@ -105,4 +105,5 @@
 
       INTEGER  ::   num_iter_max      ! Heat conservation 
+      INTEGER, DIMENSION(1)  :: imax
       !!------------------------------------------------------------------
 
@@ -182,5 +183,21 @@
             END IF
          END DO
-      END DO         
+      END DO
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max wfx_snw_sum_1d icethd_dh after inte5rnal melting = ',MAXVAL(  ABS( wfx_snw_sum_1d(:) )  )
+         imax = MAXLOC( ABS( wfx_snw_sum_1d(:) ) )
+         WRITE(9419,*) 'at location = ',imax
+         WRITE(9419,*) 'zh_s = ',zh_s(imax,:)
+         WRITE(9419,*) 'a_i_1d = ',a_i_1d(imax)
+      ENDIF
+
+      IF(narea == 68) THEN
+         WRITE(968,*) 'max wfx_snw_sum_1d icethd_dh after inte5rnal melting = ',MAXVAL(  ABS( wfx_snw_sum_1d(:) )  )
+         imax = MAXLOC( ABS( wfx_snw_sum_1d(:) ) )
+         WRITE(968,*) 'at location = ',imax
+         WRITE(968,*) 'zh_s = ',zh_s(imax,:)
+         WRITE(968,*) 'a_i_1d = ',a_i_1d(imax)
+      ENDIF
 
       ! Snow precipitation
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_do.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/icethd_do.F90	(revision 11355)
@@ -207,4 +207,12 @@
       END DO
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'npti in icetd_do = ',npti
+         WRITE(9419,*) 'max qlead = ',MAXVAL(  qlead(:,:) )
+         WRITE(9419,*) 'min qlead = ',MINVAL(  qlead(:,:) )
+         WRITE(9419,*) 'qlead(48,31) = ', qlead(48,31)
+         WRITE(9419,*) 'tau_icebfr(48,31) = ',tau_icebfr(48,31)
+      ENDIF
+
       ! Move from 2-D to 1-D vectors
       IF ( npti > 0 ) THEN
@@ -230,4 +238,9 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
+
+         IF(narea == 419) THEN
+            WRITE(9419,*) 'wfx_opw(48,31) = ',wfx_opw(48,31)
+            WRITE(9419,*) 'max wfx_opw_1d icethd_do after to 1d = ',MAXVAL(  ABS( wfx_opw_1d(:) )  )
+         ENDIF
 
          ! Convert units for ice internal energy
@@ -295,4 +308,14 @@
             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_rdtice
          END DO
+
+         IF(narea == 419) THEN
+            WRITE(9419,*) 'max wfx_opw_1d icethd_do after calculation = ',MAXVAL(  ABS( wfx_opw_1d(:) )  )
+            WRITE(9419,*) 'max zv_newice icethd_do after calculation = ',MAXVAL(  ABS( zv_newice(:) )  )
+            WRITE(9419,*) 'rhoi = ',rhoi
+            WRITE(9419,*) 'r1_rdtice = ',r1_rdtice
+            WRITE(9419,*) 'max qlead_1d icethd_do after calculation = ',MAXVAL(  ABS( qlead_1d(:) )  )
+            WRITE(9419,*) 'max ze_newice icethd_do after calculation = ',MAXVAL(  ABS( ze_newice(:) )  )
+            WRITE(9419,*) 'max t_bo_1d icethd_do after calculation = ',MAXVAL(  ABS( t_bo_1d(:) )  )
+         ENDIF
          
          zv_frazb(1:npti) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/iceupdate.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/ICE/iceupdate.F90	(revision 11355)
@@ -96,4 +96,6 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
+      INTEGER, DIMENSION(2)  ::   imax                  ! min/max loc indices
+
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('ice_update')
@@ -112,4 +114,52 @@
          qemp_ice   (:,:)   = 0._wp
          qevap_ice  (:,:,:) = 0._wp
+      ENDIF
+
+      IF(narea == 419) THEN
+         ! emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)
+         ! wfx_ice = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
+         !      &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
+         ! wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
+         WRITE(9419,*) 'max emp_oce iceupdate before loop = ',MAXVAL(  ABS( emp_oce(:,:) )  )
+         WRITE(9419,*) 'max wfx_err_sub iceupdate before loop = ',MAXVAL(  ABS( wfx_err_sub(:,:) )  )
+         WRITE(9419,*) 'max wfx_bog iceupdate before loop = ',MAXVAL(  ABS( wfx_bog(:,:) )  )
+         WRITE(9419,*) 'max wfx_bom iceupdate before loop = ',MAXVAL(  ABS( wfx_bom(:,:) )  )
+         WRITE(9419,*) 'max wfx_sum iceupdate before loop = ',MAXVAL(  ABS( wfx_sum(:,:) )  )
+         WRITE(9419,*) 'max wfx_sni iceupdate before loop = ',MAXVAL(  ABS( wfx_sni(:,:) )  )
+         WRITE(9419,*) 'max wfx_opw iceupdate before loop = ',MAXVAL(  ABS( wfx_opw(:,:) )  )
+         WRITE(9419,*) 'max wfx_dyn iceupdate before loop = ',MAXVAL(  ABS( wfx_dyn(:,:) )  )
+         WRITE(9419,*) 'max wfx_res iceupdate before loop = ',MAXVAL(  ABS( wfx_res(:,:) )  )
+         WRITE(9419,*) 'max wfx_lam iceupdate before loop = ',MAXVAL(  ABS( wfx_lam(:,:) )  )
+         WRITE(9419,*) 'max wfx_pnd iceupdate before loop = ',MAXVAL(  ABS( wfx_pnd(:,:) )  )
+         WRITE(9419,*) 'max wfx_snw_sni iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_sni(:,:) )  )
+         WRITE(9419,*) 'max wfx_snw_dyn iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_dyn(:,:) )  )
+         WRITE(9419,*) 'max wfx_snw_sum iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_sum(:,:) )  )
+         imax(:) = MAXLOC( ABS( wfx_snw_sum(:,:) ) )
+         WRITE(9419,*) 'wfx_snw_sum at location = ',imax
+         WRITE(9419,*) 'wfx_snw_sum(imax) = ',wfx_snw_sum(imax(1),imax(2))
+      ENDIF
+
+      IF(narea == 68) THEN
+         ! emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)
+         ! wfx_ice = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
+         !      &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
+         ! wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
+         WRITE(968,*) 'max emp_oce iceupdate before loop = ',MAXVAL(  ABS( emp_oce(:,:) )  )
+         WRITE(968,*) 'max wfx_err_sub iceupdate before loop = ',MAXVAL(  ABS( wfx_err_sub(:,:) )  )
+         WRITE(968,*) 'max wfx_bog iceupdate before loop = ',MAXVAL(  ABS( wfx_bog(:,:) )  )
+         WRITE(968,*) 'max wfx_bom iceupdate before loop = ',MAXVAL(  ABS( wfx_bom(:,:) )  )
+         WRITE(968,*) 'max wfx_sum iceupdate before loop = ',MAXVAL(  ABS( wfx_sum(:,:) )  )
+         WRITE(968,*) 'max wfx_sni iceupdate before loop = ',MAXVAL(  ABS( wfx_sni(:,:) )  )
+         WRITE(968,*) 'max wfx_opw iceupdate before loop = ',MAXVAL(  ABS( wfx_opw(:,:) )  )
+         WRITE(968,*) 'max wfx_dyn iceupdate before loop = ',MAXVAL(  ABS( wfx_dyn(:,:) )  )
+         WRITE(968,*) 'max wfx_res iceupdate before loop = ',MAXVAL(  ABS( wfx_res(:,:) )  )
+         WRITE(968,*) 'max wfx_lam iceupdate before loop = ',MAXVAL(  ABS( wfx_lam(:,:) )  )
+         WRITE(968,*) 'max wfx_pnd iceupdate before loop = ',MAXVAL(  ABS( wfx_pnd(:,:) )  )
+         WRITE(968,*) 'max wfx_snw_sni iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_sni(:,:) )  )
+         WRITE(968,*) 'max wfx_snw_dyn iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_dyn(:,:) )  )
+         WRITE(968,*) 'max wfx_snw_sum iceupdate before loop = ',MAXVAL(  ABS( wfx_snw_sum(:,:) )  )
+         imax(:) = MAXLOC( ABS( wfx_snw_sum(:,:) ) )
+         WRITE(968,*) 'wfx_snw_sum at location = ',imax
+         WRITE(968,*) 'wfx_snw_sum(imax) = ',wfx_snw_sum(imax(1),imax(2))
       ENDIF
       
@@ -177,4 +227,10 @@
          END DO
       END DO
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp iceupdate after loop = ',MAXVAL(  ABS( emp(:,:) )  )
+         imax(:) = MAXLOC( ABS( emp(:,:) ) )
+         WRITE(9419,*) 'at location = ',imax
+      ENDIF
 
       ! Storing the transmitted variables
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/DYN/sshwzv.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/DYN/sshwzv.F90	(revision 11355)
@@ -103,6 +103,20 @@
       ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha ssh_nxt before maths = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max sshb ssh_nxt before maths = ',MAXVAL(  ABS( sshb(:,:) )  )
+         WRITE(9419,*) 'max emp_b ssh_nxt before maths = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp ssh_nxt before maths = ',MAXVAL(  ABS( emp(:,:) )  )
+         WRITE(9419,*) 'max zhdiv ssh_nxt before maths = ',MAXVAL(  ABS( zhdiv(:,:) )  )
+         WRITE(9419,*) 'max ssmask ssh_nxt before maths = ',MAXVAL(  ABS( ssmask(:,:) )  )
+      ENDIF
+
       ! 
       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha ssh_nxt after maths = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
       !
 #if defined key_agrif
@@ -116,4 +130,9 @@
          ENDIF
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha ssh_nxt after bdy_ssh = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       !                                           !------------------------------!
       !                                           !           outputs            !
@@ -122,4 +141,9 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask )
       !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha ssh_nxt after prt_ctl = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_timing )   CALL timing_stop('ssh_nxt')
       !
@@ -246,4 +270,10 @@
       ENDIF
       !              !==  Euler time-stepping: no filter, just swap  ==!
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn in sshwzv = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha in sshwzv = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF ( neuler == 0 .AND. kt == nit000 ) THEN
          sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbccpl.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbccpl.F90	(revision 11355)
@@ -1378,4 +1378,9 @@
          CALL iom_put( 'frq_m', frq_m )
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'k_ice, kt, k_fsbc = ',k_ice, kt, k_fsbc
+         WRITE(9419,*) 'MOD( kt-1, k_fsbc ) = ',MOD( kt-1, k_fsbc )
+      ENDIF
       
       !                                                      ! ========================= !
@@ -1396,8 +1401,18 @@
             zemp(:,:) = 0._wp
          ENDIF
+
+         IF(narea == 419) THEN
+            WRITE(9419,*) 'max zemp after coupling = ',MAXVAL(  ABS( zemp(:,:) )  )
+         ENDIF
+
          !
          !                                                        ! runoffs and calving (added in emp)
          IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
          IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
+
+         IF(narea == 419) THEN
+            WRITE(9419,*) 'max zemp after runoff = ',MAXVAL(  ABS( zemp(:,:) )  )
+            WRITE(9419,*) 'max emp after runoff = ',MAXVAL(  ABS( emp(:,:) )  )
+         ENDIF
  
          IF( srcv(jpr_icb)%laction )  THEN 
@@ -1410,4 +1425,9 @@
          ELSE                   ;   emp(:,:) =                              zemp(:,:)
          ENDIF
+
+         IF(narea == 419) THEN
+            WRITE(9419,*) 'max emp after conversion from zemp = ',MAXVAL(  ABS( emp(:,:) )  )
+         ENDIF
+
          !
          !                                                       ! non solar heat flux over the ocean (qns)
@@ -1733,4 +1753,14 @@
       ENDIF
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'emp_tot(48,31) = ',emp_tot(48,31)
+         WRITE(9419,*) 'emp_ice(48,31) = ',emp_ice(48,31)
+         WRITE(9419,*) 'emp_oce(48,31) = ',emp_oce(48,31)
+         WRITE(9419,*) 'sprecip(48,31) = ',sprecip(48,31)
+         WRITE(9419,*) 'tprecip(48,31) = ',tprecip(48,31)
+         WRITE(9419,*) 'evap_ice(48,31) = ',evap_ice(48,31,:)
+         WRITE(9419,*) 'devap_ice(48,31) = ',devap_ice(48,31,:)
+      ENDIF
+
 #else
       zsnw(:,:) = picefr(:,:)
@@ -1874,4 +1904,14 @@
       ENDIF
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'qns_tot(48,31) = ',qns_tot(48,31)
+         WRITE(9419,*) 'qns_oce(48,31) = ',qns_oce(48,31)
+         WRITE(9419,*) 'qns_ice(48,31) = ',qns_ice(48,31,:)
+         WRITE(9419,*) 'qevap_ice(48,31) = ',qevap_ice(48,31,:)
+         WRITE(9419,*) 'qprec_ice(48,31) = ',qprec_ice(48,31)
+         WRITE(9419,*) 'qemp_oce(48,31) = ',qemp_oce(48,31)
+         WRITE(9419,*) 'qemp_ice(48,31) = ',qemp_ice(48,31)
+      ENDIF
+
 #else
       zcptsnw (:,:) = zcptn(:,:)
@@ -1914,4 +1954,5 @@
       CASE( 'oce only' )
          zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
+         zqsr_ice(:,:,:) = 0.0
       CASE( 'conservative' )
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1939,4 +1980,9 @@
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             END DO
+
+            IF(narea == 419) THEN
+               WRITE(9419,*) 'zqsr_ice(48,31) oce and ice = ',zqsr_ice(48,31,:)
+            ENDIF
+
          ENDIF
       CASE( 'mixed oce-ice' )
@@ -1980,4 +2026,10 @@
       ENDIF
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'qsr_tot(48,31) = ',qsr_tot(48,31)
+         WRITE(9419,*) 'qsr_ice(48,31) = ',qsr_ice(48,31,:)
+      ENDIF
+
+
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
@@ -1999,6 +2051,11 @@
          ENDDO
       ELSE
-         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
-      ENDIF
+         dqns_ice(:,:,:) = 0.0
+      ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'dqns_ice(48,31) = ',dqns_ice(48,31,:)
+      ENDIF
+
 
 #if defined key_si3      
@@ -2010,4 +2067,10 @@
          qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
       END SELECT
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'qml_ice(48,31) = ',qml_ice(48,31,:)
+         WRITE(9419,*) 'qcn_ice(48,31) = ',qcn_ice(48,31,:)
+      ENDIF
+
       !
       !                                                      ! ========================= !
@@ -2030,4 +2093,9 @@
          !
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'qtr_ice_top(48,31) = ',qtr_ice_top(48,31,:)
+      ENDIF
+
       !
 #endif
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbcmod.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/SBC/sbcmod.F90	(revision 11355)
@@ -414,4 +414,9 @@
       IF( ln_wave     )   CALL sbc_wave( kt )            ! surface waves
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'nsbc, jp_purecpl = ',nsbc, jp_purecpl
+         WRITE(9419,*) 'max emp before sbc_cpl_rcv = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       !
       !                                            !==  sbc formulation  ==!
@@ -435,4 +440,8 @@
       !
       !                                            !==  Misc. Options  ==!
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'nn_ice = ',nn_ice
+         WRITE(9419,*) 'max emp before ice_stp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
       !
       SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
@@ -444,4 +453,8 @@
       END SELECT
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after ice_stp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       IF( ln_icebergs    )   THEN
                                      CALL icb_stp( kt )           ! compute icebergs
@@ -455,4 +468,8 @@
       IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after sbc_rnf = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
 
@@ -461,5 +478,9 @@
       ! Special treatment of freshwater fluxes over closed seas in the model domain
       ! Should not be run if ln_diurnal_only
-      IF( l_sbc_clo .AND. (.NOT. ln_diurnal_only) )   CALL sbc_clo( kt )   
+      IF( l_sbc_clo .AND. (.NOT. ln_diurnal_only) )   CALL sbc_clo( kt )  
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after sbc_clo = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF 
 
 !!$!RBbug do not understand why see ticket 667
@@ -490,4 +511,9 @@
          END WHERE
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after north fold adjustments = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF 
+
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
@@ -517,4 +543,9 @@
          ENDIF
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after setting forcing fields = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF 
+
       !                                                ! ---------------------------------------- !
       IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
@@ -552,4 +583,9 @@
       CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
       CALL iom_put( "vtau", vtau )   ! j-wind stress
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max emp after outputs = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF 
+
       !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/nemogcm.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/nemogcm.F90	(revision 11355)
@@ -150,4 +150,14 @@
       CALL mpp_max( 'nemogcm', nstop )
 
+      IF(narea == 419) THEN
+         OPEN(9419, file='crash_stats_419.out', status='replace')
+         WRITE(9419,*) 'Opening output file'
+      ENDIF
+
+      IF(narea == 68) THEN
+         OPEN(968, file='crash_stats_68.out', status='replace')
+         WRITE(968,*) 'Opening output file'
+      ENDIF
+
       IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
 
@@ -156,4 +166,15 @@
       !                            !-----------------------!
       istp = nit000
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn nemogcm before any steps = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max emp_b = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
+      ! Initialise emp to zero
+      emp(:,:) = 0.0
+
       !
 #if defined key_c1d
@@ -218,4 +239,6 @@
       !                            !------------------------!
       IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA
+      CLOSE(unit=9419)
+      CLOSE(unit=968)
       !
       IF( nstop /= 0 .AND. lwp ) THEN        ! error print
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/step.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/step.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/step.F90	(revision 11355)
@@ -93,4 +93,12 @@
       !
       IF( ln_timing )   CALL timing_start('stp')
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before update IO = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max emp_b = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -107,4 +115,11 @@
       IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before external forcing = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max emp_b = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
@@ -115,4 +130,7 @@
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
 
+
+
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
@@ -120,4 +138,11 @@
       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before ocean physics = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max emp_b = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -155,4 +180,11 @@
       IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff. 
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before ocean dynamics = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+         WRITE(9419,*) 'max emp_b = ',MAXVAL(  ABS( emp_b(:,:) )  )
+         WRITE(9419,*) 'max emp = ',MAXVAL(  ABS( emp(:,:) )  )
+      ENDIF
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       !  Ocean dynamics : hdiv, ssh, e3, u, v, w
@@ -160,8 +192,34 @@
 
                             CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dom_vvl_sf_nxt = ',MAXVAL(  ABS( ssha(:,:) )  )
+         CALL FLUSH(9419)
+      ENDIF
+
+      IF(narea == 68) THEN
+         WRITE(968,*) 'max ssha step before dom_vvl_sf_nxt = ',MAXVAL(  ABS( ssha(:,:) )  )
+         CALL FLUSH(968)
+      ENDIF
+
+
       IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before wzv = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                             CALL wzv           ( kstp )  ! now cross-level velocity 
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before wAimp = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning
                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before simplification = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                             
 !!jc: fs simplification
@@ -182,4 +240,9 @@
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
                &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before bdy_dyn3d_dmp = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends
 #if defined key_agrif
@@ -187,10 +250,45 @@
                &         CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dom_adv = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_adv       ( kstp )  ! advection (vector or flux form)
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_vor = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_ldft = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_osm = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_hpg = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_spg = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
+
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before div_hor = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
 
                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well
@@ -201,6 +299,16 @@
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max ssha step before dyn_zdf = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       
                          CALL dyn_zdf       ( kstp )  ! vertical diffusion
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before cool skin = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -208,4 +316,9 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
       IF ( ln_diurnal )  CALL stp_diurnal( kstp )
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before diagnostics = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
       
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -230,4 +343,9 @@
 #endif
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before active tracers = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Active tracers                              
@@ -237,5 +355,17 @@
       IF(  lk_asminc .AND. ln_asmiau .AND. &
          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before tra_sbc = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL tra_sbc       ( kstp )  ! surface boundary condition
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step after tra_sbc = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
       IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
@@ -258,4 +388,9 @@
                          CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields
       IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before boundary conditions = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -277,8 +412,32 @@
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
                          CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before dyn_nxt = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before ssh_swp = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before dom_vvl_sf_swp = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before dia_hsb = ',MAXVAL(  ABS( sshn(:,:) )  )
+         WRITE(9419,*) 'max ssha = ',MAXVAL(  ABS( ssha(:,:) )  )
+      ENDIF
+
       IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics
 
@@ -291,4 +450,8 @@
       IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
 
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before agrif = ',MAXVAL(  ABS( sshn(:,:) )  )
+      ENDIF
+
 #if defined key_agrif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -304,4 +467,9 @@
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn step before stp_ctl = ',MAXVAL(  ABS( sshn(:,:) )  )
+      ENDIF
+
                          CALL stp_ctl      ( kstp, indic )
                          
Index: /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/stpctl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/stpctl.F90	(revision 11354)
+++ /NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/stpctl.F90	(revision 11355)
@@ -108,4 +108,8 @@
          REWIND( numstp )
       ENDIF
+
+      IF(narea == 419) THEN
+         WRITE(9419,*) 'max sshn stpctl before crash = ',MAXVAL(  ABS( sshn(:,:) )  )
+      ENDIF
       !
       !                                   !==  test of extrema  ==!
@@ -174,4 +178,7 @@
          
          CALL dia_wri_state( 'output.abort' )     ! create an output.abort file
+
+         CLOSE(unit=9419)
+         CLOSE(unit=968)
          
          IF( .NOT. ln_ctl ) THEN
