Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12475)
@@ -221,4 +221,5 @@
           <field id="dmibog"       long_name="sea-ice mass change through basal growth"                standard_name="tendency_of_sea_ice_amount_due_to_congelation_ice_accumulation"          unit="kg/m2/s" />
           <field id="dmisni"       long_name="sea-ice mass change through snow-to-ice conversion"      standard_name="tendency_of_sea_ice_amount_due_to_snow_conversion"                       unit="kg/m2/s" />
+          <field id="dmisnp"       long_name="sea-ice mass change through snow melting into ponds"     standard_name="tendency_of_sea_ice_amount_due_to_snow_melt_into_ponds"                  unit="kg/m2/s" />
           <field id="dmtsub"       long_name="snow and sea ice mass change through sublimation"        standard_name="tendency_of_snow_and_ice_amount_due_to_sublimation"                      unit="kg/m2/s" />
           <field id="dmssub"       long_name="sea-ice mass change through evaporation and sublimation" standard_name="water_evaporation_flux"                                                  unit="kg/m2/s" />
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90	(revision 12475)
@@ -245,4 +245,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: mass flux from ice-ocean mass exchange                   [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni     !: mass flux from snow ice growth component of wfx_ice      [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snp     !: mass flux from snow melting into melt ponds              [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw     !: mass flux from lateral ice growth component of wfx_ice   [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog     !: mass flux from bottom ice growth component of wfx_ice    [kg.m-2.s-1]
@@ -415,5 +416,5 @@
          &      wfx_bog    (jpi,jpj) , wfx_dyn   (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,           &
          &      wfx_res    (jpi,jpj) , wfx_sni   (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,           &
-         &      afx_tot    (jpi,jpj) , rn_amax_2d(jpi,jpj) ,                                                 &
+         &      afx_tot    (jpi,jpj) , rn_amax_2d(jpi,jpj) , wfx_snp(jpi,jpj) ,                              &
          &      qsb_ice_bot(jpi,jpj) , qlead     (jpi,jpj) ,                                                 &
          &      sfx_res    (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90	(revision 12475)
@@ -77,4 +77,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sum_1d  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sni_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snp_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_opw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d 
@@ -203,4 +204,5 @@
          &      fhld_1d       (jpij) , wfx_sub_1d    (jpij) , wfx_bog_1d    (jpij) , wfx_bom_1d    (jpij) ,  &
          &      wfx_sum_1d    (jpij) , wfx_sni_1d    (jpij) , wfx_opw_1d    (jpij) , wfx_res_1d    (jpij) ,  &
+         &      wfx_snp_1d    (jpij) ,                                                                       &
          &      wfx_snw_sub_1d(jpij) , wfx_snw_dyn_1d(jpij) , wfx_ice_sub_1d(jpij) , wfx_err_sub_1d(jpij) ,  &
          &      wfx_lam_1d    (jpij) , wfx_dyn_1d    (jpij) , wfx_pnd_1d    (jpij) , dqns_ice_1d   (jpij) , evap_ice_1d (jpij) , &
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn.F90	(revision 12475)
@@ -105,5 +105,18 @@
       !
 
-      write(numout,*)'ice_dyn 0: e_i = ',e_i (42,26,1,1)
+      write(numout,*)'ice_dyn 0:sv_i, e_i, t_i = ',sv_i(3,4,1), ' ', e_i(3,4,1,1), ' ', t_i(3,4,1,1)
+      write(numout,*)'ice_dyn 0:a_i, h_i, v_i = ',a_i(3,4,1), ' ', h_i(3,4,1), ' ', v_i(3,4,1)
+      write(numout,*)'ice_dyn 0:s_i = ',s_i(3,4,1)
+      write(numout,*)'ice_dyn 0: MEAN(sv_i) = ',SUM(sv_i)/SIZE(sv_i)
+      write(numout,*)'ice_dyn 0: MEAN(e_i) = ',SUM(e_i)/SIZE(e_i)
+      write(numout,*)'ice_dyn 0: MEAN(t_i) = ',SUM(t_i)/SIZE(t_i)
+      write(numout,*)'ice_dyn 0: MEAN(a_i) = ',SUM(a_i)/SIZE(a_i)
+      write(numout,*)'ice_dyn 0: MEAN(h_i) = ',SUM(h_i)/SIZE(h_i)
+      write(numout,*)'ice_dyn 0: MEAN(v_i) = ',SUM(v_i)/SIZE(v_i)
+      write(numout,*)'ice_dyn 0: MEAN(s_i) = ',SUM(s_i)/SIZE(s_i)
+      write(numout,*)'ice_dyn 0: MEAN(u_ice) = ',SUM(u_ice)/SIZE(u_ice)
+      write(numout,*)'ice_dyn 0: MEAN(e2u) = ',SUM(e2u)/SIZE(e2u)
+      write(numout,*)'ice_dyn 0: u_ice, e2u = ',u_ice(3,4), e2u(3,4)
+      
       !
       SELECT CASE( nice_dyn )          !-- Set which dynamics is running
@@ -111,7 +124,11 @@
       CASE ( np_dynALL )           !==  all dynamical processes  ==!
          !
-         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology 
+          write(numout,*)'ice_dyn 1:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
+          write(numout,*)'ice_dyn 1: u_ice, e2u = ',u_ice(3,4), e2u(3,4)
          CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+          write(numout,*)'ice_dyn 2:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
          CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting 
+          write(numout,*)'ice_dyn 3:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
          CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
          !
@@ -150,5 +167,5 @@
       END SELECT
 
-      write(numout,*)'ice_dyn 1: e_i = ',e_i (42,26,1,1)
+      write(numout,*)'ice_dyn 4: e_i = ',e_i (42,26,1,1)
       !
       !
@@ -174,5 +191,5 @@
          !
       ENDIF
-      write(numout,*)'ice_dyn 2: e_i = ',e_i (42,26,1,1)
+      write(numout,*)'ice_dyn 5: e_i = ',e_i (42,26,1,1)
       !
       ! controls
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_umx.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_umx.F90	(revision 12475)
@@ -105,4 +105,6 @@
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
+
+      write(numout,*)'ice_dyn_adv_umx 1:psv_i, pe_s = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1)
       !
       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
@@ -142,7 +144,12 @@
       zdt = rdt_ice / REAL(icycle)
 
+      write(numout,*)'ice_dyn_adv_umx 2: pu_ice, e2u = ',pu_ice(3,4), e2u(3,4)
+
       ! --- transport --- !
       zudy(:,:) = pu_ice(:,:) * e2u(:,:)
       zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
+
+      write(numout,*)'ice_dyn_adv_umx 2: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
+
       !
       ! setup transport for each ice cat 
@@ -151,4 +158,7 @@
          zv_cat(:,:,jl) = zvdx(:,:)
       END DO
+
+      write(numout,*)'ice_dyn_adv_umx 3: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
+
       !
       ! --- define velocity for advection: u*grad(H) --- !
@@ -166,4 +176,6 @@
          END DO
       END DO
+
+      write(numout,*)'ice_dyn_adv_umx 4: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
 
       !---------------!
@@ -209,4 +221,7 @@
          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
             &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
+
+         write(numout,*)'ice_dyn_adv_umx 5: zudy, zvdx, zu_cat, zv_cat, zcu_box, zcv_box = ',zudy(3,4), zvdx(3,4), zu_cat(3,4,1), zv_cat(3,4,1), zcu_box(3,4), zcv_box(3,4)
+
          !
          !                             ! --------------------------------- !
@@ -227,4 +242,5 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
                &                                      psv_i, psv_i )
+            write(numout,*)'ice_dyn_adv_umx 6:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
             !== Ice heat content ==!
             DO jk = 1, nlay_i
@@ -254,4 +270,5 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, psv_i, zua_ups, zva_ups )
+            write(numout,*)'ice_dyn_adv_umx 7:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
             !== Ice heat content ==!
             DO jk = 1, nlay_i
@@ -296,4 +313,6 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
                &                                      zhvar, psv_i, zuv_ups, zvv_ups )
+
+            write(numout,*)'ice_dyn_adv_umx 8:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
             !== Ice heat content ==!
             DO jk = 1, nlay_i
@@ -360,4 +379,6 @@
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
+
+         write(numout,*)'ice_dyn_adv_umx 9:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
          !
          ! Make sure ice thickness is not too big
@@ -365,5 +386,9 @@
          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
 
-      END DO
+         write(numout,*)'ice_dyn_adv_umx 10:psv_i, pe_s, jt = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1), ' ', jt
+
+      END DO
+
+      write(numout,*)'ice_dyn_adv_umx 11:psv_i, pe_s = ',psv_i(3,4,1), ' ', pe_i (3,4,1,1)
       !
    END SUBROUTINE ice_dyn_adv_umx
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg.F90	(revision 12475)
@@ -63,4 +63,5 @@
       ! controls
       IF( ln_timing    )   CALL timing_start('icedyn_rhg')                                                             ! timing
+      write(numout,*)'ice_dyn_rhg 1: u_ice = ',u_ice(3,4)
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
@@ -77,4 +78,5 @@
          END DO
       ENDIF
+      write(numout,*)'ice_dyn_rhg 2: u_ice = ',u_ice(3,4)
       !
       !--------------!
@@ -86,4 +88,5 @@
          !                             !------------------------!
          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
+         write(numout,*)'ice_dyn_rhg 3: u_ice = ',u_ice(3,4)
          !         
       END SELECT
@@ -91,10 +94,14 @@
       IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
          IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
+         write(numout,*)'ice_dyn_rhg 4: u_ice = ',u_ice(3,4)
       ENDIF
       !
       ! controls
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      write(numout,*)'ice_dyn_rhg 5: u_ice = ',u_ice(3,4)
       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints
-      IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')                                                             ! timing
+      write(numout,*)'ice_dyn_rhg 6: u_ice = ',u_ice(3,4)
+      IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')  
+      write(numout,*)'ice_dyn_rhg 7: u_ice = ',u_ice(3,4)                                                           ! timing
       !
    END SUBROUTINE ice_dyn_rhg
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg_evp.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg_evp.F90	(revision 12475)
@@ -185,4 +185,6 @@
 
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
+
+      write(numout,*)'ice_dyn_rhg_evp 1: u_ice = ',u_ice(3,4)
       !
 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
@@ -225,4 +227,6 @@
       CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
 
+      write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
+
       !------------------------------------------------------------------------------!
       ! 1) define some variables and initialize arrays
@@ -267,4 +271,6 @@
       ELSE                                               ;   zkt = 0._wp
       ENDIF
+
+      write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
       !
       !------------------------------------------------------------------------------!
@@ -328,4 +334,6 @@
       END DO
       CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
+
+      write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
       !
       !                                  !== Landfast ice parameterization ==!
@@ -368,4 +376,6 @@
       IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
 
+      write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
+
       !------------------------------------------------------------------------------!
       ! 3) Solution of the momentum equation, iterative procedure
@@ -397,4 +407,6 @@
          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
 
+         write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
+
          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
             DO ji = 2, jpi ! no vector loop
@@ -440,4 +452,6 @@
          END DO
          CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
+
+         write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
 
          DO jj = 1, jpjm1
@@ -488,4 +502,6 @@
             END DO
          END DO
+         write(numout,*)'ice_dyn_rhg_evp 3: u_ice = ',u_ice(3,4)
+
          !
          ! --- Computation of ice velocity --- !
@@ -558,4 +574,11 @@
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+
+                  IF (ji == 3 .AND. jj == 4) THEN
+                    write(numout,*)'ice_dyn_rhg_evp 6: u_ice(ji,jj), u_oce(ji,jj), v_iceU(ji,jj), v_oceU(ji,jj) = ',u_ice(ji,jj), u_oce(ji,jj), v_iceU(ji,jj), v_oceU(ji,jj)
+                    write(numout,*)'ice_dyn_rhg_evp 6: zTauO, zTauU_ib(ji,jj), zvel =',zTauO, zTauU_ib(ji,jj), zvel
+                    write(numout,*)'ice_dyn_rhg_evp 6: r1_e1u(ji,jj), zmf(ji,jj), e1v(ji,jj) = ',r1_e1u(ji,jj), zmf(ji,jj), e1v(ji,jj)
+                  ENDIF
+
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -580,4 +603,11 @@
                      &            ) * zmaskU(ji,jj)
                   ENDIF
+
+                  IF (ji == 3 .AND. jj == 4) THEN
+                    write(numout,*)'ice_dyn_rhg_evp 7: u_ice(ji,jj), zmU_t(ji,jj), zTauE, zTauO = ',u_ice(ji,jj), zmU_t(ji,jj), zTauE, zTauO
+                    write(numout,*)'ice_dyn_rhg_evp 7: zepsi, zTauO, zTauB, rswitch = ', zepsi, zTauO, zTauB, rswitch
+                    write(numout,*)'ice_dyn_rhg_evp 7: zdtevp, rn_lfrelax, zswitchU(ji,jj), zmaskU(ji,jj) = ', zdtevp, rn_lfrelax, zswitchU(ji,jj), zmaskU(ji,jj)
+                  ENDIF
+
                END DO
             END DO
@@ -608,4 +638,11 @@
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+
+                  IF (ji == 3 .AND. jj == 4) THEN
+                    write(numout,*)'ice_dyn_rhg_evp 8: u_ice(ji,jj), u_oce(ji,jj), v_iceU(ji,jj), v_oceU(ji,jj) = ',u_ice(ji,jj), u_oce(ji,jj), v_iceU(ji,jj), v_oceU(ji,jj)
+                    write(numout,*)'ice_dyn_rhg_evp 8: zTauO, zTauU_ib(ji,jj), zvel =',zTauO, zTauU_ib(ji,jj), zvel
+                    write(numout,*)'ice_dyn_rhg_evp 8: r1_e1u(ji,jj), zmf(ji,jj), e1v(ji,jj) = ',r1_e1u(ji,jj), zmf(ji,jj), e1v(ji,jj)
+                  ENDIF
+
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
@@ -630,7 +667,15 @@
                      &            ) * zmaskU(ji,jj)
                   ENDIF
+
+                  IF (ji == 3 .AND. jj == 4) THEN
+                    write(numout,*)'ice_dyn_rhg_evp 9: u_ice(ji,jj), zmU_t(ji,jj), zTauE, zTauO = ',u_ice(ji,jj), zmU_t(ji,jj), zTauE, zTauO
+                    write(numout,*)'ice_dyn_rhg_evp 9: zepsi, zTauO, zTauB, rswitch = ', zepsi, zTauO, zTauB, rswitch
+                    write(numout,*)'ice_dyn_rhg_evp 9: zdtevp, rn_lfrelax, zswitchU(ji,jj), zmaskU(ji,jj) = ', zdtevp, rn_lfrelax, zswitchU(ji,jj), zmaskU(ji,jj)
+                  ENDIF
                END DO
             END DO
             CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
+
+            write(numout,*)'ice_dyn_rhg_evp 10: u_ice = ',u_ice(3,4)
             !
 #if defined key_agrif
@@ -706,4 +751,5 @@
          CALL bdy_ice_dyn( 'V' )
       ENDIF
+      write(numout,*)'ice_dyn_rhg_evp 5: u_ice = ',u_ice(3,4)
       !
       !------------------------------------------------------------------------------!
@@ -757,4 +803,5 @@
       pstress2_i (:,:) = zs2 (:,:)
       pstress12_i(:,:) = zs12(:,:)
+      write(numout,*)'ice_dyn_rhg_evp 4: u_ice = ',u_ice(3,4)
       !
 
@@ -853,4 +900,6 @@
             END DO
          END DO
+
+         write(numout,*)'ice_dyn_rhg_evp 11: u_ice = ',u_ice(3,4)
          
          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
@@ -887,4 +936,6 @@
 
       ENDIF
+
+      write(numout,*)'ice_dyn_rhg_evp 2: u_ice = ',u_ice(3,4)
       !
    END SUBROUTINE ice_dyn_rhg_evp
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90	(revision 12475)
@@ -188,5 +188,5 @@
          kt_ice = kt                              ! -- Ice model time step
 
-         write(numout,*)'icestp 0: e_s = ',e_s (42,26,1,1)
+         write(numout,*)'icestp 0: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
          !
          u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
@@ -215,5 +215,5 @@
                                         CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )      
 
-         write(numout,*)'icestp 1: e_s = ',e_s (42,26,1,1)
+         write(numout,*)'icestp 1: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
     
          !-------------------------------------!
@@ -226,23 +226,23 @@
             &                           CALL ice_dyn( kt )            ! -- Ice dynamics
 
-         write(numout,*)'icestp 2: e_s = ',e_s (42,26,1,1)
+         write(numout,*)'icestp 2: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
          !
          !                          !==  lateral boundary conditions  ==!
          IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
 
-         write(numout,*)'icestp 3: e_s = ',e_s (42,26,1,1)
+         write(numout,*)'icestp 3: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
          !
          !                          !==  previous lead fraction and ice volume for flux calculations
                                         CALL ice_var_glo2eqv          ! h_i and h_s for ice albedo calculation
 
-         write(numout,*)'icestp 4: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
+         write(numout,*)'icestp 4: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
 
                                         CALL ice_var_agg(1)           ! at_i for coupling 
 
-         write(numout,*)'icestp 5: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
+         write(numout,*)'icestp 5: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
 
                                         CALL store_fields             ! Store now ice values
 
-         write(numout,*)'icestp 6: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
+         write(numout,*)'icestp 6: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
          !
          !------------------------------------------------------!
@@ -260,13 +260,13 @@
          !------------------------------------------------------!
                                         CALL ice_sbc_flx( kt, ksbc )
-         write(numout,*)'icestp 7: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
+         write(numout,*)'icestp 7: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)
          !----------------------------!
          ! --- ice thermodynamics --- !
          !----------------------------!
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics   
-         write(numout,*)'icestp 8: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)   
+         write(numout,*)'icestp 8: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)   
          !
                                         CALL ice_cor( kt , 2 )        ! -- Corrections
-         write(numout,*)'icestp 9: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
+         write(numout,*)'icestp 9: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)
          !
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
@@ -490,4 +490,5 @@
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
       wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
+      wfx_snp(:,:) = 0._wp
       wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
       wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90	(revision 12475)
@@ -442,4 +442,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snp_1d (1:npti), wfx_snp          )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res          )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
@@ -536,4 +537,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snp_1d (1:npti), wfx_snp        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90	(revision 12475)
@@ -296,6 +296,5 @@
                ! Snow melt trapped in ponds contrbutes to ice thickness and snow to ice conversion
                h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah_pnd(ji,jk) * rhos * r1_rhoi ) 
-               wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice
-               wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice
+               wfx_snp_1d(ji)     = wfx_snp_1d(ji)     + a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice  ! As zdeltah_pnd is negative this comes out as negative
                !
             ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90	(revision 12475)
@@ -102,5 +102,4 @@
       !
       REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only) 
-      LOGICAL,  DIMENSION(jpij) ::   l_snow       ! is snow above the threshold
       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
@@ -157,12 +156,9 @@
          ! If the snow thickness drops below zhs_min then reduce the snow fraction instead 
          IF( h_s_1d(ji) < zhs_min ) THEN 
-           ! isnow(ji) = h_s_1d(ji) / zhs_min 
-           isnow(ji) = 0.0_wp
+           isnow(ji) = h_s_1d(ji) / zhs_min 
            zh_s(ji) = zhs_min * r1_nlay_s 
-           l_snow(ji) = .false.
          ELSE 
            isnow(ji) = 1.0_wp 
            zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
-           l_snow(ji) = .true.
          END IF 
 
@@ -322,5 +318,4 @@
                IF ( .NOT. l_T_converged(ji) ) &
                   zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
-               IF ( .NOT. l_snow(ji) ) zkappa_s(ji,jk) = 0.0_wp
             END DO
          END DO
@@ -370,9 +365,5 @@
          DO jk = 1, nlay_s
             DO ji = 1, npti
-               IF (l_snow(ji)) then
-                  zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
-               ELSE
-                  zeta_s(ji,jk) = 0.0_wp
-               END IF
+               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
             END DO
          END DO
@@ -722,5 +713,5 @@
             DO ji = 1, npti
                !                               !---------------------!
-               IF( l_snow(ji) ) THEN   !  snow-covered cells !
+               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
                   !                            !---------------------!
                   ! snow interior terms (bottom equation has the same form as the others)
@@ -779,7 +770,12 @@
                      zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
                         &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
-                  ENDIF
+                  ENDIF               
                   
                ENDIF
+
+                  IF (to_print(ji) == 10) THEN
+                     write(numout,*)'after minmax: h_s_1d(ji), jm_min(ji), nlay_s = ',h_s_1d(ji), ' ', jm_min(ji), ' ', nlay_s
+                  ENDIF   
+
                !
                zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
@@ -802,5 +798,12 @@
             DO jk = jm_mint+1, jm_maxt
                DO ji = 1, npti
+
+ 
                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
+                  IF (to_print(ji) == 10) THEN
+                     write(numout,*)'before calc zdiagbis: jk, jm, zdiagbis(ji,jm-1) = ',jk, ' ', jm, ' ', zdiagbis(ji,jm-1)
+                     write(numout,*)'before calc zdiagbis:ztrid(ji,jm,2), ztrid(ji,jm,1), ztrid(ji,jm-1,3) = ', ztrid(ji,jm,2), ' ', ztrid(ji,jm,1), ' ', ztrid(ji,jm-1,3)
+                  ENDIF
+
                   zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
                   zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
@@ -838,4 +841,5 @@
                 write(numout,*)'icethd_zdf_bl99 2: t_s_1d(ji,1), e_s_1d(ji,1) = ',t_s_1d(ji,1), ' ', e_s_1d(ji,1)
                 write(numout,*)'icethd_zdf_bl99 2: zindtbis(ji,nlay_s+1), ztrid(ji,nlay_s+1,3), t_i_1d(ji,1) = ', zindtbis(ji,nlay_s+1), ' ', ztrid(ji,nlay_s+1,3), ' ', t_i_1d(ji,1)
+                write(numout,*)'icethd_zdf_bl99 2: zdiagbis(ji,nlay_s+1), h_s_1d(ji) = ',zdiagbis(ji,nlay_s+1), ' ', h_s_1d(ji)
               END IF
             END DO
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90	(revision 12475)
@@ -154,4 +154,6 @@
                &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
 
+
+
             ! add the snow melt water to snow mass flux to the ocean
             wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
@@ -161,4 +163,6 @@
             emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
 
+
+            
 
             ! Salt flux at the ocean surface      
@@ -166,4 +170,6 @@
             sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
                &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) + sfx_pnd(ji,jj)
+
+
             
             ! Mass of snow and ice per unit area   
@@ -177,4 +183,31 @@
          END DO
       END DO
+
+            write(numout,*)'iceupdate: MEAN(wfx_pnd) = ',SUM(wfx_pnd)/SIZE(wfx_pnd)
+            write(numout,*)'iceupdate: MEAN(wfx_sum) = ',SUM(wfx_sum)/SIZE(wfx_sum)
+            write(numout,*)'iceupdate: MEAN(wfx_bom) = ',SUM(wfx_bom)/SIZE(wfx_bom)
+            write(numout,*)'iceupdate: MEAN(wfx_ice) = ',SUM(wfx_ice)/SIZE(wfx_ice)
+            write(numout,*)'iceupdate: MEAN(wfx_sni) = ',SUM(wfx_sni)/SIZE(wfx_sni)
+            write(numout,*)'iceupdate: MEAN(wfx_snp) = ',SUM(wfx_snp)/SIZE(wfx_snp)
+            write(numout,*)'iceupdate: MEAN(fmmflx) = ',SUM(fmmflx)/SIZE(fmmflx)
+            write(numout,*)'iceupdate: MEAN(wfx_snw) = ',SUM(wfx_snw)/SIZE(wfx_snw)
+            write(numout,*)'iceupdate: MEAN(emp) = ',SUM(emp)/SIZE(emp)
+            write(numout,*)'iceupdate: MEAN(emp_oce) = ',SUM(emp_oce)/SIZE(emp_oce)
+            write(numout,*)'iceupdate: MEAN(sfx) = ',SUM(sfx)/SIZE(sfx)
+            write(numout,*)'iceupdate: MEAN(sfx_sni) = ',SUM(sfx_sni)/SIZE(sfx_sni)
+            write(numout,*)'iceupdate: MEAN(sfx_pnd) = ',SUM(sfx_pnd)/SIZE(sfx_pnd)
+            write(numout,*)'iceupdate: MEAN(sfx_sum) = ',SUM(sfx_sum)/SIZE(sfx_sum)
+            write(numout,*)'iceupdate: MEAN(sfx_bom) = ',SUM(sfx_sum)/SIZE(sfx_bom)
+            write(numout,*)'iceupdate: MAXVAL(wfx_pnd) = ',MAXVAL(wfx_pnd)
+            write(numout,*)'iceupdate: MINVAL(wfx_pnd) = ',MINVAL(wfx_pnd)
+            write(numout,*)'iceupdate: MAXVAL(wfx_bom) = ',MAXVAL(wfx_bom)
+            write(numout,*)'iceupdate: MINVAL(wfx_bom) = ',MINVAL(wfx_bom)
+            write(numout,*)'iceupdate: MAXVAL(wfx_sni) = ',MAXVAL(wfx_sni)
+            write(numout,*)'iceupdate: MINVAL(wfx_sni) = ',MINVAL(wfx_sni)
+            write(numout,*)'iceupdate: MAXVAL(sfx_pnd) = ',MAXVAL(sfx_pnd)
+            write(numout,*)'iceupdate: MINVAL(sfx_pnd) = ',MINVAL(sfx_pnd)
+            write(numout,*)'iceupdate: MAXVAL(sfx_bom) = ',MAXVAL(sfx_bom)
+            write(numout,*)'iceupdate: MINVAL(sfx_bom) = ',MINVAL(sfx_bom)
+
 
       ! Storing the transmitted variables
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icewri.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icewri.F90	(revision 12475)
@@ -172,9 +172,10 @@
       !------------------
       ! trends
-      IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res - wfx_pnd ) ! Sea-ice mass change from thermodynamics
+      IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res - wfx_pnd - wfx_snp ) ! Sea-ice mass change from thermodynamics
       IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s)
       IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water
       IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth
       IF( iom_use('dmisni') )   CALL iom_put( "dmisni", - wfx_sni                           )   ! Sea-ice mass change through snow-to-ice conversion
+      IF( iom_use('dmisnp') )   CALL iom_put( "dmisnp", - wfx_snp                           )   ! Sea-ice mass change through snow melt into ponds
       IF( iom_use('dmisum') )   CALL iom_put( "dmisum", - wfx_sum                           )   ! Sea-ice mass change through surface melting
       IF( iom_use('dmipnd') )   CALL iom_put( "dmipnd", - wfx_pnd                           )   ! Sea-ice mass change through meltponds overflowing to ocean
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/IOM/restart.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/IOM/restart.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/IOM/restart.F90	(revision 12475)
@@ -281,4 +281,7 @@
          neuler = 0
       ENDIF
+
+      write(numout,*)'rst_read 2: ub(3,4,1) = ', ub(3,4,1)
+
       !
       CALL iom_get( numror, jpdom_autoglo, 'un'     , un, ldxios = lrxios )   ! now    fields
@@ -306,4 +309,5 @@
          !
       ENDIF
+      write(numout,*)'rst_read 3: ub(3,4,1) = ', ub(3,4,1)
       !
    END SUBROUTINE rst_read
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcmod.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcmod.F90	(revision 12475)
@@ -407,4 +407,7 @@
       !                                            !        forcing field computation         !
       !                                            ! ---------------------------------------- !
+
+      write(numout,*)'sbcmod 1: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       !
       ll_sas = nn_components == jp_iam_sas               ! component flags
@@ -413,4 +416,6 @@
       IF( .NOT.ll_sas )   CALL sbc_ssm ( kt )            ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
       IF( ln_wave     )   CALL sbc_wave( kt )            ! surface waves
+
+      write(numout,*)'sbcmod 2: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
 
       !
@@ -436,4 +441,6 @@
       !                                            !==  Misc. Options  ==!
       !
+      write(numout,*)'sbcmod 3: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
       CASE(  1 )   ;         CALL sbc_ice_if   ( kt )             ! Ice-cover climatology ("Ice-if" model)
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcssm.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcssm.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbcssm.F90	(revision 12475)
@@ -59,4 +59,5 @@
       !!---------------------------------------------------------------------
       !
+      write(numout,*)'sbcssm 1: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
       DO jj = 1, jpj
@@ -83,4 +84,5 @@
          !
          frq_m(:,:) = fraqsr_1lev(:,:)
+         write(numout,*)'sbcssm 2: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
          !
       ELSE
@@ -106,4 +108,5 @@
             !
             frq_m(:,:) = zcoef * fraqsr_1lev(:,:)
+            write(numout,*)'sbcssm 3: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
             !                                             ! ---------------------------------------- !
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
@@ -116,4 +119,5 @@
             e3t_m(:,:) = 0._wp
             frq_m(:,:) = 0._wp
+            write(numout,*)'sbcssm 4: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
          ENDIF
          !                                                ! ---------------------------------------- !
@@ -121,4 +125,5 @@
          !                                                ! ---------------------------------------- !
          ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1)
+         write(numout,*)'sbcssm 5: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,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) )
@@ -146,4 +151,5 @@
             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,*)'sbcssm 6: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
             !
          ENDIF
@@ -167,4 +173,5 @@
             !
             IF( lwxios ) CALL iom_swap(      cxios_context          )
+            write(numout,*)'sbcssm 7: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
          ENDIF
          !
@@ -180,4 +187,7 @@
          CALL iom_put( 'frq_m', frq_m )
       ENDIF
+
+      write(numout,*)'sbcssm 8: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       !
    END SUBROUTINE sbc_ssm
@@ -195,4 +205,7 @@
       !!----------------------------------------------------------------------
       !
+
+      write(numout,*)'sbcssm 1: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       IF( nn_fsbc == 1 ) THEN
          !
@@ -216,4 +229,7 @@
             CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m, ldxios = lrxios  )    !   "         "    height      (T-point)
             CALL iom_get( numror, jpdom_autoglo, 'e3t_m'  , e3t_m, ldxios = lrxios  )    ! 1st level thickness          (T-point)
+
+            write(numout,*)'sbcssm 2: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
             ! fraction of solar net radiation absorbed in 1st T level
             IF( iom_varid( numror, 'frq_m', ldstop = .FALSE. ) > 0 ) THEN
@@ -233,4 +249,5 @@
                e3t_m(:,:) = zcoef * e3t_m(:,:)
                frq_m(:,:) = zcoef * frq_m(:,:)
+               write(numout,*)'sbcssm 3: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
             ELSE
                IF(lwp) WRITE(numout,*) '   mean fields read in the ocean restart file'
@@ -238,4 +255,7 @@
          ENDIF
       ENDIF
+
+      write(numout,*)'sbcssm 4: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       !
       IF( .NOT.l_ssm_mean ) THEN   ! default initialisation. needed by iceistate
@@ -253,4 +273,7 @@
          !
       ENDIF
+
+      write(numout,*)'sbcssm 5: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
       !
       IF( .NOT. ln_traqsr )   fraqsr_1lev(:,:) = 1._wp   ! default definition: qsr 100% in the fisrt level 
@@ -267,4 +290,6 @@
       ENDIF
 
+      write(numout,*)'sbcssm 6: ssu_m(3,4), ub(3,4,1) = ', ssu_m(3,4), ' ', ub(3,4,1)
+
    END SUBROUTINE sbc_ssm_init
 
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/step.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/step.F90	(revision 12474)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/step.F90	(revision 12475)
@@ -93,4 +93,5 @@
       !
       IF( ln_timing )   CALL timing_start('stp')
+      write(numout,*)'stp 1: ub(3,4,1) = ', ub(3,4,1)
       !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -110,7 +111,13 @@
       ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      write(numout,*)'stp 2: ub(3,4,1) = ', ub(3,4,1)
       IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
+      write(numout,*)'stp 3: ub(3,4,1) = ', ub(3,4,1)
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
+      write(numout,*)'stp 4: ub(3,4,1) = ', ub(3,4,1)
       IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+
+      write(numout,*)'stp 5: ub(3,4,1) = ', ub(3,4,1)
+
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
 
