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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90	(revision 12369)
@@ -171,4 +171,5 @@
    INTEGER, PUBLIC, PARAMETER ::   np_cnd_ON  = 1  !: forcing from conduction flux (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
    INTEGER, PUBLIC, PARAMETER ::   np_cnd_EMU = 2  !: emulate conduction flux via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
+   INTEGER , PUBLIC :: cat
 
    !                                     !!** ice-vertical diffusion namelist (namthd_zdf) **
@@ -296,4 +297,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_i       !: Ice thickness                           (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i       !: Ice fractional areas (concentration)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple    !: Ice fractional area at last coupling time
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_i       !: Ice volume per unit area                (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_s       !: Snow volume per unit area               (m)
@@ -379,4 +381,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top   !: Surface conduction flux (W/m2)
 
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   to_print_2d
    !
    !!----------------------------------------------------------------------
@@ -393,5 +396,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(15), ii
+      INTEGER :: ierr(16), ii
       !!-----------------------------------------------------------------
       ierr(:) = 0
@@ -428,4 +431,7 @@
 
       ii = ii + 1
+      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) )
+
+      ii = ii + 1
       ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                   &
          &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  &
@@ -469,4 +475,6 @@
       ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
 
+      ALLOCATE( to_print_2d(jpi,jpj), STAT = ierr(ii) )
+
       ice_alloc = MAXVAL( ierr(:) )
       IF( ice_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90	(revision 12369)
@@ -35,5 +35,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d     
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d  
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   to_print
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d
    
@@ -186,5 +187,5 @@
          &      hfx_sum_1d(jpij) , hfx_bom_1d    (jpij) , hfx_bog_1d(jpij) ,   & 
          &      hfx_dif_1d(jpij) , hfx_opw_1d    (jpij) , hfx_dyn_1d(jpij) ,   &
-         &      rn_amax_1d(jpij) ,                                             &
+         &      rn_amax_1d(jpij) , to_print(jpij) ,                            &
          &      hfx_thd_1d(jpij) , hfx_spr_1d    (jpij) ,                      &
          &      hfx_snw_1d(jpij) , hfx_sub_1d    (jpij) ,                      &
Index: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icecor.F90	(revision 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icecor.F90	(revision 12369)
@@ -85,4 +85,7 @@
       IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
       !                             !-----------------------------------------------------
+
+      WRITE(numout,*) 'icecor: sv_i before bounding = ',sv_i(3,4,1)
+
          zzc = rhoi * r1_rdtice
          DO jl = 1, jpl
@@ -92,12 +95,22 @@
                   sv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
+
+                  IF ( ji == 3 .AND. jj == 4 .AND. jl == 1 ) THEN
+                     WRITE(numout,*) 'icecor: zsal, v_i(ji,jj,jl), sv_i(ji,jj,jl) = ',zsal, v_i(ji,jj,jl), sv_i(ji,jj,jl)
+                  ENDIF
                END DO
             END DO
          END DO
       ENDIF
+
+
+      WRITE(numout,*) 'icecor before redist: v_i(42,26,1), a_i(42,26,1), h_i = ',v_i(42,26,1), ' ', a_i(42,26,1), ' ', v_i(42,26,1)/a_i(42,26,1)
+
       !                             !-----------------------------------------------------
       !                             !  Rebin categories with thickness out of bounds     !
       !                             !-----------------------------------------------------
       IF ( jpl > 1 )   CALL ice_itd_reb( kt )
+
+      WRITE(numout,*) 'icecor after redist: v_i(42,26,1), a_i(42,26,1), h_i = ',v_i(42,26,1), ' ', a_i(42,26,1), ' ', v_i(42,26,1)/a_i(42,26,1)
 
       !                             !-----------------------------------------------------
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn.F90	(revision 12369)
@@ -103,4 +103,6 @@
       END WHERE
       !
+
+      write(numout,*)'ice_dyn 0: e_i = ',e_i (42,26,1,1)
       !
       SELECT CASE( nice_dyn )          !-- Set which dynamics is running
@@ -146,4 +148,6 @@
 
       END SELECT
+
+      write(numout,*)'ice_dyn 1: e_i = ',e_i (42,26,1,1)
       !
       !
@@ -169,4 +173,5 @@
          !
       ENDIF
+      write(numout,*)'ice_dyn 2: e_i = ',e_i (42,26,1,1)
       !
       ! controls
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90	(revision 12369)
@@ -80,4 +80,6 @@
    USE prtctl         ! Print control
 
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+
    IMPLICIT NONE
    PRIVATE
@@ -118,8 +120,65 @@
       INTEGER, INTENT(in) ::   ksbc    ! flux formulation (user defined, bulk, or Pure Coupled)
       !
-      INTEGER ::   jl   ! dummy loop index
+      INTEGER  :: ji, jj, jk, jl, jn, gi, gj   ! dummy loop indices
+      REAL(wp) :: delta_lon, delta_lat, lon_west, lon_east, lat_south, lat_north, lon_print, lat_print
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('ice_stp')
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            gi = nimpp + ji - 1
+            gj = njmpp + jj - 1
+            delta_lon = (glamt(2,jj) - glamt(1,jj))/2.0
+            delta_lat = (gphit(ji,2) - gphit(ji,1))/2.0
+
+            lon_west = glamt(ji,jj) - delta_lon
+            lon_east = glamt(ji,jj) + delta_lon
+            lat_south = gphit(ji,jj) - delta_lat
+            lat_north = gphit(ji,jj) + delta_lat
+            to_print_2d(ji,jj) = 0
+
+            lon_print = 95.0
+            lat_print = -63.0  ! For bottom melt investigation
+
+            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              to_print_2d(ji,jj) = 1
+              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
+            ENDIF
+
+            lon_print = 10.0
+            lat_print = -62.0   ! For top melt investigation
+
+            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              to_print_2d(ji,jj) = 2
+              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
+            ENDIF
+
+            lon_print = 137.81
+            lat_print = -66   ! For conductivity investigation
+
+            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              to_print_2d(ji,jj) = 3
+              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
+            ENDIF
+
+            lon_print = 137.81 - 360.0
+
+            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              to_print_2d(ji,jj) = 3
+              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
+            ENDIF
+
+            lon_print = 8
+            lat_print = -67   ! For top melt investigation
+
+            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              to_print_2d(ji,jj) = 10
+              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
+            ENDIF
+          END DO
+        END DO
+
       !
       !                                      !-----------------------!
@@ -128,4 +187,6 @@
          !
          kt_ice = kt                              ! -- Ice model time step
+
+         write(numout,*)'icestp 0: e_s = ',e_s (42,26,1,1)
          !
          u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
@@ -152,5 +213,8 @@
          !    utau_ice, vtau_ice = surface ice stress [N/m2]
          !------------------------------------------------!
-                                        CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )          
+                                        CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )      
+
+         write(numout,*)'icestp 1: e_s = ',e_s (42,26,1,1)
+    
          !-------------------------------------!
          ! --- ice dynamics and advection  --- !
@@ -161,12 +225,24 @@
          IF( ln_icedyn .AND. .NOT.lk_c1d )   &
             &                           CALL ice_dyn( kt )            ! -- Ice dynamics
+
+         write(numout,*)'icestp 2: e_s = ',e_s (42,26,1,1)
          !
          !                          !==  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)
          !
          !                          !==  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)
+
                                         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)
+
                                         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)
          !
          !------------------------------------------------------!
@@ -184,15 +260,21 @@
          !------------------------------------------------------!
                                         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)
          !----------------------------!
          ! --- ice thermodynamics --- !
          !----------------------------!
-         IF( ln_icethd )                CALL ice_thd( kt )            ! -- 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)   
          !
                                         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)
          !
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
+         write(numout,*)'icestp 10: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
+         write(numout,*)'icestp 11: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
          !
                                         CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
+         write(numout,*)'icestp 12: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
          !
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
@@ -203,4 +285,5 @@
          !
          IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
+         write(numout,*)'icestp 13: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
          !
       ENDIF   ! End sea-ice time step only
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90	(revision 12369)
@@ -85,9 +85,11 @@
       INTEGER, INTENT(in) :: kt    ! number of iteration
       !
-      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
+      INTEGER  :: ji, jj, jk, jl, jn, gi, gj   ! dummy loop indices
       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
+
+      REAL(wp) :: delta_lon, delta_lat, lon_west, lon_east, lat_south, lat_north, lon_print, lat_print
       !
       !!-------------------------------------------------------------------
@@ -152,8 +154,17 @@
             qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
 
+            If (to_print_2d(ji,jj) == 1) THEN
+              write(numout,*) 'icethd: qsb_ice_bot(ji,jj), rswitch, rau0, rcp = ',qsb_ice_bot(ji,jj), rswitch, rau0, rcp
+              write(numout,*) 'icethd: zch, zfric_u, sst_m(ji,jj), t_bo(ji,jj) = ',zch, zfric_u, sst_m(ji,jj), t_bo(ji,jj)
+            ENDIF
+
             qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
             ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
             !                              the freezing point, so that we do not have SST < T_freeze
             !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
+
+            If (to_print_2d(ji,jj) == 1) THEN
+              write(numout,*) 'icethd: qsb_ice_bot(ji,jj) after adjust = ',qsb_ice_bot(ji,jj)
+            ENDIF
 
             !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice
@@ -168,4 +179,11 @@
                fhld (ji,jj) = 0._wp
             ENDIF
+
+            If (to_print_2d(ji,jj) == 1) THEN
+              write(numout,*) 'icethd: fhld(ji,jj), rswitch, zqld, at_i(ji,jj) = ', fhld(ji,jj), rswitch, zqld, at_i(ji,jj)
+              write(numout,*) 'icethd: qsr_oce(ji,jj), frq_m(ji,jj), qns_oce(ji,jj), qemp_oce(ji,jj) = ',qsr_oce(ji,jj), frq_m(ji,jj), qns_oce(ji,jj), ' ', qemp_oce(ji,jj)
+              write(numout,*) 'icethd: at_i_b(ji,jj), 1-at_i_b(ji,jj) = ',at_i_b(ji,jj), 1-at_i_b(ji,jj)
+            ENDIF
+
             !
             ! Net heat flux on top of the ice-ocean [W.m-2]
@@ -197,5 +215,11 @@
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
       !-------------------------------------------------------------------------------------------!
+
+      write(numout,*)'icethd: sz_i(3,4,1,1) = ',sz_i(3,4,1,1)
+      write(numout,*)'icethd: t_i(3,4,1,1) = ',t_i(3,4,1,1)
+
       DO jl = 1, jpl
+
+         cat = jl
 
          ! select ice covered grid points
@@ -215,4 +239,36 @@
             !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
             !
+
+            DO jn = 1, npti
+              ji        = MOD( nptidx(jn) - 1 , jpi ) + 1
+              jj        =    ( nptidx(jn) - 1 ) / jpi + 1
+              
+              !IF ( jl == 1 ) THEN
+                to_print(jn) = to_print_2d(ji,jj)
+              !ELSE
+              !  to_print(jn) = 0
+              !END IF
+
+              !gi = nimpp + ji - 1
+              !gj = njmpp + jj - 1
+              !to_print(jn) = 0
+
+              !delta_lon = (glamt(2,jj) - glamt(1,jj))/2.0
+              !delta_lat = (gphit(ji,2) - gphit(ji,1))/2.0
+
+              !lon_west = glamt(ji,jj) - delta_lon
+              !lon_east = glamt(ji,jj) + delta_lon
+              !lat_south = gphit(ji,jj) - delta_lat
+              !lat_north = gphit(ji,jj) + delta_lat
+              !lon_print = 95.0
+              !lat_print = -63.0
+
+              !IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
+              !  to_print(jn) = 1
+              !  write(numout,*)'icethd: jn, ji, jj, gi,gj = ', jn, ji, jj, gi,gj
+              !ENDIF
+              
+            END DO
+
             s_i_new   (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp  ! --- some init --- !  (important to have them here) 
             dh_i_sum  (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm  (1:npti) = 0._wp 
@@ -230,4 +286,8 @@
             !
                               CALL ice_thd_temp                     ! --- Temperature update --- !
+
+            
+
+
             !
             IF( ln_icedH .AND. ln_virtual_itd ) &
Index: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_da.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_da.F90	(revision 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_da.F90	(revision 12369)
@@ -145,4 +145,15 @@
             ! Contribution to mass flux
             wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_rdtice * ( rhoi * h_i_1d(ji) + rhos * h_s_1d(ji) )
+
+            IF (to_print(ji) == 1) THEN
+              write(numout,*)'icethd_da: sst_1d(ji) = ',sst_1d(ji)
+            ENDIF
+
+            IF (to_print(ji) == 1) THEN
+              write(numout,*)'icethd_da: lateral melt: a_i_1d(ji), at_i_1d(ji), sst_1d(ji), t_bo_1d(ji) = ',a_i_1d(ji), at_i_1d(ji), sst_1d(ji), t_bo_1d(ji)
+              write(numout,*)'icethd_da: lateral melt: wfx_lam_1d(ji), zda, r1_rdtice, rhoi =', wfx_lam_1d(ji), zda, r1_rdtice, rhoi
+              write(numout,*)'icethd_da: lateral melt: h_i_1d(ji), rhos, h_s_1d(ji) = ',h_i_1d(ji), rhos, h_s_1d(ji)
+              write(numout,*)'icethd_da: lateral melt: zdfloe, zperi, zwlat, zda_tot(ji) = ',zdfloe, zperi, zwlat, zda_tot(ji)
+            ENDIF
             
             ! new concentration
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90	(revision 12369)
@@ -131,4 +131,8 @@
          DO ji = 1, npti
             zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+
+            IF ( to_print(ji) == 10 ) THEN
+               write(numout,*)'loading flux: zq_top(ji), qml_ice_1d(ji) = ',zq_top(ji), ' ', qml_ice_1d(ji)
+            ENDIF
          END DO
          !
@@ -139,4 +143,10 @@
             qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
             zq_top(ji)     = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+
+            IF ( to_print(ji) == 10 ) THEN
+               write(numout,*)'icethd_dh: qns_ice_1d(ji), qsr_ice_1d(ji), qtr_ice_top_1d(ji), qcn_ice_top_1d(ji) = ',qns_ice_1d(ji), ' ', qsr_ice_1d(ji), ' ', qtr_ice_top_1d(ji), ' ', qcn_ice_top_1d(ji)
+               write(numout,*)'icethd_dh: zdum, t_su_1d(ji), qml_ice_1d(ji), zq_top(ji) = ',zdum, ' ', t_su_1d(ji), qml_ice_1d(ji), zq_top(ji)
+            ENDIF
+
          END DO
          !
@@ -146,4 +156,8 @@
          zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
          zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rdt_ice )
+         IF ( to_print(ji) == 1 ) THEN
+            write(numout,*)'icethd_dh: qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji) = ',qcn_ice_bot_1d(ji), qsb_ice_bot_1d(ji)
+            write(numout,*)'icethd_dh: fhld_1d(ji), zq_bot(ji) = ',fhld_1d(ji), ', ', zq_bot(ji)
+         ENDIF
       END DO
 
@@ -181,6 +195,16 @@
                t_s_1d  (ji,jk) = rt0
             END IF
-         END DO
-      END DO         
+
+           IF (to_print(ji) == 10) THEN
+                write(numout,*)'icethd_dh 8: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
+            END IF
+         END DO
+      END DO 
+
+      DO ji = 1, npti
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_dh internal_melting: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
+         END IF
+      END DO        
 
       ! Snow precipitation
@@ -210,4 +234,8 @@
             zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )      
             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
+
+            IF ( to_print(ji) == 10 ) THEN
+               write(numout,*)'after precip: zq_top(ji), zdeltah(ji,1), zqprec(ji) = ',zq_top(ji), zdeltah(ji,1), zqprec(ji)
+            ENDIF
             
             ! update thickness
@@ -220,4 +248,9 @@
             !
          ENDIF
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_dh snow_precip: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
+         END IF
+
       END DO
 
@@ -253,5 +286,19 @@
                !
             ENDIF
-         END DO
+
+            IF ( to_print(ji) == 10 ) THEN
+               write(numout,*)'after snow melt: zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), test = ',zq_top(ji), zdeltah(ji,jk), e_s_1d(ji,jk), ( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp )
+            ENDIF
+
+            IF (to_print(ji) == 10) THEN
+               write(numout,*)'icethd_dh loop: wfx_snw_sum_1d(ji), jk = ',wfx_snw_sum_1d(ji), ' ', jk
+            END IF
+         END DO
+      END DO
+
+      DO ji = 1, npti
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_dh snow_ground: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
+         END IF
       END DO
 
@@ -301,4 +348,9 @@
               &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
+
+            IF (to_print(ji) == 10) THEN
+                write(numout,*)'icethd_dh 9: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
+                write(numout,*)'icethd_dh 9: zdh_s_pre(ji), zqprec(ji), h_s_1d(ji) = ',zdh_s_pre(ji), ' ', zqprec(ji), ' ', h_s_1d(ji)
+            END IF
          END DO
       END DO
@@ -362,4 +414,9 @@
                
             END IF
+
+            IF ( to_print(ji) == 10 ) THEN
+               write(numout,*)'after top melt: zq_top(ji), test = ',zq_top(ji), ( t_i_1d(ji,jk) >= (ztmelts+rt0) ) 
+               write(numout,*)'after top melt: wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE = ', wfx_sum_1d(ji), zdeltah(ji,jk), zfmdt, zdE
+            ENDIF
             
             ! Ice sublimation
@@ -510,4 +567,8 @@
                   zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change
 
+                  IF ( to_print(ji) == 1 ) THEN
+                    write(numout,*)'icethd_dh: zfmdt, zq_bot(ji), zdE, zdeltah(ji,jk) = ',zfmdt,', ', zq_bot(ji), zdE, zdeltah(ji,jk)
+                  ENDIF
+
                   zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
                   
@@ -527,4 +588,8 @@
                   
                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux
+
+                  IF (to_print(ji) == 1) THEN
+                     write(numout,*)'icethd_dh: wfx_bom_1d(ji), in meters, at_i_1d(ji) = ',wfx_bom_1d(ji),wfx_bom_1d(ji)*rdt_ice/(rhoi*at_i_1d(ji)), at_i_1d(ji)
+                  ENDIF
 
                   ! update heat content (J.m-2) and layer thickness
@@ -536,4 +601,6 @@
          END DO
       END DO
+
+      
 
       ! Update temperature, energy
@@ -564,4 +631,9 @@
 
          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_dh snow_extra: wfx_snw_sum_1d(ji) = ',wfx_snw_sum_1d(ji)
+         END IF
+
       END DO
 
@@ -610,4 +682,8 @@
          rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
          t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_dh: t_su_1d(ji), t_i_1d(ji,1) = ',t_su_1d(ji), ' ', t_i_1d(ji,1)
+         END IF
       END DO
 
@@ -625,4 +701,8 @@
             ! recalculate t_s_1d from e_s_1d
             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
+
+            IF (to_print(ji) == 10) THEN
+                write(numout,*)'icethd_dh 10: t_s_1d(ji,jk), e_s_1d(ji,jk) = ',t_s_1d(ji,jk), ' ', e_s_1d(ji,jk)
+            END IF
          END DO
       END DO
Index: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90	(revision 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90	(revision 12369)
@@ -144,4 +144,9 @@
 
       DO ji = 1, npti
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_pnd_start: h_ip_1d(ji), a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), ' ', a_ip_frac_1d(ji), ' ', a_ip_1d(ji)
+         END IF
+
          !                                                        !--------------------------------!
          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
@@ -186,4 +191,11 @@
             !
          ENDIF
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_pnd: h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji)
+            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt
+            write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
+         END IF
+
       END DO
       !
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 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90	(revision 12369)
@@ -143,4 +143,10 @@
       END DO
 
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 0: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
+
       !------------------
       ! 1) Initialization
@@ -151,4 +157,7 @@
          zh_i(ji) = h_i_1d(ji) * r1_nlay_i
          zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+         IF( to_print(ji) == 10 ) THEN
+           write(numout,*)'icethd_zdf_bl99: v_i_1d(ji), a_i_1d(ji), h_i_1d(ji) = ',v_i_1d(ji), ' ', a_i_1d(ji), ' ', h_i_1d(ji)
+         END IF
       END DO
       !
@@ -256,6 +265,20 @@
          DO ji = 1, npti
             IF ( .NOT. l_T_converged(ji) ) &
-               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )        
-         END DO
+               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )   
+
+            IF (to_print(ji) == 10) THEN
+               write(numout,*)'icethd_zdf_bl99: ln_cndi_U64, ln_cndi_P07, rcnd_i = ',ln_cndi_U64, ln_cndi_P07, rcnd_i
+               write(numout,*)'icethd_zdf_bl99: sz_i_1d(ji,1), epsi10, t_i_1d(ji,1) = ',sz_i_1d(ji,1), epsi10, t_i_1d(ji,1)
+               write(numout,*)'icethd_zdf_bl99: rt0, nlay_i, sz_i_1d(ji,nlay_i), t_bo_1d(ji) = ', rt0, nlay_i, sz_i_1d(ji,nlay_i), t_bo_1d(ji)
+               write(numout,*)'icethd_zdf_bl99: sz_i_1d(ji,:) = ',sz_i_1d(ji,:)
+               write(numout,*)'icethd_zdf_bl99: t_i_1d(ji,:) = ',t_i_1d(ji,:)
+               write(numout,*)'icethd_zdf_bl99: ztcond_i_cp(ji,:) = ',ztcond_i_cp(ji,:)
+               write(numout,*)'icethd_zdf_bl99: zkimin = ',zkimin
+               write(numout,*)'icethd_zdf_bl99: ztcond_i = ',ztcond_i(ji,:)
+            ENDIF
+         END DO
+
+         
+
          !
          !--- G(he) : enhancement of thermal conductivity in mono-category case
@@ -312,4 +335,8 @@
             IF ( .NOT. l_T_converged(ji) ) &
                zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+            
+            IF (to_print(ji) == 10) THEN
+               write(numout,*) 'icethd_zdf_bl99: zkappa_i(ji,0), ztcond_i(ji,0), z1_h_i(ji), isnow(ji) = ',zkappa_i(ji,0), ' ', ztcond_i(ji,0), ' ', z1_h_i(ji), ' ', isnow(ji)
+            END IF
          END DO
          !
@@ -321,4 +348,10 @@
                zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
                zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+
+               IF (to_print(ji) == 10 .AND. jk == 1 .AND. cat == 1) THEN
+                  write(numout,*) 'zeta_i(ji,jk), z1_h_i(ji), zcpi = ',zeta_i(ji,jk), ' ', z1_h_i(ji), ' ', zcpi
+                  write(numout,*) 'sz_i_1d(ji,jk), t_i_1d(ji,jk), ztiold(ji,jk) = ',sz_i_1d(ji,jk), ' ', t_i_1d(ji,jk), ' ', ztiold(ji,jk)
+                  write(numout,*) 'iconv, l_T_converged(ji) = ',iconv, ' ',l_T_converged(ji)
+               END IF
             END DO
          END DO
@@ -329,4 +362,5 @@
             END DO
          END DO
+
          !
          !----------------------------------------!
@@ -355,4 +389,10 @@
                zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
             END DO
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 16 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
             !
             !----------------------------
@@ -389,4 +429,10 @@
                   &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
             END DO
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 17.5 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
 
             DO ji = 1, npti
@@ -496,4 +542,10 @@
                !
             END DO
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'tri diagonal matrix terms zetai = ',zeta_i(ji,:)
+           END IF
+         END DO
             !
             !------------------------------
@@ -548,5 +600,16 @@
                   ENDIF
                ENDIF
-            END DO
+
+            IF (to_print(ji) == 10 ) THEN
+                write(numout,*)'icethd_zdf_bl99 1: 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 1: 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)
+            END IF
+
+            END DO
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 18 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
             !
             !--------------------------------------------------------------
@@ -579,4 +642,10 @@
             END DO
 
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 19 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
+
          !----------------------------------------!
          !                                        !
@@ -588,4 +657,9 @@
             !
             ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 7.1 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
             !
             !----------------------------
@@ -622,4 +696,12 @@
                   &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
             ENDDO
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10 .AND. iconv == 1) THEN
+             jm =  nlay_s + nlay_i + 1
+             write(numout,*) 'icethd_zdf_bl99: zindterm(ji,jm), ztiold(ji,nlay_i), zeta_i(ji,nlay_i), zradab_i(ji,nlay_i) = ',zindterm(ji,jm), ' ', ztiold(ji,nlay_i), ' ', zeta_i(ji,nlay_i), ' ', zradab_i(ji,nlay_i)
+             write(numout,*) 'icethd_zdf_bl99: zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji) = ',zkappa_i(ji,nlay_i), ' ', zg1, ' ', t_bo_1d(ji)
+           END IF
+         END DO
 
             DO ji = 1, npti
@@ -650,4 +732,12 @@
                   ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
                   zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
+
+                  IF (to_print(ji) == 10) THEN
+                     write(numout,*) 'ztrid for snow = ',ztrid   (ji,2,:)
+                     write(numout,*) 'zindterm for snow = ',zindterm(ji,2)
+                     write(numout,*) 'zindterm using degC = ',zindterm(ji,2) - rt0
+                     write(numout,*) 'nlay_s, nlay_i = ',nlay_s, ' ', nlay_i
+                  ENDIF
+
                   
                   !                            !---------------------!
@@ -662,4 +752,8 @@
                   ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
                   zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+
+                  IF (to_print(ji) == 10 .AND. iconv == 1) THEN
+                     write(numout,*) 'zeta_i(ji,1), zkappa_i(ji,1), ztiold(ji,1), qcn_ice_1d(ji) = ',zeta_i(ji,1), ' ', zkappa_i(ji,1), ' ', ztiold(ji,1), ' ', qcn_ice_1d(ji)
+                  ENDIF                  
                   
                   ! case of only one layer in the ice (surface & ice equations are altered)
@@ -725,5 +819,15 @@
                   ENDIF
                ENDIF
-            END DO
+
+              IF (to_print(ji) == 10) THEN
+                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)
+              END IF
+            END DO
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'icethd_zdf_bl99 8 iter: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
             !
             !--------------------------------------------------------------
@@ -755,6 +859,28 @@
 
          ENDIF ! k_cnd
+
+         DO ji = 1, npti
+           IF (iconv == 1 .AND. to_print(ji) == 10) THEN
+             write(numout,*) 'matrix1 befroe diagonal = ',ztrid(ji,1,1),' ',ztrid(ji,2,1),' ',ztrid(ji,3,1),' ',ztrid(ji,4,1),' ',ztrid(ji,5,1),' ',ztrid(ji,6,1)
+             write(numout,*) 'matrix1 on diagonal = ',ztrid(ji,1,2),' ',ztrid(ji,2,2),' ',ztrid(ji,3,2),' ',ztrid(ji,4,2),' ',ztrid(ji,5,2),' ',ztrid(ji,6,2)
+             write(numout,*) 'matrix1 after diagonal = ',ztrid(ji,1,3),' ',ztrid(ji,2,3),' ',ztrid(ji,3,3),' ',ztrid(ji,4,3),' ',ztrid(ji,5,3),' ',ztrid(ji,6,3)
+             write(numout,*) 'rhs1 for all = ',zindterm(ji,1),' ',zindterm(ji,2),' ',zindterm(ji,3),' ',zindterm(ji,4),' ',zindterm(ji,5),' ',zindterm(ji,6)
+             write(numout,*) 'rhs1 using degC = ',zindterm(ji,1)-rt0,' ',zindterm(ji,2)-rt0,' ',zindterm(ji,3)-rt0,' ',zindterm(ji,4)-rt0,' ',zindterm(ji,5)-rt0,' ',zindterm(ji,6)-rt0
+           END IF
+         END DO
          
       END DO  ! End of the do while iterative procedure
+
+      DO ji = 1, npti
+        IF (to_print(ji) == 10) THEN
+          write(numout,*) 'matrix2 befroe diagonal = ',ztrid(ji,1,1),' ',ztrid(ji,2,1),' ',ztrid(ji,3,1),' ',ztrid(ji,4,1),' ',ztrid(ji,5,1),' ',ztrid(ji,6,1)
+          write(numout,*) 'matrix2 on diagonal = ',ztrid(ji,1,2),' ',ztrid(ji,2,2),' ',ztrid(ji,3,2),' ',ztrid(ji,4,2),' ',ztrid(ji,5,2),' ',ztrid(ji,6,2)
+          write(numout,*) 'matrix2 after diagonal = ',ztrid(ji,1,3),' ',ztrid(ji,2,3),' ',ztrid(ji,3,3),' ',ztrid(ji,4,3),' ',ztrid(ji,5,3),' ',ztrid(ji,6,3)
+          write(numout,*) 'rhs2 for all = ',zindterm(ji,1),' ',zindterm(ji,2),' ',zindterm(ji,3),' ',zindterm(ji,4),' ',zindterm(ji,5),' ',zindterm(ji,6)
+          write(numout,*) 'rhs2 using degC = ',zindterm(ji,1)-rt0,' ',zindterm(ji,2)-rt0,' ',zindterm(ji,3)-rt0,' ',zindterm(ji,4)-rt0,' ',zindterm(ji,5)-rt0,' ',zindterm(ji,6)-rt0
+          write(numout,*) 'number of iterations iconv = ',iconv
+        END IF
+      END DO
+
       
       IF( ln_icectl .AND. lwp ) THEN
@@ -799,6 +925,18 @@
       !
       IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'before ice_var_enthalpy: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO
          
-         CALL ice_var_enthalpy       
+         CALL ice_var_enthalpy    
+
+         DO ji = 1, npti
+           IF (to_print(ji) == 10) THEN
+             write(numout,*) 'after ice_var_enthalpy: t_i_1d(ji,1), e_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', e_i_1d(ji,1)
+           END IF
+         END DO   
          
          ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
@@ -829,4 +967,12 @@
             ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
             hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
+
+            IF (to_print(ji) == 10) THEN
+               write(numout,*)'icethd_zdf_bl99: qcn_ice_top_1d(ji) = ',qcn_ice_top_1d(ji)
+               write(numout,*)'icethd_zdf_bl99: qcn_ice_1d(ji) = ',qcn_ice_1d(ji)
+               write(numout,*)'icethd_zdf_bl99: hfx_err_dif_1d(ji) = ', hfx_err_dif_1d(ji)
+               write(numout,*)'icethd_zdf_bl99: hfx_dif_1d(ji) = ',hfx_dif_1d(ji)
+            ENDIF
+
             !
          END DO
@@ -849,4 +995,10 @@
          ENDIF
          t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
+
+         IF (to_print(ji) == 10) THEN
+            write(numout,*)'icethd_zdf_bl99: cnd_ice_1d(ji), h_s_1d(ji), zkappa_s(ji,0) = ',cnd_ice_1d(ji), ' ', h_s_1d(ji), ' ', zkappa_s(ji,0)
+            write(numout,*)'icethd_zdf_bl99: ztcond_i(ji,0), h_i_1d(ji), zkappa_i(ji,0) = ',ztcond_i(ji,0), ' ', h_i_1d(ji), ' ', zkappa_i(ji,0)
+            write(numout,*)'icethd_zdf_bl99: t1_ice_1d(ji), isnow(ji), t_s_1d(ji,1), t_i_1d(ji,1) = ',t1_ice_1d(ji), ' ', isnow(ji), ' ', t_s_1d(ji,1), ' ', t_i_1d(ji,1)
+         ENDIF
       END DO
       !
Index: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90	(revision 12368)
+++ NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90	(revision 12369)
@@ -219,4 +219,8 @@
          END WHERE
       ENDIF
+
+      write(numout,*) 'icevar: sv_i(3,4,1) = ', sv_i(3,4,1)
+      write(numout,*) 'icevar: z1_v_i(3,4,1) = ', z1_v_i(3,4,1)
+
       CALL ice_var_salprof   ! salinity profile
 
@@ -241,4 +245,9 @@
                      t_i(ji,jj,jk,jl) = rt0
                   ENDIF
+
+                  If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
+                     write(numout,*)'ice_var_glo2eqv 0: t_i, e_i = ',t_i(ji,jj,jk,jl) - rt0, ' ', e_i (ji,jj,jk,jl)
+                  END IF
+
                END DO
             END DO
@@ -257,4 +266,8 @@
             t_s(:,:,jk,:) = rt0
          END WHERE
+
+         If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
+             write(numout,*)'ice_var_glo2eqv 0: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
+         END IF
       END DO
       !
@@ -327,4 +340,5 @@
             END DO
          END DO
+
          !                                      ! Slope of the linear profile 
          WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
@@ -342,4 +356,5 @@
             END DO
          END DO
+         
          !
          ! Computation of the profile
@@ -352,4 +367,8 @@
                      zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
                      sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
+
+                     IF ( jl == 1 .AND. jj == 26 .AND. ji == 42 ) THEN
+                       write(numout,*) 'icevar: jk, sz_i(ji,jj,jk,jl), s_i(ji,jj,jl), weights(jk) = ',jk, ' ',sz_i(ji,jj,jk,jl), ' ',s_i(ji,jj,jl), ' ',weights(jk)
+                     ENDIF
                   END DO
                END DO
@@ -497,4 +516,8 @@
                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
                   t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
+
+                  If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
+                    write(numout,*)'ice_var_glo2eqv 1: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
+                  END IF
                END DO
             END DO
