Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 3609)
@@ -32,4 +32,6 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
+   USE icb_oce         ! Icebergs
+   USE icbdia          ! Iceberg budgets
    USE sbcssr          ! restoring term toward SST/SSS climatology
    USE phycst          ! physical constants
@@ -59,4 +61,5 @@
 
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
+   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
    INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
@@ -65,4 +68,5 @@
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
 
    !! * Substitutions
@@ -234,4 +238,5 @@
       INTEGER  ::   ierr                                     ! error code return from allocation
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
       !!
@@ -320,4 +325,27 @@
          CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
          CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
+         !
+         IF( ln_icebergs ) THEN
+            !
+            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after 
+            !! that routine is called from nemogcm, so do it here immediately before its needed
+            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
+            IF( lk_mpp )   CALL mpp_sum( ierror )
+            IF( ierror /= 0 ) THEN
+               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
+               RETURN
+            ENDIF
+            !
+            !! iceberg vertical coordinate is class number
+            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
+               &           "number", nclasses, class_num, nb_T )
+            !
+            !! each class just needs the surface index pattern
+            ndim_bT = 3
+            DO jn = 1,nclasses
+               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
+            ENDDO
+            !
+         ENDIF
 
          ! Define the U grid FILE ( nid_U )
@@ -401,4 +429,38 @@
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+!
+         IF( ln_icebergs ) THEN
+            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
+               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
+            IF( ln_bergdia ) THEN
+               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
+                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
+                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
+            ENDIF
+         ENDIF
+
 #if ! defined key_coupled 
          CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
@@ -555,4 +617,28 @@
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
+!
+      IF( ln_icebergs ) THEN
+         !
+         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT )  
+         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
+         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT )  
+         !
+         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
+         !
+         IF( ln_bergdia ) THEN
+            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   )  
+            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   )  
+            !
+            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
+         ENDIF
+      ENDIF
+
 #if ! defined key_coupled
       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 3609)
@@ -256,5 +256,5 @@
       END DO
       ! applied the lateral boundary conditions
-      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
+      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
 
 #if defined key_agrif
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 3609)
@@ -0,0 +1,202 @@
+MODULE icb_oce
+
+   !!======================================================================
+   !!                       ***  MODULE  icb_oce  ***
+   !! Icebergs:  declare variables for iceberg tracking
+   !!======================================================================
+   !! History :  3.3  !  2010-01  (T. Martin & A. Adcroft)  Original code
+   !!             -   !  2011-03  (G. Madec)  Part conversion to NEMO form
+   !!             -   !                       Removal of mapping from another grid
+   !!             -   !  2011-04  (S. Alderson) Extensive rewrite ; Split into separate modules
+   !!----------------------------------------------------------------------
+   !!
+   !! Track Icebergs as Lagrangian objects within the model domain
+   !! Interaction with the other model variables through 'icebergs_gridded'
+   !!
+   !! A single iceberg is held as an instance of type 'iceberg'
+   !! This type defines a linked list, so each instance contains a pointer 
+   !! to the previous and next icebergs in the list
+   !!
+   !! Type 'icebergs' is a convenience container for all relevant arrays
+   !! It contains one pointer to an 'iceberg' instance representing all icebergs in the processor
+   !!
+   !! Each iceberg has a position represented as a real cartesian coordinate which is 
+   !! fractional grid cell, centred on T-points; so an iceberg position of (1.0,1.0) lies 
+   !! exactly on the first T-point and the T-cell spans 0.5 to 1.5 in each direction
+   !!
+   !! Each iceberg is assigned a unique id even in MPI
+   !! This consists of an array of integers: the first element is used to label, the second
+   !! and subsequent elements are used to count the number of times the first element wraps
+   !! around all possible values within the valid size for this datatype.
+   !! Labelling is done by starting the first label in each processor (even when only one)
+   !! as narea, and then incrementing by jpnij (i.e. the total number of processors.
+   !! This means that the source processor for each iceberg can be identified by arithmetic
+   !! modulo jpnij.
+   !! 
+   !!----------------------------------------------------------------------
+   USE par_oce   ! ocean parameters
+   USE lib_mpp   ! MPP library
+   USE fldread   ! read input fields (FLD type)
+
+   IMPLICIT NONE
+   PUBLIC
+
+   PUBLIC   icb_alloc   ! routine called by icb_init in icbini.F90 module
+
+   INTEGER, PUBLIC, PARAMETER ::   nclasses = 10   !: Number of ice bergs classes
+   INTEGER, PUBLIC, PARAMETER ::   nkounts  =  3   !: Number of integers combined for unique naming
+
+   TYPE, PUBLIC ::   icebergs_gridded   !: various icebergs properties on model grid
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   calving         ! Calving mass rate  (into stored ice)         [kg/s]
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   calving_hflx    ! Calving heat flux [heat content of calving]  [W/m2]
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   floating_melt   ! Net melting rate to icebergs + bits      [kg/s/m^2]
+      INTEGER , DIMENSION(:,:)  , ALLOCATABLE ::   maxclass        ! maximum class number at calving source point
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   tmp             ! Temporary work space
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   stored_ice      ! Accumulated ice mass flux at calving locations [kg]
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   stored_heat     ! Heat content of stored ice                      [J]
+   END TYPE icebergs_gridded
+
+   TYPE, PUBLIC ::   point              !: properties of an individual iceberg (position, mass, size, etc...)
+      INTEGER  ::   year
+      REAL(wp) ::   xi , yj                                              ! iceberg coordinates in the (i,j) referential (global)
+      REAL(wp) ::   e1 , e2                                              ! horizontal scale factors at the iceberg position
+      REAL(wp) ::   lon, lat, day                                        ! geographic position
+      REAL(wp) ::   mass, thickness, width, length, uvel, vvel           ! iceberg physical properties
+      REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment 
+      REAL(wp) ::   mass_of_bits, heat_density
+   END TYPE point
+
+   TYPE, PUBLIC ::   iceberg            !: linked list defining all the icebergs present in the model domain
+      TYPE(iceberg), POINTER      ::   prev=>NULL(), next=>NULL()   ! pointers to previous and next unique icebergs in linked list
+      INTEGER, DIMENSION(nkounts) ::   number                       ! variables which do not change for this iceberg
+      REAL(wp)                    ::   mass_scaling                 !    -        -            -                -  
+      TYPE(point), POINTER        ::   current_point => NULL()      ! variables which change with time are held in a separate type
+   END TYPE iceberg
+
+
+   TYPE(icebergs_gridded), POINTER ::   berg_grid                 !: master instance of gridded iceberg type
+   TYPE(iceberg)         , POINTER ::   first_berg => NULL()      !: master instance of linked list iceberg type
+
+   !                                                             !!! parameters controlling iceberg characteristics and modelling
+   REAL(wp)                            ::   berg_dt                   !: Time-step between iceberg CALLs (should make adaptive?)
+   REAL(wp), DIMENSION(:), ALLOCATABLE ::   first_width, first_length !: 
+   LOGICAL                             ::   l_restarted_bergs=.FALSE.  ! Indicate whether we read state from a restart or not
+   !                                                               ! arbitrary numbers for diawri entry
+   REAL(wp), DIMENSION(nclasses), PUBLIC ::   class_num=(/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 /)
+
+   ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position
+   ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e
+#if defined key_lim2 || defined key_lim3
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ui_e, vi_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e
+#endif
+
+   !!gm almost all those PARAM ARE defined in NEMO
+   REAL(wp), PUBLIC, PARAMETER :: pp_rho_ice      = 916.7_wp   !: Density of fresh ice   @ 0oC [kg/m^3]
+   REAL(wp), PUBLIC, PARAMETER :: pp_rho_water    = 999.8_wp   !: Density of fresh water @ 0oC [kg/m^3]
+   REAL(wp), PUBLIC, PARAMETER :: pp_rho_air      = 1.1_wp     !: Density of air         @ 0oC [kg/m^3]
+   REAL(wp), PUBLIC, PARAMETER :: pp_rho_seawater = 1025._wp   !: Approx. density of surface sea water @ 0oC [kg/m^3]
+   !!gm end
+   REAL(wp), PUBLIC, PARAMETER :: pp_Cd_av = 1.3_wp      !: (Vertical) Drag coefficient between bergs and atmos 
+   REAL(wp), PUBLIC, PARAMETER :: pp_Cd_ah = 0.0055_wp   !: (lateral ) Drag coefficient between bergs and atmos 
+   REAL(wp), PUBLIC, PARAMETER :: pp_Cd_wv = 0.9_wp      !: (Vertical) Drag coefficient between bergs and ocean
+   REAL(wp), PUBLIC, PARAMETER :: pp_Cd_wh = 0.0012_wp   !: (lateral ) Drag coefficient between bergs and ocean 
+   REAL(wp), PUBLIC, PARAMETER :: pp_Cd_iv = 0.9_wp      !: (Vertical) Drag coefficient between bergs and sea-ice
+!TOM> no horizontal drag for sea ice! real, PARAMETER :: pp_Cd_ih=0.0012 ! (lateral) Drag coeff. between bergs and sea-ice
+
+   !                                                         !!* namberg namelist parameters (and defaults) **
+   LOGICAL , PUBLIC ::   ln_bergdia               = .true.    !: Calculate budgets
+   INTEGER , PUBLIC ::   nn_verbose_level         = 0         !: Turn on debugging when level > 0
+   INTEGER , PUBLIC ::   nn_test_icebergs         = 0         !: Create icebergs in absence of a restart file from the supplied class nb
+   REAL(wp), PUBLIC, DIMENSION(4) ::   rn_test_box = (/ 0._wp, 1._wp, 0._wp, 1._wp /)   !: lon1,lon2,lat1,lat2 box to create them in
+   INTEGER , PUBLIC ::   nn_sample_rate           = 0         !: Timesteps between sampling of position for trajectory storage
+   INTEGER , PUBLIC ::   nn_verbose_write         = 15        !: timesteps between verbose messages
+   REAL(wp), PUBLIC ::   rn_rho_bergs             = 850._wp   !: Density of icebergs
+   REAL(wp), PUBLIC ::   rn_LoW_ratio             = 1.5_wp    !: Initial ratio L/W for newly calved icebergs
+   REAL(wp), PUBLIC ::   rn_bits_erosion_fraction = 0.        !: Fraction of erosion melt flux to divert to bergy bits
+   REAL(wp), PUBLIC ::   rn_sicn_shift            = 0._wp     !: Shift of sea-ice concentration in erosion flux modulation (0<sicn_shift<1)
+   LOGICAL , PUBLIC ::   ln_operator_splitting    = .true.    !: Use first order operator splitting for thermodynamics
+   LOGICAL , PUBLIC ::   ln_passive_mode          = .false.   !: iceberg - ocean decoupling
+   LOGICAL , PUBLIC ::   ln_time_average_weight   = .false.   !: Time average the weight on the ocean    !!gm I don't understand that !
+   REAL(wp), PUBLIC ::   rn_speed_limit           = 0._wp     !: CFL speed limit for a berg
+   !
+   !                                     ! Mass thresholds between iceberg classes [kg]
+   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_initial_mass      = (/8.8e07, 4.1e08, 3.3e09, 1.8e10, 3.8e10,   &
+      &                                                                7.5e10, 1.2e11, 2.2e11, 3.9e11, 7.4e11/)
+   !                                     ! Fraction of calving to apply to this class [non-dim]
+   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_distribution      = (/0.25  , 0.12  , 0.15  , 0.18  , 0.12  ,   &
+      &                                                                0.07  , 0.03  , 0.03  , 0.03  , 0.02  /)
+   !                                     ! Ratio between effective and real iceberg mass (non-dim)
+   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_mass_scaling      = (/2000. ,  200. ,   50. ,   20. ,   10. ,   &
+      &                                                                   5. ,    2. ,    1. ,    1. ,    1. /)
+   !                                     ! Total thickness of newly calved bergs [m]
+   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_initial_thickness = (/  40. ,   67. ,  133. ,  175. ,  250. ,   &
+      &                                                                 250. ,  250. ,  250. ,  250. ,  250./)
+
+   ! Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   src_calving, src_calving_hflx    !: accumulate input ice
+   INTEGER , PUBLIC             , SAVE                     ::   numicb                           !: iceberg IO
+   INTEGER , PUBLIC             , SAVE, DIMENSION(nkounts) ::   num_bergs                        !: iceberg counter
+   INTEGER , PUBLIC             , SAVE                     ::   nicbdi, nicbei, nicbdj, nicbej   !: processor bounds
+   REAL(wp), PUBLIC             , SAVE                     ::   ricb_left, ricb_right            !: cyclical bounds
+   INTEGER , PUBLIC             , SAVE                     ::   nicbpack                         !: packing integer
+   INTEGER , PUBLIC             , SAVE                     ::   nktberg, nknberg                 !: helpers
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldpts                       !: nfold packed points
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbflddest                      !: nfold destination proc
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   nicbfldproc                      !: nfold destination proc
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: griddata                           !: work array for icbrst
+
+   TYPE(FLD), PUBLIC, ALLOCATABLE     , DIMENSION(:)       ::   sf_icb   !: structure: file information, fields read
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id: sbc_oce.F90 3340 2012-04-02 11:05:35Z sga $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+   
+   INTEGER FUNCTION icb_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_alloc  ***
+      !!----------------------------------------------------------------------
+      INTEGER ::   ill
+      !!----------------------------------------------------------------------
+      !
+      icb_alloc = 0
+      ALLOCATE( berg_grid                      ,                                               &
+         &      berg_grid%calving    (jpi,jpj) , berg_grid%calving_hflx (jpi,jpj)          ,   &
+         &      berg_grid%stored_heat(jpi,jpj) , berg_grid%floating_melt(jpi,jpj)          ,   &
+         &      berg_grid%maxclass   (jpi,jpj) , berg_grid%stored_ice   (jpi,jpj,nclasses) ,   &
+         &      berg_grid%tmp        (jpi,jpj) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      !
+      ! expanded arrays for bilinear interpolation
+      ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   &
+         &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   &
+#if defined key_lim2 || defined key_lim3
+         &      ui_e(0:jpi+1,0:jpj+1) ,                            &
+         &      vi_e(0:jpi+1,0:jpj+1) ,                            &
+#endif
+         &      ff_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   &
+         &      first_width(nclasses) , first_length(nclasses) ,   &
+         &      src_calving (jpi,jpj) ,                            &
+         &      src_calving_hflx(jpi,jpj) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+
+      ALLOCATE( nicbfldpts(jpi) , nicbflddest(jpi) , nicbfldproc(jpni) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+
+      ALLOCATE( griddata(jpi,jpj,1), STAT=ill )
+      icb_alloc = icb_alloc + ill
+
+      IF( lk_mpp        )   CALL mpp_sum ( icb_alloc )
+      IF( icb_alloc > 0 )   CALL ctl_warn('icb_alloc: allocation of arrays failed')
+      !
+   END FUNCTION icb_alloc
+
+   !!======================================================================
+END MODULE icb_oce
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 3609)
@@ -0,0 +1,183 @@
+MODULE icbclv
+
+   !!======================================================================
+   !!                       ***  MODULE  icbclv  ***
+   !! Icebergs:  calving routines for iceberg calving
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !  2011-05  (Alderson)       budgets into separate module
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_clv_flx   : transfer input flux of ice into iceberg classes
+   !!   icb_clv       : calve icebergs from stored ice
+   !!----------------------------------------------------------------------
+   USE par_oce        ! NEMO parameters
+   USE dom_oce        ! NEMO ocean domain
+   USE phycst         ! NEMO physical constants
+   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
+   USE lbclnk         ! NEMO boundary exchanges for gridded data
+
+   USE icb_oce        ! iceberg variables
+   USE icbdia         ! iceberg diagnostics
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_clv_flx  ! routine called in icbstp.F90 module
+   PUBLIC   icb_clv      ! routine called in icbstp.F90 module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_clv_flx( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_clv_flx  ***
+      !!
+      !! ** Purpose :   accumulate ice available for calving into class arrays
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)             :: kt
+      !
+      REAL(wp)                        :: zcalving_used, zdist, zfact
+      INTEGER                         :: jn, ji, jj                    ! loop counters
+      INTEGER                         :: imx                           ! temporary integer for max berg class
+      LOGICAL, SAVE                   :: ll_first_call = .TRUE.
+      !!----------------------------------------------------------------------
+      !
+      ! Adapt calving flux and calving heat flux from coupler for use here
+      ! Use interior mask: so no bergs in overlap areas and convert from km^3/year to kg/s
+      ! this assumes that input is given as equivalent water flux so that pure water density is appropriate
+
+      zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 1000._wp
+      berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact
+
+      ! Heat in units of W/m2, and mask (just in case)
+      berg_grid%calving_hflx(:,:) = src_calving_hflx(:,:) * tmask_i(:,:)
+
+      IF( ll_first_call .AND. .NOT. l_restarted_bergs) THEN      ! This is a hack to simplify initialization
+         ll_first_call = .FALSE.
+         !do jn=1, nclasses
+         !  where (berg_grid%calving==0.) berg_grid%stored_ice(:,:,jn)=0.
+         !end do
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF( berg_grid%calving(ji,jj) /= 0._wp )                                  &    ! Need units of J
+                  berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) *         &  ! initial stored ice in kg
+                                         berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) /   &  ! J/s/m2 x m^2 = J/s
+                                         berg_grid%calving(ji,jj)                            ! /calving in kg/s
+            END DO
+         END DO
+      ENDIF
+
+      ! assume that all calving flux must be distributed even if distribution array does not sum
+      ! to one - this may not be what is intended, but it's what you've got
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            imx = berg_grid%maxclass(ji,jj)
+            zdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:imx) )
+            DO jn = 1, imx
+               berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) + &
+                                          berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist
+            END DO
+         END DO
+      END DO
+
+      ! before changing the calving, save the amount we're about to use and do budget
+      zcalving_used = SUM( berg_grid%calving(:,:) )
+      berg_grid%tmp(:,:) = berg_dt * berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:)
+      berg_grid%stored_heat (:,:) = berg_grid%stored_heat (:,:) + berg_grid%tmp(:,:)
+      CALL icb_dia_income( kt,  zcalving_used, berg_grid%tmp )
+      !
+   END SUBROUTINE icb_clv_flx
+
+   SUBROUTINE icb_clv()
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_clv  ***
+      !!
+      !! ** Purpose :   This routine takes a stored ice field and calves to the ocean,
+      !!                so the gridded array stored_ice has only non-zero entries at selected
+      !!                wet points adjacent to known land based calving points
+      !!
+      !! ** method  : - Look at each grid point and see if there's enough for each size class to calve
+      !!                If there is, a new iceberg is calved.  This happens in the order determined by
+      !!                the class definition arrays (which in the default case is smallest first)
+      !!                Note that only the non-overlapping part of the processor where icebergs are allowed
+      !!                is considered
+      !!----------------------------------------------------------------------
+      INTEGER       ::   ji, jj, jn   ! dummy loop indices
+      INTEGER       ::   icnt, icntmax
+      TYPE(iceberg) ::   newberg
+      TYPE(point)   ::   newpt
+      REAL(wp)      ::   zday, zcalved_to_berg, zheat_to_berg
+      !!----------------------------------------------------------------------
+      !
+      icntmax = 0
+      zday    = REAL(nday_year,wp) + REAL(nsec_day,wp)/86400.0_wp
+      !
+      DO jn = 1, nclasses
+         DO jj = nicbdj, nicbej
+            DO ji = nicbdi, nicbei
+               !
+               icnt = 0
+               !
+               DO WHILE (berg_grid%stored_ice(ji,jj,jn) >= rn_initial_mass(jn) * rn_mass_scaling(jn) )
+                  !
+                  newpt%lon = glamt(ji,jj)         ! at t-point (centre of the cell)
+                  newpt%lat = gphit(ji,jj)
+                  newpt%xi  = REAL( mig(ji), wp )
+                  newpt%yj  = REAL( mjg(jj), wp )
+                  !
+                  newpt%uvel = 0._wp               ! initially at rest
+                  newpt%vvel = 0._wp
+                  !                                ! set berg characteristics
+                  newpt%mass           = rn_initial_mass     (jn)
+                  newpt%thickness      = rn_initial_thickness(jn)
+                  newpt%width          = first_width         (jn)
+                  newpt%length         = first_length        (jn)
+                  newberg%mass_scaling = rn_mass_scaling     (jn)
+                  newpt%mass_of_bits   = 0._wp                          ! no bergy
+                  !
+                  newpt%year   = nyear
+                  newpt%day    = zday
+                  newpt%heat_density = berg_grid%stored_heat(ji,jj) / berg_grid%stored_ice(ji,jj,jn)   ! This is in J/kg
+                  !
+                  CALL icb_utl_incr()
+                  newberg%number(:) = num_bergs(:)
+                  !
+                  CALL icb_utl_add( newberg, newpt )
+                  !
+                  zcalved_to_berg = rn_initial_mass(jn) * rn_mass_scaling(jn)           ! Units of kg
+                  !                                ! Heat content
+                  zheat_to_berg           = zcalved_to_berg * newpt%heat_density             ! Units of J
+                  berg_grid%stored_heat(ji,jj) = berg_grid%stored_heat(ji,jj) - zheat_to_berg
+                  !                                ! Stored mass
+                  berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) - zcalved_to_berg
+                  !
+                  icnt = icnt + 1
+                  !
+                  CALL icb_dia_calve(ji, jj, jn,  zcalved_to_berg, zheat_to_berg )
+               END DO
+               icntmax = MAX( icntmax, icnt )
+            END DO
+         END DO
+      END DO
+      !
+      DO jn = 1,nclasses
+         CALL lbc_lnk( berg_grid%stored_ice(:,:,jn), 'T', 1._wp )
+      END DO
+      CALL lbc_lnk( berg_grid%stored_heat, 'T', 1._wp )
+      !
+      IF( nn_verbose_level > 0 .AND. icntmax > 1 )   WRITE(numicb,*) 'icb_clv: icnt=', icnt,' on', narea
+      !
+   END SUBROUTINE  icb_clv
+
+   !!======================================================================
+END MODULE icbclv
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 3609)
@@ -0,0 +1,580 @@
+MODULE icbdia
+
+   !!======================================================================
+   !!                       ***  MODULE  icbdia  ***
+   !! Icebergs:  initialise variables for iceberg budgets and diagnostics
+   !!======================================================================
+   !! History : 3.3 !  2010-01  (Martin, Adcroft) Original code
+   !!            -  !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -  !                            Removal of mapping from another grid
+   !!            -  !  2011-04  (Alderson)       Split into separate modules
+   !!            -  !  2011-05  (Alderson)       Budgets are now all here with lots
+   !!            -  !                            of silly routines to call to get values in
+   !!            -  !                            from the right points in the code
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !! icb_dia_init : initialise iceberg budgeting
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE dom_oce        ! ocean domain
+   USE in_out_manager ! nemo IO
+   USE lib_mpp        ! MPP library
+   USE iom            ! I/O library
+   USE icb_oce        ! iceberg variables
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_dia_init      ! routine called in icbini.F90 module
+   PUBLIC   icb_dia           ! routine called in icbstp.F90 module
+   PUBLIC   icb_dia_step      ! routine called in icbstp.F90 module
+   PUBLIC   icb_dia_put       ! routine called in icbstp.F90 module
+   PUBLIC   icb_dia_melt      ! routine called in icbthm.F90 module
+   PUBLIC   icb_dia_size      ! routine called in icbthm.F90 module
+   PUBLIC   icb_dia_speed     ! routine called in icbdyn.F90 module
+   PUBLIC   icb_dia_calve     ! routine called in icbclv.F90 module
+   PUBLIC   icb_dia_income    ! routine called in icbclv.F90 module
+
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   berg_melt       ! Melting+erosion rate of icebergs     [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   buoy_melt       ! Buoyancy component of melting rate   [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   eros_melt       ! Erosion component of melting rate    [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   conv_melt       ! Convective component of melting rate [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   bits_src        ! Mass flux from berg erosion into bergy bits [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   bits_melt       ! Melting rate of bergy bits           [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   bits_mass       ! Mass distribution of bergy bits      [kg/s/m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   virtual_area    ! Virtual surface coverage by icebergs [m2]
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE, PUBLIC  ::   berg_mass       ! Mass distribution                    [kg/m2]
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, PUBLIC  ::   real_calving    ! Calving rate into iceberg class at
+   !                                                                          ! calving locations                    [kg/s]
+   
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   tmpc                     ! Temporary work space
+   REAL(wp), DIMENSION(:)    , ALLOCATABLE ::   rsumbuf                  ! Temporary work space to reduce mpp exchanges
+   INTEGER , DIMENSION(:)    , ALLOCATABLE ::   nsumbuf                  ! Temporary work space to reduce mpp exchanges
+
+   REAL(wp)                      ::  berg_melt_net
+   REAL(wp)                      ::  bits_src_net
+   REAL(wp)                      ::  bits_melt_net
+   REAL(wp)                      ::  bits_mass_start     , bits_mass_end
+   REAL(wp)                      ::  floating_heat_start , floating_heat_end
+   REAL(wp)                      ::  floating_mass_start , floating_mass_end
+   REAL(wp)                      ::  bergs_mass_start    , bergs_mass_end
+   REAL(wp)                      ::  stored_start        , stored_heat_start
+   REAL(wp)                      ::  stored_end          , stored_heat_end
+   REAL(wp)                      ::  calving_src_net     , calving_out_net
+   REAL(wp)                      ::  calving_src_heat_net, calving_out_heat_net
+   REAL(wp)                      ::  calving_src_heat_used_net
+   REAL(wp)                      ::  calving_rcv_net  , calving_ret_net  , calving_used_net
+   REAL(wp)                      ::  heat_to_bergs_net, heat_to_ocean_net, melt_net
+   REAL(wp)                      ::  calving_to_bergs_net
+
+   INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved
+   INTEGER                       ::  nbergs_melted
+   INTEGER                       ::  nspeeding_tickets
+   INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_dia_init( )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+
+      ALLOCATE( berg_melt    (jpi,jpj)   )           ;   berg_melt   (:,:)   = 0._wp
+      ALLOCATE( buoy_melt    (jpi,jpj)   )           ;   buoy_melt   (:,:)   = 0._wp
+      ALLOCATE( eros_melt    (jpi,jpj)   )           ;   eros_melt   (:,:)   = 0._wp
+      ALLOCATE( conv_melt    (jpi,jpj)   )           ;   conv_melt   (:,:)   = 0._wp
+      ALLOCATE( bits_src     (jpi,jpj)   )           ;   bits_src    (:,:)   = 0._wp
+      ALLOCATE( bits_melt    (jpi,jpj)   )           ;   bits_melt   (:,:)   = 0._wp
+      ALLOCATE( bits_mass    (jpi,jpj)   )           ;   bits_mass   (:,:)   = 0._wp
+      ALLOCATE( virtual_area (jpi,jpj)   )           ;   virtual_area(:,:)   = 0._wp
+      ALLOCATE( berg_mass    (jpi,jpj)   )           ;   berg_mass   (:,:)   = 0._wp
+      ALLOCATE( real_calving (jpi,jpj,nclasses) )    ;   real_calving(:,:,:) = 0._wp
+      ALLOCATE( tmpc(jpi,jpj) )                      ;   tmpc        (:,:)   = 0._wp
+
+      nbergs_start              = 0
+      nbergs_end                = 0
+      stored_end                = 0._wp
+      nbergs_start              = 0._wp
+      stored_start              = 0._wp
+      nbergs_melted             = 0
+      nbergs_calved             = 0
+      nbergs_calved_by_class(:) = 0
+      nspeeding_tickets         = 0
+      stored_heat_end           = 0._wp
+      floating_heat_end         = 0._wp
+      floating_mass_end         = 0._wp
+      bergs_mass_end            = 0._wp
+      bits_mass_end             = 0._wp
+      stored_heat_start         = 0._wp
+      floating_heat_start       = 0._wp
+      floating_mass_start       = 0._wp
+      bergs_mass_start          = 0._wp
+      bits_mass_start           = 0._wp
+      bits_mass_end             = 0._wp
+      calving_used_net          = 0._wp
+      calving_to_bergs_net      = 0._wp
+      heat_to_bergs_net         = 0._wp
+      heat_to_ocean_net         = 0._wp
+      calving_rcv_net           = 0._wp
+      calving_ret_net           = 0._wp
+      calving_src_net           = 0._wp
+      calving_out_net           = 0._wp
+      calving_src_heat_net      = 0._wp
+      calving_src_heat_used_net = 0._wp
+      calving_out_heat_net      = 0._wp
+      melt_net                  = 0._wp
+      berg_melt_net             = 0._wp
+      bits_melt_net             = 0._wp
+      bits_src_net              = 0._wp
+
+      floating_mass_start       = icb_utl_mass( first_berg )
+      bergs_mass_start          = icb_utl_mass( first_berg, justbergs=.true. )
+      bits_mass_start           = icb_utl_mass( first_berg, justbits=.true. )
+      IF( lk_mpp ) THEN
+         ALLOCATE( rsumbuf(23) )          ; rsumbuf(:) = 0._wp
+         ALLOCATE( nsumbuf(4+nclasses) )  ; nsumbuf(:) = 0
+         rsumbuf(1) = floating_mass_start
+         rsumbuf(2) = bergs_mass_start
+         rsumbuf(3) = bits_mass_start
+         CALL mpp_sum( rsumbuf(1:3), 3 )
+         floating_mass_start = rsumbuf(1)
+         bergs_mass_start = rsumbuf(2)
+         bits_mass_start = rsumbuf(3)
+      ENDIF
+      !
+   END SUBROUTINE icb_dia_init
+
+
+   SUBROUTINE icb_dia( ld_budge )
+      !!----------------------------------------------------------------------
+      !! sum all the things we've accumulated so far in the current processor
+      !! in MPP case then add these sums across all processors
+      !! for this we pack variables into buffer so we only need one mpp_sum
+      !!----------------------------------------------------------------------
+      LOGICAL, INTENT(in) ::   ld_budge
+      !
+      INTEGER             ::   ik
+      REAL(wp)            ::   zunused_calving, ztmpsum, zgrdd_berg_mass, zgrdd_bits_mass
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia )   RETURN
+
+      zunused_calving      = SUM( berg_grid%calving(:,:) )
+      ztmpsum              = SUM( berg_grid%floating_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      melt_net             = melt_net + ztmpsum * berg_dt
+      calving_out_net      = calving_out_net + ( zunused_calving + ztmpsum ) * berg_dt
+      ztmpsum              = SUM( berg_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      berg_melt_net        = berg_melt_net + ztmpsum * berg_dt
+      ztmpsum              = SUM( bits_src(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      bits_src_net         = bits_src_net + ztmpsum * berg_dt
+      ztmpsum              = SUM( bits_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      bits_melt_net        = bits_melt_net + ztmpsum * berg_dt
+      ztmpsum              = SUM( src_calving(:,:) * tmask_i(:,:) )
+      calving_ret_net      = calving_ret_net + ztmpsum * berg_dt
+      ztmpsum              = SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      calving_out_heat_net = calving_out_heat_net + ztmpsum * berg_dt   ! Units of J
+
+      IF( ld_budge ) THEN
+         stored_end        = SUM( berg_grid%stored_ice(:,:,:) )
+         stored_heat_end   = SUM( berg_grid%stored_heat(:,:) )
+         floating_mass_end = icb_utl_mass( first_berg )
+         bergs_mass_end    = icb_utl_mass( first_berg,justbergs=.true. )
+         bits_mass_end     = icb_utl_mass( first_berg,justbits=.true. )
+         floating_heat_end = icb_utl_heat( first_berg )
+
+         nbergs_end        = icb_utl_count()
+         zgrdd_berg_mass   = SUM( berg_mass (:,:)*e1e2t(:,:)*tmask_i(:,:) )
+         zgrdd_bits_mass   = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) )
+
+         IF( lk_mpp ) THEN
+            rsumbuf( 1) = stored_end
+            rsumbuf( 2) = stored_heat_end
+            rsumbuf( 3) = floating_mass_end
+            rsumbuf( 4) = bergs_mass_end
+            rsumbuf( 5) = bits_mass_end
+            rsumbuf( 6) = floating_heat_end
+            rsumbuf( 7) = calving_ret_net
+            rsumbuf( 8) = calving_out_net
+            rsumbuf( 9) = calving_rcv_net
+            rsumbuf(10) = calving_src_net
+            rsumbuf(11) = calving_src_heat_net
+            rsumbuf(12) = calving_src_heat_used_net
+            rsumbuf(13) = calving_out_heat_net
+            rsumbuf(14) = calving_used_net
+            rsumbuf(15) = calving_to_bergs_net
+            rsumbuf(16) = heat_to_bergs_net
+            rsumbuf(17) = heat_to_ocean_net
+            rsumbuf(18) = melt_net
+            rsumbuf(19) = berg_melt_net
+            rsumbuf(20) = bits_src_net
+            rsumbuf(21) = bits_melt_net
+            rsumbuf(22) = zgrdd_berg_mass
+            rsumbuf(23) = zgrdd_bits_mass
+
+            CALL mpp_sum( rsumbuf(1:23), 23)
+
+            stored_end                = rsumbuf( 1)
+            stored_heat_end           = rsumbuf( 2)
+            floating_mass_end         = rsumbuf( 3)
+            bergs_mass_end            = rsumbuf( 4)
+            bits_mass_end             = rsumbuf( 5)
+            floating_heat_end         = rsumbuf( 6)
+            calving_ret_net           = rsumbuf( 7)
+            calving_out_net           = rsumbuf( 8)
+            calving_rcv_net           = rsumbuf( 9)
+            calving_src_net           = rsumbuf(10)
+            calving_src_heat_net      = rsumbuf(11)
+            calving_src_heat_used_net = rsumbuf(12)
+            calving_out_heat_net      = rsumbuf(13)
+            calving_used_net          = rsumbuf(14)
+            calving_to_bergs_net      = rsumbuf(15)
+            heat_to_bergs_net         = rsumbuf(16)
+            heat_to_ocean_net         = rsumbuf(17)
+            melt_net                  = rsumbuf(18)
+            berg_melt_net             = rsumbuf(19)
+            bits_src_net              = rsumbuf(20)
+            bits_melt_net             = rsumbuf(21)
+            zgrdd_berg_mass           = rsumbuf(22)
+            zgrdd_bits_mass           = rsumbuf(23)
+
+            nsumbuf(1) = nbergs_end
+            nsumbuf(2) = nbergs_calved
+            nsumbuf(3) = nbergs_melted
+            nsumbuf(4) = nspeeding_tickets
+            DO ik = 1, nclasses
+               nsumbuf(4+ik) = nbergs_calved_by_class(ik)
+            END DO
+
+            CALL mpp_sum( nsumbuf(1:nclasses+4), nclasses+4 )
+
+            nbergs_end        = nsumbuf(1)
+            nbergs_calved     = nsumbuf(2)
+            nbergs_melted     = nsumbuf(3)
+            nspeeding_tickets = nsumbuf(4)
+            DO ik = 1,nclasses
+               nbergs_calved_by_class(ik)= nsumbuf(4+ik)
+            ENDDO
+
+         ENDIF
+
+         CALL report_state( 'stored ice','kg','',stored_start,'',stored_end,'')
+         CALL report_state( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end)
+         CALL report_state( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'')
+         CALL report_state( 'bits','kg','',bits_mass_start,'',bits_mass_end,'')
+         CALL report_istate( 'berg #','',nbergs_start,'',nbergs_end,'')
+         CALL report_ibudget( 'berg #','calved',nbergs_calved, &
+                                       'melted',nbergs_melted, &
+                                       '#',nbergs_start,nbergs_end)
+         CALL report_budget( 'stored mass','kg','calving used',calving_used_net, &
+                                           'bergs',calving_to_bergs_net, &
+                                           'stored mass',stored_start,stored_end)
+         CALL report_budget( 'floating mass','kg','calving used',calving_to_bergs_net, &
+                                             'bergs',melt_net, &
+                                             'stored mass',floating_mass_start,floating_mass_end)
+         CALL report_budget( 'berg mass','kg','calving',calving_to_bergs_net, &
+                                              'melt+eros',berg_melt_net, &
+                                              'berg mass',bergs_mass_start,bergs_mass_end)
+         CALL report_budget( 'bits mass','kg','eros used',bits_src_net, &
+                                              'bergs',bits_melt_net, &
+                                              'stored mass',bits_mass_start,bits_mass_end)
+         CALL report_budget( 'net mass','kg','recvd',calving_rcv_net, &
+                                             'rtrnd',calving_ret_net, &
+                                             'net mass',stored_start+floating_mass_start, &
+                                                        stored_end+floating_mass_end)
+         CALL report_consistant( 'iceberg mass','kg','gridded',zgrdd_berg_mass,'bergs',bergs_mass_end)
+         CALL report_consistant( 'bits mass','kg','gridded',zgrdd_bits_mass,'bits',bits_mass_end)
+         CALL report_state( 'net heat','J','',stored_heat_start+floating_heat_start,'', &
+                                              stored_heat_end+floating_heat_end,'')
+         CALL report_state( 'stored heat','J','',stored_heat_start,'',stored_heat_end,'')
+         CALL report_state( 'floating heat','J','',floating_heat_start,'',floating_heat_end,'')
+         CALL report_budget( 'net heat','J','net heat',calving_src_heat_net, &
+                                            'net heat',calving_out_heat_net, &
+                                            'net heat',stored_heat_start+floating_heat_start, &
+                                                       stored_heat_end+floating_heat_end)
+         CALL report_budget( 'stored heat','J','calving used',calving_src_heat_used_net, &
+                                               'bergs',heat_to_bergs_net, &
+                                               'net heat',stored_heat_start,stored_heat_end)
+         CALL report_budget( 'flting heat','J','calved',heat_to_bergs_net, &
+                                               'melt',heat_to_ocean_net, &
+                                               'net heat',floating_heat_start,floating_heat_end)
+         IF (nn_verbose_level >= 1) THEN
+            CALL report_consistant( 'top interface','kg','from SIS',calving_src_net, &
+                                    'received',calving_rcv_net)
+            CALL report_consistant( 'bot interface','kg','sent',calving_out_net, &
+                                    'returned',calving_ret_net)
+         ENDIF
+         WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses)
+         IF ( nspeeding_tickets > 0 ) WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets
+
+         nbergs_start              = nbergs_end
+         stored_start              = stored_end
+         nbergs_melted             = 0
+         nbergs_calved             = 0
+         nbergs_calved_by_class(:) = 0
+         nspeeding_tickets         = 0
+         stored_heat_start         = stored_heat_end
+         floating_heat_start       = floating_heat_end
+         floating_mass_start       = floating_mass_end
+         bergs_mass_start          = bergs_mass_end
+         bits_mass_start           = bits_mass_end
+         calving_used_net          = 0._wp
+         calving_to_bergs_net      = 0._wp
+         heat_to_bergs_net         = 0._wp
+         heat_to_ocean_net         = 0._wp
+         calving_rcv_net           = 0._wp
+         calving_ret_net           = 0._wp
+         calving_src_net           = 0._wp
+         calving_out_net           = 0._wp
+         calving_src_heat_net      = 0._wp
+         calving_src_heat_used_net = 0._wp
+         calving_out_heat_net      = 0._wp
+         melt_net                  = 0._wp
+         berg_melt_net             = 0._wp
+         bits_melt_net             = 0._wp
+         bits_src_net              = 0._wp
+      ENDIF
+      !
+   END SUBROUTINE icb_dia
+
+
+   SUBROUTINE icb_dia_step
+      !!----------------------------------------------------------------------
+      !! things to reset at the beginning of each timestep
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+      berg_melt    (:,:)   = 0._wp
+      buoy_melt    (:,:)   = 0._wp
+      eros_melt    (:,:)   = 0._wp
+      conv_melt    (:,:)   = 0._wp
+      bits_src     (:,:)   = 0._wp
+      bits_melt    (:,:)   = 0._wp
+      bits_mass    (:,:)   = 0._wp
+      berg_mass    (:,:)   = 0._wp
+      virtual_area (:,:)   = 0._wp
+      real_calving (:,:,:) = 0._wp
+      !
+   END SUBROUTINE icb_dia_step
+
+
+   SUBROUTINE icb_dia_put
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia )   RETURN            !!gm useless iom will control whether it is output or not
+      !
+      CALL iom_put( "berg_melt"        , berg_melt   (:,:)   )   ! Melt rate of icebergs                     [kg/m2/s]
+      CALL iom_put( "berg_buoy_melt"   , buoy_melt   (:,:)   )   ! Buoyancy component of iceberg melt rate   [kg/m2/s]
+      CALL iom_put( "berg_eros_melt"   , eros_melt   (:,:)   )   ! Erosion component of iceberg melt rate    [kg/m2/s]
+      CALL iom_put( "berg_conv_melt"   , conv_melt   (:,:)   )   ! Convective component of iceberg melt rate [kg/m2/s]
+      CALL iom_put( "berg_virtual_area", virtual_area(:,:)   )   ! Virtual coverage by icebergs              [m2]
+      CALL iom_put( "bits_src"         , bits_src    (:,:)   )   ! Mass source of bergy bits                 [kg/m2/s]
+      CALL iom_put( "bits_melt"        , bits_melt   (:,:)   )   ! Melt rate of bergy bits                   [kg/m2/s]
+      CALL iom_put( "bits_mass"        , bits_mass   (:,:)   )   ! Bergy bit density field                   [kg/m2]
+      CALL iom_put( "berg_mass"        , berg_mass   (:,:)   )   ! Iceberg density field                     [kg/m2]
+      CALL iom_put( "berg_real_calving", real_calving(:,:,:) )   ! Calving into iceberg class                [kg/s]
+      !
+   END SUBROUTINE icb_dia_put
+
+
+   SUBROUTINE icb_dia_calve( ki, kj, kn, pcalved, pheated )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER,  INTENT(in)  ::   ki, kj, kn
+      REAL(wp), INTENT(in)  ::   pcalved
+      REAL(wp), INTENT(in)  ::   pheated
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+      real_calving(ki,kj,kn)     = real_calving(ki,kj,kn) + pcalved / berg_dt
+      nbergs_calved              = nbergs_calved              + 1
+      nbergs_calved_by_class(kn) = nbergs_calved_by_class(kn) + 1
+      calving_to_bergs_net       = calving_to_bergs_net + pcalved
+      heat_to_bergs_net          = heat_to_bergs_net    + pheated
+      !
+   END SUBROUTINE icb_dia_calve
+
+
+   SUBROUTINE icb_dia_income( kt,  pcalving_used, pheat_used )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER ,                 INTENT(in)  :: kt
+      REAL(wp),                 INTENT(in)  :: pcalving_used
+      REAL(wp), DIMENSION(:,:), INTENT(in)  :: pheat_used
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+      !
+      IF( kt == nit000 ) THEN
+         stored_start = SUM( berg_grid%stored_ice(:,:,:) )
+         IF( lk_mpp ) CALL mpp_sum( stored_start )
+         WRITE(numicb,'(a,es13.6,a)')   'icb_dia_income: initial stored mass=',stored_start,' kg'
+         !
+         stored_heat_start = SUM( berg_grid%stored_heat(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( stored_heat_start )
+         WRITE(numicb,'(a,es13.6,a)')    'icb_dia_income: initial stored heat=',stored_heat_start,' J'
+      ENDIF
+      !
+      calving_rcv_net = calving_rcv_net + SUM( berg_grid%calving(:,:) ) * berg_dt
+      calving_src_net = calving_rcv_net
+      calving_src_heat_net = calving_src_heat_net +  &
+         &                      SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) ) * berg_dt   ! Units of J
+      calving_used_net = calving_used_net + pcalving_used * berg_dt
+      calving_src_heat_used_net = calving_src_heat_used_net + SUM( pheat_used(:,:) )
+      !
+   END SUBROUTINE icb_dia_income
+
+
+   SUBROUTINE icb_dia_size(ki, kj, pWn, pLn, pAbits,   &
+      &                    pmass_scale, pMnew, pnMbits, pz1_e1e2)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER,  INTENT(in) :: ki, kj
+      REAL(wp), INTENT(in) :: pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+      virtual_area(ki,kj) = virtual_area(ki,kj) + ( pWn * pLn + pAbits ) * pmass_scale      ! m^2
+      berg_mass(ki,kj)    = berg_mass(ki,kj) + pMnew * pz1_e1e2                             ! kg/m2
+      bits_mass(ki,kj)    = bits_mass(ki,kj) + pnMbits * pz1_e1e2                           ! kg/m2
+      !
+   END SUBROUTINE icb_dia_size
+
+
+   SUBROUTINE icb_dia_speed()
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+      nspeeding_tickets = nspeeding_tickets + 1
+      !
+   END SUBROUTINE icb_dia_speed
+
+
+   SUBROUTINE icb_dia_melt(ki, kj, pmnew, pheat, pmass_scale,   &
+      &                   pdM, pdMbitsE, pdMbitsM, pdMb, pdMe,   &
+      &                   pdMv, pz1_dt_e1e2 )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in) ::   ki, kj
+      REAL(wp), INTENT(in) ::   pmnew, pheat, pmass_scale
+      REAL(wp), INTENT(in) ::   pdM, pdMbitsE, pdMbitsM, pdMb, pdMe, pdMv, pz1_dt_e1e2
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_bergdia ) RETURN
+
+      berg_melt (ki,kj) = berg_melt (ki,kj) + pdM      * pz1_dt_e1e2   ! kg/m2/s
+      bits_src  (ki,kj) = bits_src  (ki,kj) + pdMbitsE * pz1_dt_e1e2   ! mass flux into bergy bitskg/m2/s
+      bits_melt (ki,kj) = bits_melt (ki,kj) + pdMbitsM * pz1_dt_e1e2   ! melt rate of bergy bits kg/m2/s
+      buoy_melt (ki,kj) = buoy_melt (ki,kj) + pdMb     * pz1_dt_e1e2   ! kg/m2/s
+      eros_melt (ki,kj) = eros_melt (ki,kj) + pdMe     * pz1_dt_e1e2   ! erosion rate kg/m2/s
+      conv_melt (ki,kj) = conv_melt (ki,kj) + pdMv     * pz1_dt_e1e2   ! kg/m2/s
+      heat_to_ocean_net = heat_to_ocean_net + pheat * pmass_scale * berg_dt         ! J
+      IF( pmnew <= 0._wp ) nbergs_melted = nbergs_melted + 1                        ! Delete the berg if completely melted
+      !
+   END SUBROUTINE icb_dia_melt
+
+
+   SUBROUTINE report_state( cd_budgetstr, cd_budgetunits, cd_startstr, pstartval, cd_endstr,   &
+      &                     pendval, cd_delstr, kbergs )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER*(*), INTENT(in)           :: cd_budgetstr, cd_budgetunits, cd_startstr, cd_endstr, cd_delstr
+      REAL(wp),      INTENT(in)           :: pstartval, pendval
+      INTEGER,       INTENT(in), OPTIONAL :: kbergs
+      !!----------------------------------------------------------------------
+      !
+      IF ( PRESENT(kbergs) ) THEN
+         WRITE(numicb,100) cd_budgetstr // ' state:',                                    &
+                           cd_startstr  // ' start',  pstartval,         cd_budgetunits, &
+                           cd_endstr    // ' end',    pendval,           cd_budgetunits, &
+                           'Delta '     // cd_delstr, pendval-pstartval, cd_budgetunits, &
+                           '# of bergs', kbergs
+      ELSE
+         WRITE(numicb,100) cd_budgetstr // ' state:',                                   &
+                           cd_startstr  // ' start', pstartval,         cd_budgetunits, &
+                           cd_endstr    // ' end',   pendval,           cd_budgetunits, &
+                           cd_delstr    // 'Delta',  pendval-pstartval, cd_budgetunits
+      ENDIF
+      100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8)
+   END SUBROUTINE report_state
+
+
+   SUBROUTINE report_consistant( cd_budgetstr, cd_budgetunits, cd_startstr, pstartval, cd_endstr, pendval)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_budgetunits, cd_startstr, cd_endstr
+      REAL(wp),      INTENT(in) :: pstartval, pendval
+      !!----------------------------------------------------------------------
+      !
+      WRITE(numicb,200) cd_budgetstr // ' check:',                 &
+                        cd_startstr,    pstartval, cd_budgetunits, &
+                        cd_endstr,      pendval,   cd_budgetunits, &
+                        'error',        (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd'
+      200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,","))
+   END SUBROUTINE report_consistant
+
+
+   SUBROUTINE report_budget( cd_budgetstr, cd_budgetunits, cd_instr, pinval, cd_outstr,   &
+      &                      poutval, cd_delstr, pstartval, pendval)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_budgetunits, cd_instr, cd_outstr, cd_delstr
+      REAL(wp),      INTENT(in) :: pinval, poutval, pstartval, pendval
+      !
+      REAL(wp)                  :: zval
+      !!----------------------------------------------------------------------
+      !
+      zval = ( ( pendval - pstartval ) - ( pinval - poutval ) ) /   &
+         &   MAX( 1.e-30, MAX( abs( pendval - pstartval ) , ABS( pinval - poutval ) ) )
+
+      WRITE(numicb,200) cd_budgetstr // ' budget:', &
+         &              cd_instr     // ' in',      pinval,         cd_budgetunits, &
+         &              cd_outstr    // ' out',     poutval,        cd_budgetunits, &
+         &              'Delta '     // cd_delstr,  pinval-poutval, cd_budgetunits, &
+         &              'error',        zval,                       'nd'
+  200 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a8,"=",es10.3,x,a2)
+      !
+   END SUBROUTINE report_budget
+
+
+   SUBROUTINE report_istate( cd_budgetstr, cd_startstr, pstartval, cd_endstr, pendval, cd_delstr)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_startstr, cd_endstr, cd_delstr
+      INTEGER,       INTENT(in) :: pstartval, pendval
+      !
+      WRITE(numicb,100) cd_budgetstr // ' state:',           &
+         &              cd_startstr  // ' start', pstartval, &
+         &              cd_endstr    // ' end',   pendval,   &
+         &              cd_delstr    // 'Delta',  pendval-pstartval
+  100 FORMAT(a19,3(a18,"=",i14,x,:,","))
+      !
+   END SUBROUTINE report_istate
+
+
+   SUBROUTINE report_ibudget( cd_budgetstr, cd_instr, pinval, cd_outstr, poutval,   &
+      &                       cd_delstr, pstartval, pendval)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_instr, cd_outstr, cd_delstr
+      INTEGER,       INTENT(in) :: pinval, poutval, pstartval, pendval
+      !!----------------------------------------------------------------------
+      !
+      WRITE(numicb,200) cd_budgetstr // ' budget:', &
+                        cd_instr     // ' in',      pinval, &
+                        cd_outstr    // ' out',     poutval, &
+                        'Delta '     // cd_delstr,  pinval-poutval, &
+                        'error',                    ( ( pendval - pstartval ) - ( pinval - poutval ) )
+      200 FORMAT(a19,10(a18,"=",i14,x,:,","))
+   END SUBROUTINE report_ibudget
+
+   !!======================================================================
+END MODULE icbdia
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 3609)
@@ -0,0 +1,380 @@
+MODULE icbdyn
+
+   !!======================================================================
+   !!                       ***  MODULE  icbdyn  ***
+   !! Iceberg:  time stepping routine for iceberg tracking
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !  2011-05  (Alderson)       Replace broken grounding routine
+   !!            -    !                            with one of Gurvan's suggestions (just like
+   !!            -    !                            the broken one)
+   !!----------------------------------------------------------------------
+   USE par_oce        ! NEMO parameters
+   USE dom_oce        ! NEMO ocean domain
+   USE phycst         ! NEMO physical constants
+   !
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+   USE icbdia         ! iceberg budget routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_dyn()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_dyn  ***
+      !!
+      !! ** Purpose :   iceberg evolution.
+      !!
+      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
+      !!----------------------------------------------------------------------
+      REAL(wp)                        ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
+      REAL(wp)                        ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
+      REAL(wp)                        ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
+      REAL(wp)                        ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
+      REAL(wp)                        ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
+      REAL(wp)                        ::   zdt, zdt_2, zdt_6, ze1, ze2
+      LOGICAL                         ::   ll_bounced
+      TYPE(iceberg), POINTER          ::   berg
+      TYPE(point)  , POINTER          ::   pt
+      !!----------------------------------------------------------------------
+
+      ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
+      !                    with I.C.'s:   X=X1 and V=V1
+      !
+      !                                    ; A1=A(X1,V1)
+      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
+      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
+      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
+      !
+      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
+      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
+
+      ! time steps
+      zdt   = berg_dt
+      zdt_2 = zdt * 0.5_wp
+      zdt_6 = zdt / 6._wp
+
+      berg => first_berg                    ! start from the first berg
+      !
+      DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+         !
+         pt => berg%current_point
+
+         ll_bounced = .false.
+
+
+         ! STEP 1 !
+         ! ====== !
+         zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
+         zyj1 = pt%yj   ;   zvvel1 = pt%vvel
+
+
+         !                                         !**   A1 = A(X1,V1)
+         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
+            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+         !
+         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
+         zv1 = zvvel1 / ze2
+
+         ! STEP 2 !
+         ! ====== !
+         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
+         ! position using di/dt & djdt   !   V2  in m/s
+         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
+         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
+         !
+         CALL icb_ground( zxi2, zxi1, zu1,   &
+         &                zyj2, zyj1, zv1, ll_bounced )
+
+         !                                         !**   A2 = A(X2,V2)
+         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
+            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+         !
+         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
+         zv2 = zvvel2 / ze2
+         !
+         ! STEP 3 !
+         ! ====== !
+         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
+         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
+         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
+         !
+         CALL icb_ground( zxi3, zxi1, zu3,   &
+         &                zyj3, zyj1, zv3, ll_bounced )
+
+         !                                         !**   A3 = A(X3,V3)
+         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
+            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+         !
+         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
+         zv3 = zvvel3 / ze2
+
+         ! STEP 4 !
+         ! ====== !
+         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
+         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
+         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
+
+         CALL icb_ground( zxi4, zxi1, zu4,   &
+         &                zyj4, zyj1, zv4, ll_bounced )
+
+         !                                         !**   A4 = A(X4,V4)
+         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
+            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+
+         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
+         zv4 = zvvel4 / ze2
+
+         ! FINAL STEP !
+         ! ========== !
+         !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
+         !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
+         zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
+         zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
+         zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 )
+         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
+
+         CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
+         &                      zyj_n, zyj1, zvvel_n, ll_bounced )
+
+         pt%uvel = zuvel_n                        !** save in berg structure
+         pt%vvel = zvvel_n
+         pt%xi   = zxi_n
+         pt%yj   = zyj_n
+
+         ! update actual position
+         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
+         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
+
+         berg => berg%next                         ! switch to the next berg
+         !
+      END DO                                  !==  end loop over all bergs  ==!
+      !
+   END SUBROUTINE icb_dyn
+
+
+   SUBROUTINE icb_ground( pi, pi0, pu,   &
+      &                         pj, pj0, pv, ld_bounced )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_ground  ***
+      !!
+      !! ** Purpose :   iceberg grounding.
+      !!
+      !! ** Method  : - adjust velocity and then put iceberg back to start position
+      !!                NB two possibilities available one of which is hard-coded here
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
+      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
+      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
+      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
+      !
+      INTEGER  ::   ii, ii0
+      INTEGER  ::   ij, ij0
+      INTEGER  ::   ibounce_method
+      !!----------------------------------------------------------------------
+      !
+      ld_bounced = .FALSE.
+      !
+      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
+      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
+      !
+      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
+      !
+      ! map into current processor
+      ii0 = mi1( ii0 )
+      ij0 = mj1( ij0 )
+      ii  = mi1( ii  )
+      ij  = mj1( ij  )
+      !
+      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      !
+      ! From here, berg have reach land: treat grounding/bouncing
+      ! -------------------------------
+      ld_bounced = .TRUE.
+
+      !! not obvious what should happen now
+      !! if berg tries to enter a land box, the only location we can return it to is the start 
+      !! position (pi0,pj0), since it has to be in a wet box to do any melting;
+      !! first option is simply to set whole velocity to zero and move back to start point
+      !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
+      !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast
+
+      ibounce_method = 2
+      SELECT CASE ( ibounce_method )
+         CASE ( 1 )
+            pi = pi0
+            pj = pj0
+            pu = 0._wp
+            pv = 0._wp
+         CASE ( 2 )
+            IF( ii0 /= ii ) THEN
+               pi = pi0                   ! return back to the initial position
+               pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
+            ENDIF
+            IF( ij0 /= ij ) THEN
+               pj = pj0                   ! return back to the initial position
+               pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
+            ENDIF
+      END SELECT
+      !
+   END SUBROUTINE icb_ground
+
+
+   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
+      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_accel  ***
+      !!
+      !! ** Purpose :   compute the iceberg acceleration.
+      !!
+      !! ** Method  : - sum the terms in the momentum budget
+      !!----------------------------------------------------------------------
+      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
+      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
+      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
+      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
+      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
+      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
+      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
+      !
+      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
+      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
+      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
+      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
+      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
+      !
+      INTEGER  ::   itloop
+      REAL(wp) ::   zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi
+      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
+      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
+      REAL(wp) ::   z_ocn, z_atm, z_ice
+      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
+      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
+      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
+      !!----------------------------------------------------------------------
+
+      ! Interpolate gridded fields to berg
+      nknberg = berg%number(1)
+      CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     &
+         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff )
+
+      zM = berg%current_point%mass
+      zT = berg%current_point%thickness               ! total thickness
+      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
+      zF = zT - zD                                    ! freeboard
+      zW = berg%current_point%width
+      zL = berg%current_point%length
+
+      zhi   = MIN( zhi   , zD    )
+      zD_hi = MAX( 0._wp, zD-zhi )
+
+      ! Wave radiation
+      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
+      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
+      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
+      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
+      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
+      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
+      zLcutoff     = 0.125 * zLwavelength
+      zLtop        = 0.25  * zLwavelength
+      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
+      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
+      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
+      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
+      IF( zwmod /= 0._wp ) THEN
+         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
+         zvwave = zva/zwmod
+      ELSE
+         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
+      ENDIF
+
+      ! Weighted drag coefficients
+      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
+      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
+      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
+      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
+
+      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
+      !
+      DO itloop = 1, 2  ! Iterate on drag coefficients
+         !
+         zus = 0.5 * ( zuveln + puvel )
+         zvs = 0.5 * ( zvveln + pvvel )
+         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
+         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
+         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
+         !
+         ! Explicit accelerations
+         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
+         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
+         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
+         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
+         zaxe = -grav * zssh_x + zwave_rad * zuwave
+         zaye = -grav * zssh_y + zwave_rad * zvwave
+         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
+            zaxe = zaxe + zff*pvvel0
+            zaye = zaye - zff*puvel0
+         ELSE
+            zaxe = zaxe + zff*pvvel
+            zaye = zaye - zff*puvel
+         ENDIF
+         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
+            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
+            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
+         ELSE
+            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
+            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
+         endif
+
+         ! Solve for implicit accelerations
+         IF( pp_alpha + pp_beta > 0._wp ) THEN
+            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
+            zA11    = 1._wp + pp_beta *pdt*zlambda
+            zA12    =         pp_alpha*pdt*zff
+            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
+            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
+            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
+         else
+            pax = zaxe   ;   pay = zaye
+         endif
+
+         zuveln = puvel0 + pdt*pax
+         zvveln = pvvel0 + pdt*pay
+         !
+      END DO      ! itloop
+
+      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
+         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
+         IF( zspeed > 0._wp ) THEN
+            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
+            zspeed_new = zloc_dx / pdt * rn_speed_limit     ! Speed limit as a factor of dx / dt
+            IF( zspeed_new < zspeed ) THEN
+               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
+               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
+               CALL icb_dia_speed()
+            ENDIF
+         ENDIF
+      ENDIF
+      !                                      ! check the speed and acceleration limits
+      IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
+         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
+      IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
+         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
+      !
+   END SUBROUTINE icb_accel
+
+   !!======================================================================
+END MODULE icbdyn
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 3609)
@@ -0,0 +1,433 @@
+MODULE icbini
+
+   !!======================================================================
+   !!                       ***  MODULE  icbini  ***
+   !! Icebergs:  initialise variables for iceberg tracking
+   !!======================================================================
+   !! History :   -   !  2010-01  (T. Martin & A. Adcroft)  Original code
+   !!            3.3  !  2011-03  (G. Madec)  Part conversion to NEMO form ; Removal of mapping from another grid
+   !!             -   !  2011-04  (S. Alderson)  Split into separate modules ; Restore restart routines
+   !!             -   !  2011-05  (S. Alderson)  generate_test_icebergs restored ; new forcing arrays with extra halo ;
+   !!             -   !                          north fold exchange arrays added
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_init     : initialise icebergs
+   !!   icb_ini_gen  : generate test icebergs
+   !!   icb_nam      : read iceberg namelist
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean domain
+   USE in_out_manager ! IO routines and numout in particular
+   USE lib_mpp        ! mpi library and lk_mpp in particular
+   USE sbc_oce        ! ocean: surface boundary condition
+   USE iom            ! IOM library
+   USE fldread        ! field read
+   USE lbclnk         ! lateral boundary condition - MPP link
+   !
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+   USE icbrst         ! iceberg restart routines
+   USE icbtrj         ! iceberg trajectory I/O routines
+   USE icbdia         ! iceberg budget routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_init  ! routine called in nemogcm.F90 module
+
+   CHARACTER(len=100) ::   cn_dir = './'   ! Root directory for location of icb files
+   TYPE(FLD_N)        ::   sn_icb          ! information about the calving file to be read
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_init( pdt, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_init  ***
+      !!
+      !! ** Purpose :   iceberg initialization.
+      !!
+      !! ** Method  : - read the iceberg namelist
+      !!              - find non-overlapping processor interior since we can only
+      !!                have one instance of a particular iceberg
+      !!              - calculate the destinations for north fold exchanges
+      !!              - setup either test icebergs or calving file
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pdt   ! iceberg time-step (rdt*nn_fsbc)
+      INTEGER , INTENT(in) ::   kt    ! time step number
+      !
+      INTEGER ::   ji, jj, jn               ! dummy loop indices
+      INTEGER ::   i1, i2, i3               ! local integers
+      INTEGER ::   ii, inum, ivar           !   -       -
+      INTEGER ::   istat1, istat2, istat3   !   -       -
+      CHARACTER(len=300) ::   cl_sdist      ! local character
+      !!----------------------------------------------------------------------
+      !
+      CALL icb_nam               ! Read and print namelist parameters
+      !
+      IF( .NOT. ln_icebergs )   RETURN
+
+      !                          ! allocate gridded fields
+      IF( icb_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' )
+
+      !                          ! open ascii output file or files for iceberg status information
+      !                          ! note that we choose to do this on all processors since we cannot
+      !                          ! predict where icebergs will be ahead of time
+      CALL ctl_opn( numicb, 'icebergs.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+
+      ! set parameters (mostly from namelist)
+      !
+      berg_dt         = pdt
+      first_width (:) = SQRT(  rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) )  )
+      first_length(:) = rn_LoW_ratio * first_width(:)
+
+      berg_grid%calving      (:,:)   = 0._wp
+      berg_grid%calving_hflx (:,:)   = 0._wp
+      berg_grid%stored_heat  (:,:)   = 0._wp
+      berg_grid%floating_melt(:,:)   = 0._wp
+      berg_grid%maxclass     (:,:)   = nclasses
+      berg_grid%stored_ice   (:,:,:) = 0._wp
+      berg_grid%tmp          (:,:)   = 0._wp
+      src_calving            (:,:)   = 0._wp
+      src_calving_hflx       (:,:)   = 0._wp
+
+      !                          ! domain for icebergs
+      IF( lk_mpp .AND. jpni == 1 )   CALL ctl_stop( 'icbinit: having ONE processor in x currently does not work' )
+      ! NB: the issue here is simply that cyclic east-west boundary condition have not been coded in mpp case
+      ! for the north fold we work out which points communicate by asking
+      ! lbc_lnk to pass processor number (valid even in single processor case)
+      ! borrow src_calving arrays for this
+      !
+      ! pack i and j together using a scaling of a power of 10
+      nicbpack = 10000
+      IF( jpiglo >= nicbpack )   CALL ctl_stop( 'icbini: processor index packing failure' )
+      nicbfldproc(:) = -1
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            src_calving_hflx(ji,jj) = narea
+            src_calving     (ji,jj) = nicbpack * mjg(jj) + mig(ji)
+         END DO
+      END DO
+      CALL lbc_lnk( src_calving_hflx, 'T', 1._wp )
+      CALL lbc_lnk( src_calving     , 'T', 1._wp )
+
+      ! work out interior of processor from exchange array
+      ! first entry with narea for this processor is left hand interior index
+      ! last  entry                               is right hand interior index
+      jj = jpj/2
+      nicbdi = -1
+      nicbei = -1
+      DO ji = 1, jpi
+         i3 = INT( src_calving(ji,jj) )
+         i2 = INT( i3/nicbpack )
+         i1 = i3 - i2*nicbpack
+         i3 = INT( src_calving_hflx(ji,jj) )
+         IF( i1 == mig(ji) .AND. i3 == narea ) THEN
+            IF( nicbdi < 0 ) THEN   ;   nicbdi = ji
+            ELSE                    ;   nicbei = ji
+            ENDIF
+         ENDIF
+      END DO
+      !
+      ! repeat for j direction
+      ji = jpi/2
+      nicbdj = -1
+      nicbej = -1
+      DO jj = 1, jpj
+         i3 = INT( src_calving(ji,jj) )
+         i2 = INT( i3/nicbpack )
+         i1 = i3 - i2*nicbpack
+         i3 = INT( src_calving_hflx(ji,jj) )
+         IF( i2 == mjg(jj) .AND. i3 == narea ) THEN
+            IF( nicbdj < 0 ) THEN   ;   nicbdj = jj
+            ELSE                    ;   nicbej = jj
+            ENDIF
+         ENDIF
+      END DO
+      !   
+      ! special for east-west boundary exchange we save the destination index
+      i1 = MAX( nicbdi-1, 1)
+      i3 = INT( src_calving(i1,jpj/2) )
+      jj = INT( i3/nicbpack )
+      ricb_left = REAL( i3 - nicbpack*jj, wp )
+      i1 = MIN( nicbei+1, jpi )
+      i3 = INT( src_calving(i1,jpj/2) )
+      jj = INT( i3/nicbpack )
+      ricb_right = REAL( i3 - nicbpack*jj, wp )
+      
+      ! north fold
+      IF( npolj > 0 ) THEN
+         !
+         ! icebergs in row nicbej+1 get passed across fold
+         nicbfldpts(:)  = INT( src_calving(:,nicbej+1) )
+         nicbflddest(:) = INT( src_calving_hflx(:,nicbej+1) )
+         !
+         ! work out list of unique processors to talk to
+         ! pack them into a fixed size array where empty slots are marked by a -1
+         DO ji = nicbdi, nicbei
+            ii = nicbflddest(ji)
+            DO jn = 1, jpni
+               ! work along array until we find an empty slot
+               IF( nicbfldproc(jn) == -1 ) THEN
+                  nicbfldproc(jn) = ii
+                  EXIT                             !!gm EXIT should be avoided: use DO WHILE expression instead
+               ENDIF
+               ! before we find an empty slot, we may find processor number is already here so we exit
+               IF( nicbfldproc(jn) == ii ) EXIT
+            END DO
+         END DO
+      ENDIF
+      !
+      IF( nn_verbose_level > 0) THEN
+         WRITE(numicb,*) 'processor ', narea
+         WRITE(numicb,*) 'jpi, jpj   ', jpi, jpj
+         WRITE(numicb,*) 'nldi, nlei ', nldi, nlei
+         WRITE(numicb,*) 'nldj, nlej ', nldj, nlej
+         WRITE(numicb,*) 'berg i interior ', nicbdi, nicbei
+         WRITE(numicb,*) 'berg j interior ', nicbdj, nicbej
+         WRITE(numicb,*) 'berg left       ', ricb_left
+         WRITE(numicb,*) 'berg right      ', ricb_right
+         jj = jpj/2
+         WRITE(numicb,*) "central j line:"
+         WRITE(numicb,*) "i processor"
+         WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), ji=1,jpi)
+         WRITE(numicb,*) "i point"
+         WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi)
+         ji = jpi/2
+         WRITE(numicb,*) "central i line:"
+         WRITE(numicb,*) "j processor"
+         WRITE(numicb,*) (INT(src_calving_hflx(ji,jj)), jj=1,jpj)
+         WRITE(numicb,*) "j point"
+         WRITE(numicb,*) (INT(src_calving(ji,jj)), jj=1,jpj)
+         IF( npolj > 0 ) THEN
+            WRITE(numicb,*) 'north fold destination points '
+            WRITE(numicb,*) nicbfldpts
+            WRITE(numicb,*) 'north fold destination procs  '
+            WRITE(numicb,*) nicbflddest
+         ENDIF
+         CALL flush(numicb)
+      ENDIF
+      
+      src_calving     (:,:) = 0._wp
+      src_calving_hflx(:,:) = 0._wp
+
+      ! assign each new iceberg with a unique number constructed from the processor number
+      ! and incremented by the total number of processors
+      num_bergs(:) = 0
+      num_bergs(1) = narea - jpnij
+
+      ! when not generating test icebergs we need to setup calving file
+      IF( nn_test_icebergs < 0 ) THEN
+         !
+         ! maximum distribution class array does not change in time so read it once
+         cl_sdist = TRIM( cn_dir )//TRIM( sn_icb%clname )
+         CALL iom_open ( cl_sdist, inum )                              ! open file
+         ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. )
+         IF( ivar > 0 ) THEN
+            CALL iom_get  ( inum, jpdom_data, 'maxclass', src_calving )   ! read the max distribution array
+            berg_grid%maxclass(:,:) = INT( src_calving )
+            src_calving(:,:) = 0._wp
+         ENDIF
+         CALL iom_close( inum )                                     ! close file
+         !
+         WRITE(numicb,*)
+         WRITE(numicb,*) '          calving read in a file'
+         ALLOCATE( sf_icb(1), STAT=istat1 )         ! Create sf_icb structure (calving)
+         ALLOCATE( sf_icb(1)%fnow(jpi,jpj,1), STAT=istat2 )
+         ALLOCATE( sf_icb(1)%fdta(jpi,jpj,1,2), STAT=istat3 )
+         IF( istat1+istat2+istat3 > 0 ) THEN
+            CALL ctl_stop( 'sbc_icb: unable to allocate sf_icb structure' )   ;   RETURN
+         ENDIF
+         !                                          ! fill sf_icb with the namelist (sn_icb) and control print
+         CALL fld_fill( sf_icb, (/ sn_icb /), cn_dir, 'icb_init', 'read calving data', 'namicb' )
+         !
+      ENDIF
+
+      IF( .NOT.ln_rstart ) THEN
+         IF( nn_test_icebergs > 0 )   CALL icb_ini_gen()
+      ELSE
+         IF( nn_test_icebergs > 0 ) THEN
+            CALL icb_ini_gen()
+         ELSE
+            CALL icb_rst_read()
+            l_restarted_bergs = .TRUE.
+         ENDIF
+      ENDIF
+      !
+      IF( nn_sample_rate .GT. 0 ) CALL icb_trj_init( nitend )
+      !
+      CALL icb_dia_init()
+      !
+      IF( nn_verbose_level >= 2 )   CALL icb_utl_print('icb_init, initial status', nit000-1)
+      !
+   END SUBROUTINE icb_init
+
+   SUBROUTINE icb_ini_gen()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_ini_gen  ***
+      !!
+      !! ** Purpose :   iceberg generation
+      !!
+      !! ** Method  : - at each grid point of the test box supplied in the namelist
+      !!                generate an iceberg in one class determined by the value of
+      !!                parameter nn_test_icebergs
+      !!----------------------------------------------------------------------
+      INTEGER                         ::   ji, jj, ibergs
+      TYPE(iceberg)                   ::   localberg ! NOT a pointer but an actual local variable
+      TYPE(point)                     ::   localpt
+      INTEGER                         ::   iyr, imon, iday, ihr, imin, isec
+      INTEGER                         ::   iberg
+      !!----------------------------------------------------------------------
+
+      ! For convenience
+      iberg = nn_test_icebergs
+
+      ! call get_date(Time, iyr, imon, iday, ihr, imin, isec)
+      ! Convert nemo time variables from dom_oce into local versions
+      iyr  = nyear
+      imon = nmonth
+      iday = nday
+      ihr = INT(nsec_day/3600)
+      imin = INT((nsec_day-ihr*3600)/60)
+      isec = nsec_day - ihr*3600 - imin*60
+
+      ! no overlap for icebergs since we want only one instance of each across the whole domain
+      ! so restrict area of interest
+      ! use tmask here because tmask_i has been doctored on one side of the north fold line
+
+      DO jj = nicbdj, nicbej
+         DO ji = nicbdi, nicbei
+            IF( tmask(ji,jj,1) > 0._wp        .AND.                                       &
+                rn_test_box(1) < glamt(ji,jj) .AND. glamt(ji,jj) < rn_test_box(2) .AND.   &
+                rn_test_box(3) < gphit(ji,jj) .AND. gphit(ji,jj) < rn_test_box(4) ) THEN
+               localberg%mass_scaling = rn_mass_scaling(iberg)
+               localpt%xi = REAL( mig(ji), wp )
+               localpt%yj = REAL( mjg(jj), wp )
+               localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T' )
+               localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T' )
+               localpt%mass      = rn_initial_mass     (iberg)
+               localpt%thickness = rn_initial_thickness(iberg)
+               localpt%width  = first_width (iberg)
+               localpt%length = first_length(iberg)
+               localpt%year = iyr
+               localpt%day = REAL(iday,wp)+(REAL(ihr,wp)+REAL(imin,wp)/60._wp)/24._wp
+               localpt%mass_of_bits = 0._wp
+               localpt%heat_density = 0._wp
+               localpt%uvel = 0._wp
+               localpt%vvel = 0._wp
+               CALL icb_utl_incr()
+               localberg%number(:) = num_bergs(:)
+               call icb_utl_add(localberg, localpt)
+            ENDIF
+         END DO
+      END DO
+      !
+      ibergs = icb_utl_count()
+      IF( lk_mpp ) CALL mpp_sum(ibergs)
+      WRITE(numicb,'(a,i6,a)') 'diamonds, icb_ini_gen: ',ibergs,' were generated'
+      !
+   END SUBROUTINE icb_ini_gen
+
+   SUBROUTINE icb_nam
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE icb_nam  ***
+      !!
+      !! ** Purpose :   read iceberg namelist and print the variables.
+      !!
+      !! ** input   : - namberg namelist
+      !!----------------------------------------------------------------------
+      INTEGER  ::   jn      ! dummy loop indices
+      REAL(wp) ::   zfact   ! local scalar
+      !
+      NAMELIST/namberg/ ln_icebergs    , ln_bergdia     , nn_sample_rate      , rn_initial_mass      ,   &
+         &              rn_distribution, rn_mass_scaling, rn_initial_thickness, nn_verbose_write     ,   &
+         &              rn_rho_bergs   , rn_LoW_ratio   , nn_verbose_level    , ln_operator_splitting,   &
+         &              rn_bits_erosion_fraction        , rn_sicn_shift       , ln_passive_mode      ,   &
+         &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   &
+         &              rn_speed_limit , cn_dir, sn_icb
+      !!----------------------------------------------------------------------
+
+      ! (NB: frequency positive => hours, negative => months)
+      !            !   file     ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   !
+      !            !   name     !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      !
+      sn_icb = FLD_N( 'calving' ,    -1     , 'calving'  ,  .TRUE.    , .TRUE. ,   'yearly'  , ''       , ''         )
+
+      REWIND( numnam )              ! Namelist namrun : iceberg parameters
+      READ  ( numnam, namberg )
+      
+      IF( .NOT. ln_icebergs ) THEN   ! no icebergs
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'icb_nam :   ln_icebergs = F , NO icebergs used'
+            WRITE(numout,*) '~~~~~~~~ '
+         ENDIF
+         RETURN
+      ENDIF
+
+      IF( nn_test_icebergs > nclasses ) THEN
+          IF(lwp) WRITE(numout,*) 'Resetting nn_test_icebergs to ', nclasses
+          nn_test_icebergs = nclasses
+      ENDIF
+
+      zfact = SUM( rn_distribution )
+      IF( zfact < 1._wp ) THEN
+         IF( zfact <= 0._wp ) THEN
+            CALL ctl_stop( 'icb_nam: sum of berg distribution equal to zero' )
+         ELSE
+            rn_distribution(:) = rn_distribution(:) / zfact
+            CALL ctl_warn( 'icb_nam: sum of berg input distribution not equal to one and so RESCALED' )
+         ENDIF
+      ENDIF
+
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
+         WRITE(numout,*) '~~~~~~~~ '
+         WRITE(numout,*) '   Calculate budgets                                            ln_bergdia       = ', ln_bergdia
+         WRITE(numout,*) '   Period between sampling of position for trajectory storage   nn_sample_rate = ', nn_sample_rate
+         WRITE(numout,*) '   Mass thresholds between iceberg classes (kg)                 rn_initial_mass     ='
+         DO jn=1,nclasses
+            WRITE(numout,'(a,f15.2)') '                                                                ',rn_initial_mass(jn)
+         ENDDO
+         WRITE(numout,*) '   Fraction of calving to apply to this class (non-dim)         rn_distribution     ='
+         DO jn = 1, nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_distribution(jn)
+         END DO
+         WRITE(numout,*) '   Ratio between effective and real iceberg mass (non-dim)      rn_mass_scaling     = '
+         DO jn = 1, nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_mass_scaling(jn)
+         END DO
+         WRITE(numout,*) '   Total thickness of newly calved bergs (m)                    rn_initial_thickness = '
+         DO jn = 1, nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_initial_thickness(jn)
+         END DO
+         WRITE(numout,*) '   Timesteps between verbose messages                           nn_verbose_write    = ', nn_verbose_write
+
+         WRITE(numout,*) '   Density of icebergs                           rn_rho_bergs  = ', rn_rho_bergs
+         WRITE(numout,*) '   Initial ratio L/W for newly calved icebergs   rn_LoW_ratio  = ', rn_LoW_ratio
+         WRITE(numout,*) '   Turn on more verbose output                          level  = ', nn_verbose_level
+         WRITE(numout,*) '   Use first order operator splitting for thermodynamics    ',   &
+            &                    'use_operator_splitting = ', ln_operator_splitting
+         WRITE(numout,*) '   Fraction of erosion melt flux to divert to bergy bits    ',   &
+            &                    'bits_erosion_fraction = ', rn_bits_erosion_fraction
+
+         WRITE(numout,*) '   Shift of sea-ice concentration in erosion flux modulation ',   &
+            &                    '(0<sicn_shift<1)    rn_sicn_shift  = ', rn_sicn_shift
+         WRITE(numout,*) '   Do not add freshwater flux from icebergs to ocean                ',   &
+            &                    '                  passive_mode            = ', ln_passive_mode
+         WRITE(numout,*) '   Time average the weight on the ocean   time_average_weight       = ', ln_time_average_weight
+         WRITE(numout,*) '   Create icebergs in absence of a restart file   nn_test_icebergs  = ', nn_test_icebergs
+         WRITE(numout,*) '                   in lon/lat box                                   = ', rn_test_box
+         WRITE(numout,*) '   CFL speed limit for a berg            speed_limit                = ', rn_speed_limit
+         WRITE(numout,*) '   Writing Iceberg status information to icebergs.stat file        '
+      ENDIF
+      !
+   END SUBROUTINE icb_nam
+
+   !!======================================================================
+
+END MODULE icbini
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 3609)
@@ -0,0 +1,826 @@
+MODULE icblbc
+
+   !!======================================================================
+   !!                       ***  MODULE  icblbc  ***
+   !! Ocean physics:  routines to handle boundary exchanges for icebergs
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !  2011-05  (Alderson)       MPP exchanges written based on lib_mpp
+   !!            -    !  2011-05  (Alderson)       MPP and single processor boundary
+   !!            -    !                            conditions added
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_lbc       : -  Pass icebergs across cyclic boundaries
+   !!   icb_lbc_mpp   : -  In MPP pass icebergs from linked list between processors
+   !!                      as they advect around
+   !!                   -  Lagrangian processes cannot be handled by existing NEMO MPP
+   !!                      routines because they do not lie on regular jpi,jpj grids
+   !!                   -  Processor exchanges are handled as in lib_mpp whenever icebergs step 
+   !!                      across boundary of interior domain (nicbdi-nicbei, nicbdj-nicbej)
+   !!                      so that iceberg does not exist in more than one processor
+   !!                   -  North fold exchanges controlled by three arrays:
+   !!                         nicbflddest - unique processor numbers that current one exchanges with
+   !!                         nicbfldproc - processor number that current grid point exchanges with
+   !!                         nicbfldpts  - packed i,j point in exchanging processor
+   !!----------------------------------------------------------------------
+
+   USE par_oce                             ! ocean parameters
+   USE dom_oce                             ! ocean domain
+   USE in_out_manager                      ! IO parameters
+   USE lib_mpp                             ! MPI code and lk_mpp in particular
+   USE icb_oce                             ! define iceberg arrays
+   USE icbutl                              ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+#if defined key_mpp_mpi
+
+!$AGRIF_DO_NOT_TREAT
+   INCLUDE 'mpif.h'
+!$AGRIF_END_DO_NOT_TREAT
+
+   TYPE, PUBLIC :: buffer
+      INTEGER :: size=0
+      REAL(wp), DIMENSION(:,:), POINTER :: data
+   END TYPE buffer
+
+   TYPE(buffer), POINTER           ::   obuffer_n=>NULL() , ibuffer_n=>NULL()
+   TYPE(buffer), POINTER           ::   obuffer_s=>NULL() , ibuffer_s=>NULL()
+   TYPE(buffer), POINTER           ::   obuffer_e=>NULL() , ibuffer_e=>NULL()
+   TYPE(buffer), POINTER           ::   obuffer_w=>NULL() , ibuffer_w=>NULL()
+
+   ! north fold exchange buffers
+   TYPE(buffer), POINTER           ::   obuffer_f=>NULL() , ibuffer_f=>NULL()
+
+   INTEGER, PARAMETER, PRIVATE     ::   jp_delta_buf = 25             ! Size by which to increment buffers
+   INTEGER, PARAMETER, PRIVATE     ::   jp_buffer_width = 15+nkounts  ! items to store for each berg
+
+#endif
+
+   PUBLIC   icb_lbc
+   PUBLIC   icb_lbc_mpp
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_lbc()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_lbc  ***
+      !!
+      !! ** Purpose :   in non-mpp case need to deal with cyclic conditions
+      !!                including north-fold
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER ::   this
+      TYPE(point)  , POINTER ::   pt
+      INTEGER                ::   iine
+      !!----------------------------------------------------------------------
+
+      !! periodic east/west boundaries
+      !! =============================
+
+      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
+
+         this => first_berg
+         DO WHILE( ASSOCIATED(this) )
+            pt => this%current_point
+            iine = INT( pt%xi + 0.5 )
+            IF( iine > mig(nicbei) ) THEN
+               pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp
+            ELSE IF( iine < mig(nicbdi) ) THEN
+               pt%xi = ricb_left + MOD(pt%xi, 1._wp )
+            ENDIF
+            this => this%next
+         END DO
+         !
+      ENDIF
+
+      !! north/south boundaries
+      !! ======================
+      ! south symmetric
+      IF( nperio == 2 )   CALL ctl_stop(' south symmetric condition not implemented for icebergs')
+      ! north fold
+      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   CALL icb_lbc_nfld()
+      !
+   END SUBROUTINE icb_lbc
+
+
+   SUBROUTINE icb_lbc_nfld()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_lbc_nfld  ***
+      !!
+      !! ** Purpose :   single processor north fold exchange
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER ::   this
+      TYPE(point)  , POINTER ::   pt
+      INTEGER                ::   iine, ijne, ipts
+      INTEGER                ::   iiglo, ijglo
+      !!----------------------------------------------------------------------
+      !
+      this => first_berg
+      DO WHILE( ASSOCIATED(this) )
+         pt => this%current_point
+         ijne = INT( pt%yj + 0.5 )
+         IF( ijne .GT. mjg(nicbej) ) THEN
+            !
+            iine = INT( pt%xi + 0.5 )
+            ipts  = nicbfldpts (mi1(iine))
+            !
+            ! moving across the cut line means both position and
+            ! velocity must change
+            ijglo = INT( ipts/nicbpack )
+            iiglo = ipts - nicbpack*ijglo
+            pt%xi = iiglo - ( pt%xi - REAL(iine,wp) )
+            pt%yj = ijglo - ( pt%yj - REAL(ijne,wp) )
+            pt%uvel = -1._wp * pt%uvel
+            pt%vvel = -1._wp * pt%vvel
+         ENDIF
+         this => this%next
+      END DO
+      !
+   END SUBROUTINE icb_lbc_nfld
+
+#if defined key_mpp_mpi
+   !!----------------------------------------------------------------------
+   !!   'key_mpp_mpi'             MPI massively parallel processing library
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE icb_lbc_mpp()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_lbc_mpp  ***
+      !!
+      !! ** Purpose :   multi processor exchange
+      !!
+      !! ** Method  :   identify direction for exchange, pack into a buffer
+      !!                which is basically a real array and delete from linked list
+      !!                length of buffer is exchanged first with receiving processor
+      !!                then buffer is sent if necessary
+      !!----------------------------------------------------------------------
+      TYPE(iceberg)         , POINTER     ::   tmpberg, this
+      TYPE(point)           , POINTER     ::   pt
+      INTEGER                             ::   ibergs_to_send_e, ibergs_to_send_w
+      INTEGER                             ::   ibergs_to_send_n, ibergs_to_send_s
+      INTEGER                             ::   ibergs_rcvd_from_e, ibergs_rcvd_from_w
+      INTEGER                             ::   ibergs_rcvd_from_n, ibergs_rcvd_from_s
+      INTEGER                             ::   i, ibergs_start, ibergs_end
+      INTEGER                             ::   iine, ijne
+      INTEGER                             ::   ipe_N, ipe_S, ipe_W, ipe_E
+      REAL(wp), DIMENSION(2)              ::   zewbergs, zwebergs, znsbergs, zsnbergs
+      INTEGER                             ::   iml_req1, iml_req2, iml_req3, iml_req4
+      INTEGER                             ::   iml_req5, iml_req6, iml_req7, iml_req8, iml_err
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   iml_stat
+
+      ! set up indices of neighbouring processors
+      ipe_N = -1
+      ipe_S = -1
+      ipe_W = -1
+      ipe_E = -1
+      IF( nbondi .EQ.  0 .OR. nbondi .EQ. 1) ipe_W = nowe
+      IF( nbondi .EQ. -1 .OR. nbondi .EQ. 0) ipe_E = noea
+      IF( nbondj .EQ.  0 .OR. nbondj .EQ. 1) ipe_S = noso
+      IF( nbondj .EQ. -1 .OR. nbondj .EQ. 0) ipe_N = nono
+      !
+      ! at northern line of processors with north fold handle bergs differently
+      IF( npolj > 0 ) ipe_N = -1
+
+      ! if there's only one processor in x direction then don't let mpp try to handle periodicity
+      IF( jpni == 1 ) THEN
+         ipe_E = -1
+         ipe_W = -1
+      ENDIF
+
+      IF( nn_verbose_level >= 2 ) THEN
+         WRITE(numicb,*) 'processor west  : ', ipe_W
+         WRITE(numicb,*) 'processor east  : ', ipe_E
+         WRITE(numicb,*) 'processor north : ', ipe_N
+         WRITE(numicb,*) 'processor south : ', ipe_S
+         WRITE(numicb,*) 'processor nimpp : ', nimpp
+         WRITE(numicb,*) 'processor njmpp : ', njmpp
+         WRITE(numicb,*) 'processor nbondi: ', nbondi
+         WRITE(numicb,*) 'processor nbondj: ', nbondj
+         CALL flush( numicb )
+      ENDIF
+
+      ! periodicity is handled here when using mpp when there is more than one processor in
+      ! the i direction, but it also has to happen when jpni=1 case so this is dealt with
+      ! in icb_lbc and called here
+
+      IF( jpni == 1 ) CALL icb_lbc()
+
+      ! Note that xi is adjusted when swapping because of periodic condition
+
+      IF( nn_verbose_level > 0 ) THEN
+         ! store the number of icebergs on this processor at start
+         ibergs_start = icb_utl_count()
+      ENDIF
+
+      ibergs_to_send_e   = 0
+      ibergs_to_send_w   = 0
+      ibergs_to_send_n   = 0
+      ibergs_to_send_s   = 0
+      ibergs_rcvd_from_e = 0
+      ibergs_rcvd_from_w = 0
+      ibergs_rcvd_from_n = 0
+      ibergs_rcvd_from_s = 0
+
+      IF( ASSOCIATED(first_berg) ) THEN      ! Find number of bergs that headed east/west
+         this => first_berg
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            iine = INT( pt%xi + 0.5 )
+            IF( ipe_E >= 0 .AND. iine > mig(nicbei) ) THEN
+               tmpberg => this
+               this => this%next
+               ibergs_to_send_e = ibergs_to_send_e + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to east'
+                  CALL flush( numicb )
+               ENDIF
+               ! deal with periodic case
+               tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp
+               ! now pack it into buffer and delete from list
+               CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e)
+               CALL icb_utl_delete(first_berg, tmpberg)
+            ELSE IF( ipe_W >= 0 .AND. iine < mig(nicbdi) ) THEN
+               tmpberg => this
+               this => this%next
+               ibergs_to_send_w = ibergs_to_send_w + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to west'
+                  CALL flush( numicb )
+               ENDIF
+               ! deal with periodic case
+               tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp )
+               ! now pack it into buffer and delete from list
+               CALL icb_pack_into_buffer( tmpberg, obuffer_w, ibergs_to_send_w)
+               CALL icb_utl_delete(first_berg, tmpberg)
+            ELSE
+               this => this%next
+            ENDIF
+         END DO
+      ENDIF
+      IF( nn_verbose_level >= 3) THEN
+         WRITE(numicb,*) 'bergstep ',nktberg,' send ew: ', ibergs_to_send_e, ibergs_to_send_w
+         CALL flush(numicb)
+      ENDIF
+
+      ! send bergs east and receive bergs from west (ie ones that were sent east) and vice versa
+
+      ! pattern here is copied from lib_mpp code
+
+      SELECT CASE ( nbondi )
+      CASE( -1 )
+         zwebergs(1) = ibergs_to_send_e
+         CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req1)
+         CALL mpprecv( 11, zewbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
+         ibergs_rcvd_from_e = INT( zewbergs(2) )
+      CASE(  0 )
+         zewbergs(1) = ibergs_to_send_w
+         zwebergs(1) = ibergs_to_send_e
+         CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req2)
+         CALL mppsend( 12, zwebergs(1), 1, ipe_E, iml_req3)
+         CALL mpprecv( 11, zewbergs(2), 1 )
+         CALL mpprecv( 12, zwebergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+         IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
+         ibergs_rcvd_from_e = INT( zewbergs(2) )
+         ibergs_rcvd_from_w = INT( zwebergs(2) )
+      CASE(  1 )
+         zewbergs(1) = ibergs_to_send_w
+         CALL mppsend( 11, zewbergs(1), 1, ipe_W, iml_req4)
+         CALL mpprecv( 12, zwebergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
+         ibergs_rcvd_from_w = INT( zwebergs(2) )
+      END SELECT
+      IF( nn_verbose_level >= 3) THEN
+         WRITE(numicb,*) 'bergstep ',nktberg,' recv ew: ', ibergs_rcvd_from_w, ibergs_rcvd_from_e
+         CALL flush(numicb)
+      ENDIF
+
+      SELECT CASE ( nbondi )
+      CASE( -1 )
+         IF( ibergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, ibergs_to_send_e*jp_buffer_width, ipe_E, iml_req1 )
+         IF( ibergs_rcvd_from_e > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_e, ibergs_rcvd_from_e)
+            CALL mpprecv( 13, ibuffer_e%data, ibergs_rcvd_from_e*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_e
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_e, i)
+         ENDDO
+      CASE(  0 )
+         IF( ibergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, ibergs_to_send_w*jp_buffer_width, ipe_W, iml_req2 )
+         IF( ibergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, ibergs_to_send_e*jp_buffer_width, ipe_E, iml_req3 )
+         IF( ibergs_rcvd_from_e > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_e, ibergs_rcvd_from_e)
+            CALL mpprecv( 13, ibuffer_e%data, ibergs_rcvd_from_e*jp_buffer_width )
+         ENDIF
+         IF( ibergs_rcvd_from_w > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_w, ibergs_rcvd_from_w)
+            CALL mpprecv( 14, ibuffer_w%data, ibergs_rcvd_from_w*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+         IF( ibergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_e
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_e, i)
+         END DO
+         DO i = 1, ibergs_rcvd_from_w
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_w, i)
+         ENDDO
+      CASE(  1 )
+         IF( ibergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, ibergs_to_send_w*jp_buffer_width, ipe_W, iml_req4 )
+         IF( ibergs_rcvd_from_w > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_w, ibergs_rcvd_from_w)
+            CALL mpprecv( 14, ibuffer_w%data, ibergs_rcvd_from_w*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_w
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_w, i)
+         END DO
+      END SELECT
+
+      ! Find number of bergs that headed north/south
+      ! (note: this block should technically go ahead of the E/W recv block above
+      !  to handle arbitrary orientation of PEs. But for simplicity, it is
+      !  here to accomodate diagonal transfer of bergs between PEs -AJA)
+
+      IF( ASSOCIATED(first_berg) ) THEN
+         this => first_berg
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            ijne = INT( pt%yj + 0.5 )
+            IF( ipe_N >= 0 .AND. ijne .GT. mjg(nicbej) ) THEN
+               tmpberg => this
+               this => this%next
+               ibergs_to_send_n = ibergs_to_send_n + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to north'
+                  CALL flush( numicb )
+               ENDIF
+               CALL icb_pack_into_buffer( tmpberg, obuffer_n, ibergs_to_send_n)
+               CALL icb_utl_delete(first_berg, tmpberg)
+            ELSE IF( ipe_S >= 0 .AND. ijne .LT. mjg(nicbdj) ) THEN
+               tmpberg => this
+               this => this%next
+               ibergs_to_send_s = ibergs_to_send_s + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for transfer to south'
+                  CALL flush( numicb )
+               ENDIF
+               CALL icb_pack_into_buffer( tmpberg, obuffer_s, ibergs_to_send_s)
+               CALL icb_utl_delete(first_berg, tmpberg)
+            ELSE
+               this => this%next
+            ENDIF
+         END DO
+      ENDIF
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',nktberg,' send ns: ', ibergs_to_send_n, ibergs_to_send_s
+         call flush(numicb)
+      endif
+
+      ! send bergs north
+      ! and receive bergs from south (ie ones sent north)
+
+      SELECT CASE ( nbondj )
+      CASE( -1 )
+         zsnbergs(1) = ibergs_to_send_n
+         CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req1)
+         CALL mpprecv( 15, znsbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
+         ibergs_rcvd_from_n = INT( znsbergs(2) )
+      CASE(  0 )
+         znsbergs(1) = ibergs_to_send_s
+         zsnbergs(1) = ibergs_to_send_n
+         CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req2)
+         CALL mppsend( 16, zsnbergs(1), 1, ipe_N, iml_req3)
+         CALL mpprecv( 15, znsbergs(2), 1 )
+         CALL mpprecv( 16, zsnbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+         IF( l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
+         ibergs_rcvd_from_n = INT( znsbergs(2) )
+         ibergs_rcvd_from_s = INT( zsnbergs(2) )
+      CASE(  1 )
+         znsbergs(1) = ibergs_to_send_s
+         CALL mppsend( 15, znsbergs(1), 1, ipe_S, iml_req4)
+         CALL mpprecv( 16, zsnbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
+         ibergs_rcvd_from_s = INT( zsnbergs(2) )
+      END SELECT
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',nktberg,' recv ns: ', ibergs_rcvd_from_s, ibergs_rcvd_from_n
+         call flush(numicb)
+      endif
+
+      SELECT CASE ( nbondj )
+      CASE( -1 )
+         IF( ibergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, ibergs_to_send_n*jp_buffer_width, ipe_N, iml_req1 )
+         IF( ibergs_rcvd_from_n > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_n, ibergs_rcvd_from_n)
+            CALL mpprecv( 17, ibuffer_n%data, ibergs_rcvd_from_n*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_n
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_n, i)
+         END DO
+      CASE(  0 )
+         IF( ibergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, ibergs_to_send_s*jp_buffer_width, ipe_S, iml_req2 )
+         IF( ibergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, ibergs_to_send_n*jp_buffer_width, ipe_N, iml_req3 )
+         IF( ibergs_rcvd_from_n > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_n, ibergs_rcvd_from_n)
+            CALL mpprecv( 17, ibuffer_n%data, ibergs_rcvd_from_n*jp_buffer_width )
+         ENDIF
+         IF( ibergs_rcvd_from_s > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_s, ibergs_rcvd_from_s)
+            CALL mpprecv( 18, ibuffer_s%data, ibergs_rcvd_from_s*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+         IF( ibergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( iml_req3, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_n
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_n, i)
+         END DO
+         DO i = 1, ibergs_rcvd_from_s
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_s, i)
+         ENDDO
+      CASE(  1 )
+         IF( ibergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, ibergs_to_send_s*jp_buffer_width, ipe_S, iml_req4 )
+         IF( ibergs_rcvd_from_s > 0 ) THEN
+            CALL icb_increase_ibuffer(ibuffer_s, ibergs_rcvd_from_s)
+            CALL mpprecv( 18, ibuffer_s%data, ibergs_rcvd_from_s*jp_buffer_width )
+         ENDIF
+         IF( ibergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( iml_req4, iml_stat, iml_err )
+         DO i = 1, ibergs_rcvd_from_s
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'
+               CALL flush( numicb )
+            ENDIF
+            CALL icb_unpack_from_buffer(first_berg, ibuffer_s, i)
+         END DO
+      END SELECT
+
+      IF( nn_verbose_level > 0 ) THEN
+         ! compare the number of icebergs on this processor from the start to the end
+         ibergs_end = icb_utl_count()
+         i = ( ibergs_rcvd_from_n + ibergs_rcvd_from_s + ibergs_rcvd_from_e + ibergs_rcvd_from_w ) - &
+             ( ibergs_to_send_n + ibergs_to_send_s + ibergs_to_send_e + ibergs_to_send_w )
+         IF( ibergs_end-(ibergs_start+i) .NE. 0 ) THEN
+            WRITE( numicb,*   ) 'send_bergs_to_other_pes: net change in number of icebergs'
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_end=', &
+                                ibergs_end,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_start=', &
+                                ibergs_start,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: delta=', &
+                                i,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: error=', &
+                                ibergs_end-(ibergs_start+i),' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_n=', &
+                                ibergs_to_send_n,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_s=', &
+                                ibergs_to_send_s,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_e=', &
+                                ibergs_to_send_e,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_to_send_w=', &
+                                ibergs_to_send_w,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_n=', &
+                                ibergs_rcvd_from_n,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_s=', &
+                                ibergs_rcvd_from_s,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_e=', &
+                                ibergs_rcvd_from_e,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: ibergs_rcvd_from_w=', &
+                                ibergs_rcvd_from_w,' on PE',narea
+  1000      FORMAT(a,i5,a,i4)
+            CALL ctl_stop('send_bergs_to_other_pes: lost or gained an iceberg or two')
+         ENDIF
+      ENDIF
+
+      ! deal with north fold if we necessary when there is more than one top row processor
+      ! note that for jpni=1 north fold has been dealt with above in call to icb_lbc
+      IF( npolj /= 0 .AND. jpni > 1 ) CALL icb_lbc_mpp_nfld( )
+
+      IF( nn_verbose_level > 0 ) THEN
+         i = 0
+         this => first_berg
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            iine = INT( pt%xi + 0.5 )
+            ijne = INT( pt%yj + 0.5 )
+            IF( iine .LT. mig(nicbdi) .OR. &
+                iine .GT. mig(nicbei) .OR. &
+                ijne .LT. mjg(nicbdj) .OR. &
+                ijne .GT. mjg(nicbej)) THEN
+               i = i + 1
+               WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne
+               WRITE(numicb,*) '                   ', nimpp, njmpp
+               WRITE(numicb,*) '                   ', nicbdi, nicbei, nicbdj, nicbej
+               CALL flush( numicb )
+            ENDIF
+            this => this%next
+         ENDDO ! WHILE
+         CALL mpp_sum(i)
+         IF( i .GT. 0 ) THEN
+            WRITE( numicb,'(a,i4)') 'send_bergs_to_other_pes: # of bergs outside computational domain = ',i
+            CALL ctl_stop('send_bergs_to_other_pes:  there are bergs still in halos!')
+         ENDIF ! root_pe
+      ENDIF ! debug
+      !
+      CALL mppsync()
+      !
+   END SUBROUTINE icb_lbc_mpp
+
+
+   SUBROUTINE icb_lbc_mpp_nfld()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_lbc_mpp_nfld  ***
+      !!
+      !! ** Purpose :   north fold treatment in multi processor exchange
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+      TYPE(iceberg)         , POINTER     :: tmpberg, this
+      TYPE(point)           , POINTER     :: pt
+      INTEGER                             :: ibergs_to_send
+      INTEGER                             :: ibergs_to_rcv
+      INTEGER                             :: iiglo, ijglo, jk, jn
+      INTEGER                             :: ifldproc, iproc, ipts
+      INTEGER                             :: iine, ijne
+      REAL(wp), DIMENSION(2)              :: zsbergs, znbergs
+      INTEGER                             :: iml_req1, iml_req2, iml_err
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) :: iml_stat
+
+      ! set up indices of neighbouring processors
+
+      ! nicbfldproc is a list of unique processor numbers that this processor
+      ! exchanges with (including itself), so we loop over this array; since
+      ! its of fixed size, the first -1 marks end of list of processors
+      !
+      DO jn = 1, jpni
+         IF( nicbfldproc(jn) /= -1 ) THEN
+            ifldproc = nicbfldproc(jn)
+            ibergs_to_send = 0
+   
+            ! Find number of bergs that need to be exchanged
+            ! Pick out exchanges with processor ifldproc
+            ! if ifldproc is this processor then don't send
+            !
+            IF( ASSOCIATED(first_berg) ) THEN
+               this => first_berg
+               DO WHILE (ASSOCIATED(this))
+                  pt => this%current_point
+                  iine = INT( pt%xi + 0.5 )
+                  ijne = INT( pt%yj + 0.5 )
+                  ipts  = nicbfldpts (mi1(iine))
+                  iproc = nicbflddest(mi1(iine))
+                  IF( ijne .GT. mjg(nicbej) ) THEN
+                     IF( iproc == ifldproc ) THEN
+                        !
+                        ! moving across the cut line means both position and
+                        ! velocity must change
+                        ijglo = INT( ipts/nicbpack )
+                        iiglo = ipts - nicbpack*ijglo
+                        pt%xi = iiglo - ( pt%xi - REAL(iine,wp) )
+                        pt%yj = ijglo - ( pt%yj - REAL(ijne,wp) )
+                        pt%uvel = -1._wp * pt%uvel
+                        pt%vvel = -1._wp * pt%vvel
+                        !
+                        ! now remove berg from list and pack it into a buffer
+                        IF( iproc /= narea ) THEN
+                           tmpberg => this
+                           ibergs_to_send = ibergs_to_send + 1
+                           IF( nn_verbose_level >= 4 ) THEN
+                              WRITE(numicb,*) 'bergstep ',nktberg,' packing berg ',tmpberg%number(:),' for north fold'
+                              CALL flush( numicb )
+                           ENDIF
+                           CALL icb_pack_into_buffer( tmpberg, obuffer_f, ibergs_to_send)
+                           CALL icb_utl_delete(first_berg, tmpberg)
+                        ENDIF
+                        !
+                     ENDIF
+                  ENDIF
+                  this => this%next
+               END DO
+            ENDIF
+            if( nn_verbose_level >= 3) then
+               write(numicb,*) 'bergstep ',nktberg,' send nfld: ', ibergs_to_send
+               call flush(numicb)
+            endif
+            !
+            ! if we're in this processor, then we've done everything we need to
+            ! so go on to next element of loop
+            IF( ifldproc == narea ) CYCLE
+   
+            zsbergs(1) = ibergs_to_send
+            CALL mppsend( 21, zsbergs(1), 1, ifldproc-1, iml_req1)
+            CALL mpprecv( 21, znbergs(2), 1 )
+            IF( l_isend ) CALL mpi_wait( iml_req1, iml_stat, iml_err )
+            ibergs_to_rcv = INT( znbergs(2) )
+   
+            ! send bergs
+   
+            IF( ibergs_to_send > 0 )  &
+                CALL mppsend( 12, obuffer_f%data, ibergs_to_send*jp_buffer_width, ifldproc-1, iml_req2 )
+            IF( ibergs_to_rcv  > 0 ) THEN
+               CALL icb_increase_ibuffer(ibuffer_f, ibergs_to_rcv)
+               CALL mpprecv( 12, ibuffer_f%data, ibergs_to_rcv*jp_buffer_width )
+            ENDIF
+            IF( ibergs_to_send > 0 .AND. l_isend ) CALL mpi_wait( iml_req2, iml_stat, iml_err )
+            DO jk = 1, ibergs_to_rcv
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',nktberg,' unpacking berg ',INT(ibuffer_f%data(16,jk)),' from north fold'
+                  CALL flush( numicb )
+               ENDIF
+               CALL icb_unpack_from_buffer(first_berg, ibuffer_f, jk )
+            END DO
+            !
+         ENDIF
+         !
+      END DO
+      !
+   END SUBROUTINE icb_lbc_mpp_nfld
+
+
+   SUBROUTINE icb_pack_into_buffer( berg, pbuff, kb )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER :: berg
+      TYPE(buffer) , POINTER :: pbuff
+      INTEGER               , INTENT(in) :: kb
+      ! 
+      INTEGER ::   k   ! local integer
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ASSOCIATED(pbuff) ) CALL icb_increase_buffer( pbuff, jp_delta_buf )
+      IF( kb .GT. pbuff%size ) CALL icb_increase_buffer( pbuff, jp_delta_buf )
+
+      !! pack points into buffer
+
+      pbuff%data( 1,kb) = berg%current_point%lon
+      pbuff%data( 2,kb) = berg%current_point%lat
+      pbuff%data( 3,kb) = berg%current_point%uvel
+      pbuff%data( 4,kb) = berg%current_point%vvel
+      pbuff%data( 5,kb) = berg%current_point%xi
+      pbuff%data( 6,kb) = berg%current_point%yj
+      pbuff%data( 7,kb) = float(berg%current_point%year)
+      pbuff%data( 8,kb) = berg%current_point%day
+      pbuff%data( 9,kb) = berg%current_point%mass
+      pbuff%data(10,kb) = berg%current_point%thickness
+      pbuff%data(11,kb) = berg%current_point%width
+      pbuff%data(12,kb) = berg%current_point%length
+      pbuff%data(13,kb) = berg%current_point%mass_of_bits
+      pbuff%data(14,kb) = berg%current_point%heat_density
+
+      pbuff%data(15,kb) = berg%mass_scaling
+      DO k=1,nkounts
+         pbuff%data(15+k,kb) = REAL( berg%number(k), wp )
+      END DO
+      !
+   END SUBROUTINE icb_pack_into_buffer
+
+
+   SUBROUTINE icb_unpack_from_buffer(first, pbuff, kb)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      TYPE(iceberg),             POINTER :: first
+      TYPE(buffer) ,             POINTER :: pbuff
+      INTEGER      , INTENT(in)          :: kb
+      ! 
+      TYPE(iceberg)                      :: currentberg
+      TYPE(point)                        :: pt
+      INTEGER                            :: ik
+      !!----------------------------------------------------------------------
+      !
+      pt%lon            =      pbuff%data( 1,kb)
+      pt%lat            =      pbuff%data( 2,kb)
+      pt%uvel           =      pbuff%data( 3,kb)
+      pt%vvel           =      pbuff%data( 4,kb)
+      pt%xi             =      pbuff%data( 5,kb)
+      pt%yj             =      pbuff%data( 6,kb)
+      pt%year           = INT( pbuff%data( 7,kb) )
+      pt%day            =      pbuff%data( 8,kb)
+      pt%mass           =      pbuff%data( 9,kb)
+      pt%thickness      =      pbuff%data(10,kb)
+      pt%width          =      pbuff%data(11,kb)
+      pt%length         =      pbuff%data(12,kb)
+      pt%mass_of_bits   =      pbuff%data(13,kb)
+      pt%heat_density   =      pbuff%data(14,kb)
+
+      currentberg%mass_scaling =      pbuff%data(15,kb)
+      DO ik = 1, nkounts
+         currentberg%number(ik) = INT( pbuff%data(15+ik,kb) )
+      END DO
+      !
+      CALL icb_utl_add(currentberg, pt )
+      !
+   END SUBROUTINE icb_unpack_from_buffer
+
+
+   SUBROUTINE icb_increase_buffer(old,kdelta)
+      !!----------------------------------------------------------------------
+      TYPE(buffer), POINTER    :: old
+      INTEGER     , INTENT(in) :: kdelta
+      ! 
+      TYPE(buffer), POINTER ::   new
+      INTEGER ::   inew_size
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ASSOCIATED(old) ) THEN   ;   inew_size = kdelta
+      ELSE                               ;   inew_size = old%size + kdelta
+      ENDIF
+      ALLOCATE( new )
+      ALLOCATE( new%data( jp_buffer_width, inew_size) )
+      new%size = inew_size
+      IF( ASSOCIATED(old) ) THEN
+         new%data(:,1:old%size) = old%data(:,1:old%size)
+         DEALLOCATE(old%data)
+         DEALLOCATE(old)
+      ENDIF
+      old => new
+      !
+   END SUBROUTINE icb_increase_buffer
+
+
+   SUBROUTINE icb_increase_ibuffer(old,kdelta)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      TYPE(buffer),            POINTER :: old
+      INTEGER     , INTENT(in)         :: kdelta
+      !
+      TYPE(buffer),            POINTER :: new
+      INTEGER                          :: inew_size, iold_size
+      !!----------------------------------------------------------------------
+
+      IF( .NOT. ASSOCIATED(old) ) THEN
+         inew_size = kdelta + jp_delta_buf
+         iold_size = 0
+      ELSE
+         iold_size = old%size
+         IF( kdelta .LT. old%size ) THEN
+            inew_size = old%size + kdelta
+         ELSE
+            inew_size = kdelta + jp_delta_buf
+         ENDIF
+      ENDIF
+
+      IF( iold_size .NE. inew_size ) THEN
+         ALLOCATE( new )
+         ALLOCATE( new%data( jp_buffer_width, inew_size) )
+         new%size = inew_size
+         IF( ASSOCIATED(old) ) THEN
+            new%data(:,1:old%size) = old%data(:,1:old%size)
+            DEALLOCATE(old%data)
+            DEALLOCATE(old)
+         ENDIF
+         old => new
+        !WRITE( numicb,*) 'icb_increase_ibuffer',narea,' increased to',inew_size
+      ENDIF
+      !
+   END SUBROUTINE icb_increase_ibuffer
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default case:            Dummy module        share memory computing
+   !!----------------------------------------------------------------------
+   SUBROUTINE icb_lbc_mpp()
+      WRITE(numout,*) 'icb_lbc_mpp: You should not have seen this message!!'
+   END SUBROUTINE icb_lbc_mpp
+
+#endif
+
+   !!======================================================================
+END MODULE icblbc
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 3609)
@@ -0,0 +1,414 @@
+MODULE icbrst
+
+   !!======================================================================
+   !!                       ***  MODULE  icbrst  ***
+   !! Ocean physics:  read and write iceberg restart files
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !  2011-04  (Alderson)       Restore restart routine
+   !!            -    !                            Currently needs a fixed processor
+   !!            -    !                            layout between restarts
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_rst_read    : read restart file
+   !!   icb_rst_write   : write restart file
+   !!----------------------------------------------------------------------
+   USE par_oce        ! NEMO parameters
+   USE dom_oce        ! NEMO domain
+   USE in_out_manager ! NEMO IO routines
+   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
+   USE netcdf         ! netcdf routines for IO
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_rst_read    ! routine called in icbini.F90 module
+   PUBLIC   icb_rst_write   ! routine called in icbstp.F90 module
+   
+   INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid
+   INTEGER ::   nmassid, nthicknessid, nwidthid, nlengthid
+   INTEGER ::   nyearid, ndayid
+   INTEGER ::   nscaling_id, nmass_of_bits_id, nheat_density_id, numberid
+   INTEGER ::   nsiceid, nsheatid, ncalvid, ncalvhid, nkountid
+   INTEGER ::   nret, ncid, nc_dim
+   
+   INTEGER,  DIMENSION(3)                  :: nstrt3, nlngth3
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_rst_read()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_rst_read  ***
+      !!
+      !! ** Purpose :   read a iceberg restart file
+      !!      NB: for this version, we just read back in the restart for this processor
+      !!      so we cannot change the processor layout currently with iceberg code
+      !!----------------------------------------------------------------------
+      INTEGER                      ::   idim, ivar, iatt
+      INTEGER                      ::   jn, iunlim_dim, ibergs_in_file
+      INTEGER                      ::   iclass
+      INTEGER, DIMENSION(1)        ::   istrt, ilngth, idata
+      INTEGER, DIMENSION(2)        ::   istrt2, ilngth2
+      INTEGER, DIMENSION(nkounts)  ::   idata2
+      REAL(wp), DIMENSION(1)       ::   zdata                                         ! need 1d array to read in with
+                                                                                            ! start and count arrays
+      LOGICAL                      ::   ll_found_restart
+      CHARACTER(len=80)            ::   cl_filename
+      CHARACTER(len=NF90_MAX_NAME) ::   cl_dname
+      TYPE(iceberg)                ::   localberg ! NOT a pointer but an actual local variable
+      TYPE(point)                  ::   localpt   ! NOT a pointer but an actual local variable
+      !!----------------------------------------------------------------------
+
+      ! Find a restart file
+      cl_filename = ' '
+      IF ( lk_mpp ) THEN
+         cl_filename = ' '
+         WRITE( cl_filename, '("restart_icebergs_",I4.4,".nc")' ) narea-1
+         INQUIRE( file=TRIM(cl_filename), exist=ll_found_restart )
+      ELSE
+         cl_filename = 'restart_icebergs.nc'
+         INQUIRE( file=TRIM(cl_filename), exist=ll_found_restart )
+      ENDIF
+
+      IF ( .NOT. ll_found_restart) THEN                     ! only do the following if a file was found
+         CALL ctl_stop('icebergs: no restart file found')
+      ENDIF
+
+      IF (nn_verbose_level >= 0 .AND. lwp)  &
+         WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM(cl_filename)
+
+      nret = NF90_OPEN(TRIM(cl_filename), NF90_NOWRITE, ncid)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed')
+
+      nret = nf90_inquire(ncid, idim, ivar, iatt, iunlim_dim)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed')
+
+      IF( iunlim_dim .NE. -1) THEN
+
+         nret = nf90_inquire_dimension(ncid, iunlim_dim, cl_dname, ibergs_in_file)
+         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed')
+
+         nret = NF90_INQ_VARID(ncid, 'number', numberid)
+         nret = NF90_INQ_VARID(ncid, 'mass_scaling', nscaling_id)
+         nret = NF90_INQ_VARID(ncid, 'xi', nxid)
+         nret = NF90_INQ_VARID(ncid, 'yj', nyid)
+         nret = NF90_INQ_VARID(ncid, 'lon', nlonid)
+         nret = NF90_INQ_VARID(ncid, 'lat', nlatid)
+         nret = NF90_INQ_VARID(ncid, 'uvel', nuvelid)
+         nret = NF90_INQ_VARID(ncid, 'vvel', nvvelid)
+         nret = NF90_INQ_VARID(ncid, 'mass', nmassid)
+         nret = NF90_INQ_VARID(ncid, 'thickness', nthicknessid)
+         nret = NF90_INQ_VARID(ncid, 'width', nwidthid)
+         nret = NF90_INQ_VARID(ncid, 'length', nlengthid)
+         nret = NF90_INQ_VARID(ncid, 'year', nyearid)
+         nret = NF90_INQ_VARID(ncid, 'day', ndayid)
+         nret = NF90_INQ_VARID(ncid, 'mass_of_bits', nmass_of_bits_id)
+         nret = NF90_INQ_VARID(ncid, 'heat_density', nheat_density_id)
+
+         ilngth(1) = 1
+         istrt2(1) = 1
+         ilngth2(1) = nkounts
+         ilngth2(2) = 1
+         DO jn=1, ibergs_in_file
+
+            istrt(1) = jn
+            istrt2(2) = jn
+
+            nret = NF90_GET_VAR(ncid, numberid, idata2, istrt2, ilngth2 )
+            localberg%number(:) = idata2(:)
+
+            nret = NF90_GET_VAR(ncid, nscaling_id, zdata, istrt, ilngth )
+            localberg%mass_scaling = zdata(1)
+
+            nret = NF90_GET_VAR(ncid, nlonid, zdata, istrt, ilngth)
+            localpt%lon = zdata(1)
+            nret = NF90_GET_VAR(ncid, nlatid, zdata, istrt, ilngth)
+            localpt%lat = zdata(1)
+            IF (nn_verbose_level >= 2 .AND. lwp) THEN
+               WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ',jn,' is at ', &
+                                              localpt%lon,localpt%lat,' on PE ',narea-1
+            ENDIF
+            nret = NF90_GET_VAR(ncid, nxid, zdata, istrt, ilngth)
+            localpt%xi = zdata(1)
+            nret = NF90_GET_VAR(ncid, nyid, zdata, istrt, ilngth)
+            localpt%yj = zdata(1)
+            nret = NF90_GET_VAR(ncid, nuvelid, zdata, istrt, ilngth )
+            localpt%uvel = zdata(1)
+            nret = NF90_GET_VAR(ncid, nvvelid, zdata, istrt, ilngth )
+            localpt%vvel = zdata(1)
+            nret = NF90_GET_VAR(ncid, nmassid, zdata, istrt, ilngth )
+            localpt%mass = zdata(1)
+            nret = NF90_GET_VAR(ncid, nthicknessid, zdata, istrt, ilngth )
+            localpt%thickness = zdata(1)
+            nret = NF90_GET_VAR(ncid, nwidthid, zdata, istrt, ilngth )
+            localpt%width = zdata(1)
+            nret = NF90_GET_VAR(ncid, nlengthid, zdata, istrt, ilngth )
+            localpt%length = zdata(1)
+            nret = NF90_GET_VAR(ncid, nyearid, idata, istrt, ilngth )
+            localpt%year = idata(1)
+            nret = NF90_GET_VAR(ncid, ndayid, zdata, istrt, ilngth )
+            localpt%day = zdata(1)
+            nret = NF90_GET_VAR(ncid, nmass_of_bits_id, zdata, istrt, ilngth )
+            localpt%mass_of_bits = zdata(1)
+            nret = NF90_GET_VAR(ncid, nheat_density_id, zdata, istrt, ilngth )
+            localpt%heat_density = zdata(1)
+            !
+            CALL icb_utl_add( localberg, localpt )
+         END DO
+         !
+      ENDIF
+
+      nret = NF90_INQ_DIMID( ncid, 'c', nc_dim )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed')
+
+      nret = NF90_INQUIRE_DIMENSION( ncid, nc_dim, cl_dname, iclass )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed')
+
+      nret = NF90_INQ_VARID(ncid, 'kount'       , nkountid)
+      nret = NF90_INQ_VARID(ncid, 'calving'     , ncalvid)
+      nret = NF90_INQ_VARID(ncid, 'calving_hflx', ncalvhid)
+      nret = NF90_INQ_VARID(ncid, 'stored_ice'  , nsiceid)
+      nret = NF90_INQ_VARID(ncid, 'stored_heat' , nsheatid)
+
+      nstrt3(1) = 1
+      nstrt3(2) = 1
+      nlngth3(1) = jpi
+      nlngth3(2) = jpj
+      nlngth3(3) = 1
+
+      DO jn = 1, iclass
+         nstrt3(3) = jn
+         nret      = NF90_GET_VAR( ncid, nsiceid , griddata, nstrt3, nlngth3 )
+         berg_grid%stored_ice(:,:,jn) = griddata(:,:,1)
+      END DO
+
+      nret = NF90_GET_VAR( ncid, ncalvid , src_calving          (:,:) )
+      nret = NF90_GET_VAR( ncid, ncalvhid, src_calving_hflx     (:,:) )
+      nret = NF90_GET_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) )
+      nret = NF90_GET_VAR( ncid, nkountid, idata2(:) )
+      num_bergs(:) = idata2(:)
+
+      ! Finish up
+      nret = NF90_CLOSE(ncid)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed')
+
+      ! Sanity check
+      jn = icb_utl_count()
+      IF (nn_verbose_level >= 0)   &
+         WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',jn,' on PE',narea-1
+      IF( lk_mpp ) THEN
+         CALL mpp_sum(ibergs_in_file)
+         CALL mpp_sum(jn)
+      ENDIF
+      IF(lwp)   WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, read_restart_bergs: there were',ibergs_in_file,   &
+         &                                    ' bergs in the restart file and', jn,' bergs have been read'
+      !
+      IF( lwp .and. nn_verbose_level >= 0)  WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'
+      !
+   END SUBROUTINE icb_rst_read
+
+
+   SUBROUTINE icb_rst_write( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE icb_rst_write  ***
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt
+      !
+      INTEGER ::   jn   ! dummy loop index
+      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
+      CHARACTER(len=80)      :: cl_filename
+      TYPE(iceberg), POINTER :: this
+      TYPE(point)  , POINTER :: pt
+      !!----------------------------------------------------------------------
+
+      IF( lk_mpp ) THEN
+         WRITE(cl_filename,'("icebergs_",I8.8,"_restart_",I4.4,".nc")') kt, narea-1
+      ELSE
+         WRITE(cl_filename,'("icebergs_",I8.8,"_restart.nc")') kt
+      ENDIF
+      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_filename)
+
+      nret = NF90_CREATE(TRIM(cl_filename), NF90_CLOBBER, ncid)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed')
+
+      ! Dimensions
+      nret = NF90_DEF_DIM(ncid, 'x', jpi, ix_dim)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim x failed')
+
+      nret = NF90_DEF_DIM(ncid, 'y', jpj, iy_dim)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim y failed')
+
+      nret = NF90_DEF_DIM(ncid, 'c', nclasses, nc_dim)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim c failed')
+
+      nret = NF90_DEF_DIM(ncid, 'k', nkounts, ik_dim)
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed')
+
+      IF (associated(first_berg)) then
+         nret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED, in_dim)
+         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim n failed')
+      ENDIF
+
+      ! Variables
+      nret = NF90_DEF_VAR(ncid, 'kount'       , NF90_INT   , (/ ik_dim /), nkountid)
+      nret = NF90_DEF_VAR(ncid, 'calving'     , NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvid)
+      nret = NF90_DEF_VAR(ncid, 'calving_hflx', NF90_DOUBLE, (/ ix_dim, iy_dim /), ncalvhid)
+      nret = NF90_DEF_VAR(ncid, 'stored_ice'  , NF90_DOUBLE, (/ ix_dim, iy_dim, nc_dim /), nsiceid)
+      nret = NF90_DEF_VAR(ncid, 'stored_heat' , NF90_DOUBLE, (/ ix_dim, iy_dim /), nsheatid)
+
+      ! Attributes
+      nret = NF90_PUT_ATT(ncid, ncalvid , 'long_name', 'iceberg calving')
+      nret = NF90_PUT_ATT(ncid, ncalvid , 'units', 'some')
+      nret = NF90_PUT_ATT(ncid, ncalvhid, 'long_name', 'heat flux associated with iceberg calving')
+      nret = NF90_PUT_ATT(ncid, ncalvhid, 'units', 'some')
+      nret = NF90_PUT_ATT(ncid, nsiceid , 'long_name', 'stored ice used to calve icebergs')
+      nret = NF90_PUT_ATT(ncid, nsiceid , 'units', 'kg/s')
+      nret = NF90_PUT_ATT(ncid, nsheatid, 'long_name', 'heat in stored ice used to calve icebergs')
+      nret = NF90_PUT_ATT(ncid, nsheatid, 'units', 'J/kg/s')
+
+      IF ( ASSOCIATED(first_berg) ) THEN
+
+         ! Only add berg variables for this PE if we have anything to say
+
+         ! Variables
+         nret = NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, in_dim, nlonid)
+         nret = NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, in_dim, nlatid)
+         nret = NF90_DEF_VAR(ncid, 'xi', NF90_DOUBLE, in_dim, nxid)
+         nret = NF90_DEF_VAR(ncid, 'yj', NF90_DOUBLE, in_dim, nyid)
+         nret = NF90_DEF_VAR(ncid, 'uvel', NF90_DOUBLE, in_dim, nuvelid)
+         nret = NF90_DEF_VAR(ncid, 'vvel', NF90_DOUBLE, in_dim, nvvelid)
+         nret = NF90_DEF_VAR(ncid, 'mass', NF90_DOUBLE, in_dim, nmassid)
+         nret = NF90_DEF_VAR(ncid, 'thickness', NF90_DOUBLE, in_dim, nthicknessid)
+         nret = NF90_DEF_VAR(ncid, 'width', NF90_DOUBLE, in_dim, nwidthid)
+         nret = NF90_DEF_VAR(ncid, 'length', NF90_DOUBLE, in_dim, nlengthid)
+         nret = NF90_DEF_VAR(ncid, 'number', NF90_INT, (/ik_dim,in_dim/), numberid)
+         nret = NF90_DEF_VAR(ncid, 'year', NF90_INT, in_dim, nyearid)
+         nret = NF90_DEF_VAR(ncid, 'day', NF90_DOUBLE, in_dim, ndayid)
+         nret = NF90_DEF_VAR(ncid, 'mass_scaling', NF90_DOUBLE, in_dim, nscaling_id)
+         nret = NF90_DEF_VAR(ncid, 'mass_of_bits', NF90_DOUBLE, in_dim, nmass_of_bits_id)
+         nret = NF90_DEF_VAR(ncid, 'heat_density', NF90_DOUBLE, in_dim, nheat_density_id)
+
+         ! Attributes
+         nret = NF90_PUT_ATT(ncid, nlonid, 'long_name', 'longitude')
+         nret = NF90_PUT_ATT(ncid, nlonid, 'units', 'degrees_E')
+         nret = NF90_PUT_ATT(ncid, nlatid, 'long_name', 'latitude')
+         nret = NF90_PUT_ATT(ncid, nlatid, 'units', 'degrees_N')
+         nret = NF90_PUT_ATT(ncid, nxid, 'long_name', 'x grid box position')
+         nret = NF90_PUT_ATT(ncid, nxid, 'units', 'fractional')
+         nret = NF90_PUT_ATT(ncid, nyid, 'long_name', 'y grid box position')
+         nret = NF90_PUT_ATT(ncid, nyid, 'units', 'fractional')
+         nret = NF90_PUT_ATT(ncid, nuvelid, 'long_name', 'zonal velocity')
+         nret = NF90_PUT_ATT(ncid, nuvelid, 'units', 'm/s')
+         nret = NF90_PUT_ATT(ncid, nvvelid, 'long_name', 'meridional velocity')
+         nret = NF90_PUT_ATT(ncid, nvvelid, 'units', 'm/s')
+         nret = NF90_PUT_ATT(ncid, nmassid, 'long_name', 'mass')
+         nret = NF90_PUT_ATT(ncid, nmassid, 'units', 'kg')
+         nret = NF90_PUT_ATT(ncid, nthicknessid, 'long_name', 'thickness')
+         nret = NF90_PUT_ATT(ncid, nthicknessid, 'units', 'm')
+         nret = NF90_PUT_ATT(ncid, nwidthid, 'long_name', 'width')
+         nret = NF90_PUT_ATT(ncid, nwidthid, 'units', 'm')
+         nret = NF90_PUT_ATT(ncid, nlengthid, 'long_name', 'length')
+         nret = NF90_PUT_ATT(ncid, nlengthid, 'units', 'm')
+         nret = NF90_PUT_ATT(ncid, numberid, 'long_name', 'iceberg number on this processor')
+         nret = NF90_PUT_ATT(ncid, numberid, 'units', 'count')
+         nret = NF90_PUT_ATT(ncid, nyearid, 'long_name', 'calendar year of calving event')
+         nret = NF90_PUT_ATT(ncid, nyearid, 'units', 'years')
+         nret = NF90_PUT_ATT(ncid, ndayid, 'long_name', 'year day of calving event')
+         nret = NF90_PUT_ATT(ncid, ndayid, 'units', 'days')
+         nret = NF90_PUT_ATT(ncid, nscaling_id, 'long_name', 'scaling factor for mass of calving berg')
+         nret = NF90_PUT_ATT(ncid, nscaling_id, 'units', 'none')
+         nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'long_name', 'mass of bergy bits')
+         nret = NF90_PUT_ATT(ncid, nmass_of_bits_id, 'units', 'kg')
+         nret = NF90_PUT_ATT(ncid, nheat_density_id, 'long_name', 'heat density')
+         nret = NF90_PUT_ATT(ncid, nheat_density_id, 'units', 'J/kg')
+
+      ENDIF ! associated(first_berg)
+
+      ! End define mode
+      nret = NF90_ENDDEF(ncid)
+
+      ! --------------------------------
+      ! now write some data
+
+      nstrt3(1) = 1
+      nstrt3(2) = 1
+      nlngth3(1) = jpi
+      nlngth3(2) = jpj
+      nlngth3(3) = 1
+
+      DO jn=1,nclasses
+         griddata(:,:,1) = berg_grid%stored_ice(:,:,jn)
+         nstrt3(3) = jn
+         nret = NF90_PUT_VAR( ncid, nsiceid, griddata, nstrt3, nlngth3 )
+         IF (nret .ne. NF90_NOERR) THEN
+            IF( lwp ) WRITE(numout,*) TRIM(NF90_STRERROR( nret ))
+            CALL ctl_stop('icebergs, write_restart: nf_put_var stored_ice failed')
+         ENDIF
+      ENDDO
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_ice  written'
+
+      nret = NF90_PUT_VAR( ncid, nkountid, num_bergs(:) )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var kount failed')
+
+      nret = NF90_PUT_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed')
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_heat written'
+
+      nret = NF90_PUT_VAR( ncid, ncalvid , src_calving(:,:) )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving failed')
+      nret = NF90_PUT_VAR( ncid, ncalvhid, src_calving_hflx(:,:) )
+      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed')
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: calving written'
+
+      IF ( ASSOCIATED(first_berg) ) THEN
+
+         ! Write variables
+         ! just write out the current point of the trajectory
+
+         this => first_berg
+         jn = 0
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            jn=jn+1
+
+            nret = NF90_PUT_VAR(ncid, numberid, this%number, (/1,jn/), (/nkounts,1/) )
+            nret = NF90_PUT_VAR(ncid, nscaling_id, this%mass_scaling, (/ jn /) )
+
+            nret = NF90_PUT_VAR(ncid, nlonid, pt%lon, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nlatid, pt%lat, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nxid, pt%xi, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nyid, pt%yj, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nuvelid, pt%uvel, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nvvelid, pt%vvel, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nmassid, pt%mass, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nthicknessid, pt%thickness, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nwidthid, pt%width, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nlengthid, pt%length, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nyearid, pt%year, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, ndayid, pt%day, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nmass_of_bits_id, pt%mass_of_bits, (/ jn /) )
+            nret = NF90_PUT_VAR(ncid, nheat_density_id, pt%heat_density, (/ jn /) )
+
+            this=>this%next
+         END DO
+         !
+      ENDIF ! associated(first_berg)
+
+      ! Finish up
+      nret = NF90_CLOSE(ncid)
+      IF (nret /= NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_close failed')
+      !
+   END SUBROUTINE icb_rst_write
+   !
+END MODULE icbrst
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 3609)
@@ -0,0 +1,175 @@
+MODULE icbstp
+
+   !!======================================================================
+   !!                       ***  MODULE  icbstp  ***
+   !! Icebergs:  initialise variables for iceberg tracking
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !                            Move budgets to icbdia routine
+   !!            -    !  2011-05  (Alderson)       Add call to copy forcing arrays
+   !!            -    !                            into icb copies with haloes
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_stp       : start iceberg tracking
+   !!   icb_end       : end   iceberg tracking
+   !!----------------------------------------------------------------------
+   USE par_oce        ! nemo parameters
+   USE dom_oce        ! ocean domain
+   USE sbc_oce        ! ocean surface forcing
+   USE phycst
+   USE in_out_manager ! nemo IO
+   USE lib_mpp
+   USE iom
+   USE timing         ! timing
+
+   USE icb_oce        ! define iceberg arrays
+   USE icbini         ! iceberg initialisation routines
+   USE icbutl         ! iceberg utility routines
+   USE icbrst         ! iceberg restart routines
+   USE icbdyn         ! iceberg dynamics (ie advection) routines
+   USE icbclv         ! iceberg calving routines
+   USE icbthm         ! iceberg thermodynamics routines
+   USE icblbc         ! iceberg lateral boundary routines (including mpp)
+   USE icbtrj         ! iceberg trajectory I/O routines
+   USE icbdia         ! iceberg budget
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_stp        ! routine called in sbcmod.F90 module
+   PUBLIC   icb_end        ! routine called in nemogcm.F90 module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE icb_stp( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_stp  ***
+      !!
+      !! ** Purpose :   iceberg time stepping.
+      !!
+      !! ** Method  : - top level routine to do things in the correct order
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! time step index
+      !
+      LOGICAL ::   ll_sample_traj, ll_budget, ll_verbose   ! local logical
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('icb_stp')
+
+      !! start of timestep housekeeping
+
+      nktberg = kt
+
+      IF( nn_test_icebergs < 0 ) THEN      ! read calving data
+         !
+         CALL fld_read ( kt, 1, sf_icb )
+         src_calving(:,:)      = sf_icb(1)%fnow(:,:,1)    ! calving in km^3/year (water equivalent)
+         src_calving_hflx(:,:) = 0._wp                    ! NO heat flux for now
+         !
+      ENDIF
+
+      berg_grid%floating_melt(:,:) = 0._wp
+
+      ! anything that needs to be reset to zero each timestep for budgets is dealt with here
+      CALL icb_dia_step()
+
+      ll_verbose = .FALSE.
+      IF( nn_verbose_write > 0 .AND. &
+          MOD(kt-1,nn_verbose_write ) == 0 )   ll_verbose = nn_verbose_level >= 0
+
+      ! write out time
+      IF( ll_verbose ) WRITE(numicb,9100) nktberg, ndastp, nsec_day
+ 9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8)
+
+      ! copy nemo forcing arrays into iceberg versions with extra halo
+      ! only necessary for variables not on T points
+      CALL icb_utl_copy()
+
+      !!----------------------------------------------------------------------
+      !! process icebergs
+
+                                     CALL icb_clv_flx( kt )   ! Accumulate ice from calving
+
+                                     CALL icb_clv()           ! Calve excess stored ice into icebergs
+
+
+!                               !==  For each berg, evolve  ==!
+      !
+      IF( ASSOCIATED(first_berg) )   CALL icb_dyn()           ! ice berg dynamics
+
+      IF( lk_mpp ) THEN          ;   CALL icb_lbc_mpp()       ! Send bergs to other PEs
+      ELSE                       ;   CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+
+      IF( ASSOCIATED(first_berg) )   CALL icb_thm( kt )       ! Ice berg thermodynamics (melting) + rolling
+
+      !!----------------------------------------------------------------------
+      !! end of timestep housekeeping
+
+      ll_sample_traj = .FALSE.
+      IF( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 )   ll_sample_traj = .TRUE.
+      IF( ll_sample_traj .AND.   &
+          ASSOCIATED(first_berg) )   CALL icb_trj_write( kt )  ! For each berg, record trajectory
+
+      ! Gridded diagnostics
+      ! To get these iom_put's and those preceding to actually do something
+      ! use key_iomput in cpp file and create content for XML file
+
+      CALL iom_put( "calving"           , berg_grid%calving      (:,:)   )  ! 'calving mass input'
+      CALL iom_put( "berg_floating_melt", berg_grid%floating_melt(:,:)   )  ! 'Melt rate of icebergs + bits' , 'kg/m2/s'
+      CALL iom_put( "berg_stored_ice"   , berg_grid%stored_ice   (:,:,:) )  ! 'Accumulated ice mass by class', 'kg'
+
+      ! store mean budgets
+      CALL icb_dia_put()
+
+      ! Dump icebergs to screen
+      if ( nn_verbose_level >= 2 )   CALL icb_utl_print( 'icb_stp, status', kt )
+
+      ! Diagnose budgets
+      ll_budget = .FALSE.
+      IF( nn_verbose_write > 0 .AND. MOD(kt-1,nn_verbose_write) == 0 )   ll_budget = ln_bergdia
+      CALL icb_dia( ll_budget )
+
+      IF( MOD(kt,nn_stock) == 0 ) THEN
+         CALL icb_rst_write( kt )
+         IF( nn_sample_rate > 0 )   CALL icb_trj_sync()
+      ENDIF
+
+      IF( nn_timing == 1 ) CALL timing_stop('icb_stp')
+      !
+   END SUBROUTINE icb_stp
+
+
+   SUBROUTINE icb_end( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_end  ***
+      !!
+      !! ** Purpose :   close iceberg files
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )  :: kt
+      !!----------------------------------------------------------------------
+
+      ! only write a restart if not done in icb_stp
+      IF( MOD(kt,nn_stock) .NE. 0 ) CALL icb_rst_write( kt )
+
+      ! finish with trajectories if they were written
+      IF( nn_sample_rate .GT. 0 ) CALL icb_trj_end()
+
+      IF(lwp)   WRITE(numout,'(a,i6)') 'icebergs: icb_end complete', narea
+      CALL flush( numicb )
+      CLOSE( numicb )
+      !
+   END SUBROUTINE icb_end
+
+   !!-------------------------------------------------------------------------
+
+END MODULE icbstp
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 3609)
@@ -0,0 +1,216 @@
+MODULE icbthm
+
+   !!======================================================================
+   !!                       ***  MODULE  icbthm  ***
+   !! Icebergs:  thermodynamics routines for icebergs
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!            -    !  2011-05  (Alderson)       Use tmask instead of tmask_i
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_thm : initialise
+   !!             reference for equations - M = Martin + Adcroft, OM 34, 2010
+   !!----------------------------------------------------------------------
+   USE par_oce        ! NEMO parameters
+   USE dom_oce        ! NEMO domain
+   USE in_out_manager ! NEMO IO routines, numout in particular
+   USE lib_mpp        ! NEMO MPI routines, ctl_stop in particular
+   USE phycst         ! NEMO physical constants
+   USE sbc_oce
+
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+   USE icbdia         ! iceberg budget routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_thm ! routine called in icbstp.F90 module
+
+CONTAINS
+
+   SUBROUTINE icb_thm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_thm  ***
+      !!
+      !! ** Purpose :   compute the iceberg thermodynamics.
+      !!
+      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg
+      !
+      INTEGER  ::   ii, ij
+      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn
+      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv
+      REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat
+      REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
+      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
+      TYPE(iceberg), POINTER ::   this, next
+      TYPE(point)  , POINTER ::   pt
+      !!----------------------------------------------------------------------
+      !
+      z1_rday = 1._wp / rday
+      
+      ! we're either going to ignore berg fresh water melt flux and associated heat
+      ! or we pass it into the ocean, so at this point we set them both to zero,
+      ! accumulate the contributions to them from each iceberg in the while loop following
+      ! and then pass them (or not) to the ocean
+      !
+      berg_grid%floating_melt(:,:) = 0._wp
+      berg_grid%calving_hflx(:,:)  = 0._wp
+    
+      this => first_berg
+      DO WHILE( associated(this) )
+         !
+         pt => this%current_point
+         nknberg = this%number(1)
+         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, &
+         &                    pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, &
+         &                 pt%sst, pt%cn, pt%hi, zff )
+         !
+         zSST = pt%sst
+         zIC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
+         zM   = pt%mass
+         zT   = pt%thickness                               ! total thickness
+       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
+       ! F   = zT - D ! freeboard
+         zW   = pt%width
+         zL   = pt%length
+         zxi  = pt%xi                                      ! position in (i,j) referential
+         zyj  = pt%yj
+         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg
+         ii  = mi1( ii )
+         ij  = INT( zyj + 0.5 )              
+         ij  = mj1( ij )
+         zVol = zT * zW * zL
+         zdt = berg_dt   ;   z1_dt = 1._wp / zdt
+
+         ! Environment
+         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
+         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
+         zSs  = 1.5 * SQRT( zdva ) + 0.1 * zdva                ! Sea state      (eqn M.A9)
+
+         ! Melt rates in m/s (i.e. division by rday)
+         zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2)            , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10)
+         zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2)      , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 )
+         zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+cos(rpi*(zIC**3))) , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
+
+         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
+            zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
+            znVol  = zTn * zW * zL                        ! new volume (m^3)
+            zMnew1 = (znVol/zVol) * zM                    ! new mass (kg)
+            zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
+            !
+            zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
+            zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
+            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
+            zMnew2 = (znVol/zVol) * zM                    ! new mass (kg)
+            zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
+            !
+            zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
+            zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
+            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
+            zMnew  = ( znVol / zVol ) * zM                ! new mass (kg)
+            zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
+            zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
+            !
+         ELSE                                         ! Update dimensions of berg
+            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )         ! (m)
+            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )         ! (m)
+            zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
+            ! Update volume and mass of berg
+            znVol = zTn*zWn*zLn                           ! (m^3)
+            zMnew = (znVol/zVol)*zM                       ! (kg)
+            zdM   = zM - zMnew                           ! (kg)
+            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt         ! approx. mass loss to basal melting (kg)
+            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt         ! approx. mass lost to erosion (kg)
+            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt         ! approx. mass loss to buoyant convection (kg)
+         ENDIF
+
+         IF( rn_bits_erosion_fraction > 0._wp ) THEN      ! Bergy bits
+            !
+            zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
+            zdMbitsE = rn_bits_erosion_fraction * zdMe                        ! change in mass of bits (kg)
+            znMbits  = zMbits + zdMbitsE                                               ! add new bergy bits to mass (kg)
+            zLbits   = MIN( zL, zW, zT, 40._wp )                                        ! assume bergy bits are smallest dimension or 40 meters
+            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                           ! Effective bottom area (assuming T=Lbits)
+            zMbb     = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday    ! Basal turbulent melting (for bits)
+            zMbb     = rn_rho_bergs * zAbits * zMbb                                 ! in kg/s
+            zdMbitsM = MIN( zMbb*zdt , znMbits )                                       ! bergy bits mass lost to melting (kg)
+            znMbits  = znMbits-zdMbitsM                                                ! remove mass lost to bergy bits melt
+            IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
+               zdMbitsM = zdMbitsM + znMbits                                           ! instantly melt all the bergy bits
+               znMbits  = 0._wp
+            ENDIF
+         ELSE                                                     ! No bergy bits
+            zAbits   = 0._wp
+            zdMbitsE = 0._wp
+            zdMbitsM = 0._wp
+            znMbits  = pt%mass_of_bits                             ! retain previous value incase non-zero
+         ENDIF
+
+         ! use tmask rather than tmask_i when dealing with icebergs
+         IF( tmask(ii,ij,1) /= 0._wp ) THEN    ! Add melting to the grid and field diagnostics
+            z1_e1e2    = 1._wp / e1e2t(ii,ij) * this%mass_scaling
+            z1_dt_e1e2 = z1_dt * z1_e1e2
+            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
+            berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt    * z1_e1e2    ! kg/m2/s
+            zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s
+            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat    * z1_e1e2    ! W/m2
+            CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling,       &
+            &                          zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
+            &                          zdMv, z1_dt_e1e2 )
+         ELSE
+            WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
+            CALL icb_utl_print_berg( this, kt )
+            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
+            CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
+         ENDIF
+
+         ! Rolling
+         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
+         IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
+            zT  = zTn
+            zTn = zWn
+            zWn = zT
+         endif
+
+         ! Store the new state of iceberg (with L>W)
+         pt%mass         = zMnew
+         pt%mass_of_bits = znMbits
+         pt%thickness    = zTn
+         pt%width        = min(zWn,zLn)
+         pt%length       = max(zWn,zLn)
+
+         next=>this%next
+
+!!gm  add a test to avoid over melting ?
+
+         IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
+            CALL icb_utl_delete( first_berg, this )
+            !
+         ELSE                            ! Diagnose mass distribution on grid
+            z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling
+            CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
+            &                  this%mass_scaling, zMnew, znMbits, z1_e1e2)
+         ENDIF
+         !
+         this=>next
+         !
+      END DO
+      
+      ! now use melt and associated heat flux in ocean (or not)
+      !
+      IF(.NOT. ln_passive_mode ) THEN
+         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
+         emps(:,:) = emps(:,:) - berg_grid%floating_melt(:,:)
+!!       qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:)  !!gm heat flux not yet properly coded ==>> need it, SOLVE that!
+      ENDIF
+      !
+   END SUBROUTINE icb_thm
+
+   !!======================================================================
+END MODULE icbthm
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 3609)
@@ -0,0 +1,267 @@
+MODULE icbtrj
+
+   !!======================================================================
+   !!                       ***  MODULE  icbtrj  ***
+   !! Ocean physics:  trajectory I/O routines
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-05  (Alderson)       New module to handle trajectory output
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_trj_init          :
+   !!----------------------------------------------------------------------
+   USE par_oce        ! NEMO parameters
+   USE dom_oce        ! NEMO ocean domain
+   USE phycst         ! NEMO physical constants
+   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
+   USE in_out_manager ! NEMO IO, numout in particular
+   USE netcdf
+   !
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_trj_init    ! routine called in icbini.F90 module
+   PUBLIC   icb_trj_write   ! routine called in icbstp.F90 module
+   PUBLIC   icb_trj_sync    ! routine called in icbstp.F90 module
+   PUBLIC   icb_trj_end     ! routine called in icbstp.F90 module
+
+   INTEGER ::   num_traj
+   INTEGER ::   n_dim, m_dim
+   INTEGER ::   ntrajid
+   INTEGER ::   numberid, nstepid, nscaling_id
+   INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid
+   INTEGER ::   nuoid, nvoid, nuaid, nvaid, nuiid, nviid
+   INTEGER ::   nsshxid, nsshyid, nsstid, ncntid, nthkid
+   INTEGER ::   nthicknessid, nwidthid, nlengthid
+   INTEGER ::   nyearid, ndayid
+   INTEGER ::   nmass_of_bits_id, nheat_density_id
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_trj_init( ktend )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_trj_init  ***
+      !!
+      !! ** Purpose :   initialise iceberg trajectory output files
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )                 :: ktend
+      !
+      INTEGER                               :: iret
+      CHARACTER(len=80)                     :: cl_filename
+      TYPE(iceberg), POINTER                :: this
+      TYPE(point)  , POINTER                :: pt
+      !!----------------------------------------------------------------------
+
+      IF( lk_mpp ) THEN   ;   WRITE(cl_filename,'("trajectory_icebergs_",I6.6,"_",I4.4,".nc")') ktend, narea-1
+      ELSE                ;   WRITE(cl_filename,'("trajectory_icebergs_",I6.6         ,".nc")') ktend
+      ENDIF
+      IF ( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, icb_trj_init: creating ',TRIM(cl_filename)
+
+      iret = NF90_CREATE(TRIM(cl_filename), NF90_CLOBBER, ntrajid)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, icb_trj_init: nf_create failed')
+
+      ! Dimensions
+      iret = NF90_DEF_DIM(ntrajid, 'n', NF90_UNLIMITED, n_dim)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, icb_trj_init: nf_def_dim n failed')
+      iret = NF90_DEF_DIM(ntrajid, 'k', nkounts, m_dim)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, icb_trj_init: nf_def_dim k failed')
+
+      ! Variables
+      iret = NF90_DEF_VAR(ntrajid, 'iceberg_number', NF90_INT, (/m_dim,n_dim/), numberid)
+      iret = NF90_DEF_VAR(ntrajid, 'timestep', NF90_INT, n_dim, nstepid)
+      iret = NF90_DEF_VAR(ntrajid, 'mass_scaling', NF90_DOUBLE, n_dim, nscaling_id)
+      iret = NF90_DEF_VAR(ntrajid, 'lon', NF90_DOUBLE, n_dim, nlonid)
+      iret = NF90_DEF_VAR(ntrajid, 'lat', NF90_DOUBLE, n_dim, nlatid)
+      iret = NF90_DEF_VAR(ntrajid, 'xi', NF90_DOUBLE, n_dim, nxid)
+      iret = NF90_DEF_VAR(ntrajid, 'yj', NF90_DOUBLE, n_dim, nyid)
+      iret = NF90_DEF_VAR(ntrajid, 'uvel', NF90_DOUBLE, n_dim, nuvelid)
+      iret = NF90_DEF_VAR(ntrajid, 'vvel', NF90_DOUBLE, n_dim, nvvelid)
+      iret = NF90_DEF_VAR(ntrajid, 'uto', NF90_DOUBLE, n_dim, nuoid)
+      iret = NF90_DEF_VAR(ntrajid, 'vto', NF90_DOUBLE, n_dim, nvoid)
+      iret = NF90_DEF_VAR(ntrajid, 'uta', NF90_DOUBLE, n_dim, nuaid)
+      iret = NF90_DEF_VAR(ntrajid, 'vta', NF90_DOUBLE, n_dim, nvaid)
+      iret = NF90_DEF_VAR(ntrajid, 'uti', NF90_DOUBLE, n_dim, nuiid)
+      iret = NF90_DEF_VAR(ntrajid, 'vti', NF90_DOUBLE, n_dim, nviid)
+      iret = NF90_DEF_VAR(ntrajid, 'ssh_x', NF90_DOUBLE, n_dim, nsshyid)
+      iret = NF90_DEF_VAR(ntrajid, 'ssh_y', NF90_DOUBLE, n_dim, nsshxid)
+      iret = NF90_DEF_VAR(ntrajid, 'sst', NF90_DOUBLE, n_dim, nsstid)
+      iret = NF90_DEF_VAR(ntrajid, 'icnt', NF90_DOUBLE, n_dim, ncntid)
+      iret = NF90_DEF_VAR(ntrajid, 'ithk', NF90_DOUBLE, n_dim, nthkid)
+      iret = NF90_DEF_VAR(ntrajid, 'mass', NF90_DOUBLE, n_dim, nmassid)
+      iret = NF90_DEF_VAR(ntrajid, 'thickness', NF90_DOUBLE, n_dim, nthicknessid)
+      iret = NF90_DEF_VAR(ntrajid, 'width', NF90_DOUBLE, n_dim, nwidthid)
+      iret = NF90_DEF_VAR(ntrajid, 'length', NF90_DOUBLE, n_dim, nlengthid)
+      iret = NF90_DEF_VAR(ntrajid, 'year', NF90_INT, n_dim, nyearid)
+      iret = NF90_DEF_VAR(ntrajid, 'day', NF90_DOUBLE, n_dim, ndayid)
+      iret = NF90_DEF_VAR(ntrajid, 'mass_of_bits', NF90_DOUBLE, n_dim, nmass_of_bits_id)
+      iret = NF90_DEF_VAR(ntrajid, 'heat_density', NF90_DOUBLE, n_dim, nheat_density_id)
+
+      ! Attributes
+      iret = NF90_PUT_ATT(ntrajid, numberid, 'long_name', 'iceberg number on this processor')
+      iret = NF90_PUT_ATT(ntrajid, numberid, 'units', 'count')
+      iret = NF90_PUT_ATT(ntrajid, nstepid, 'long_name', 'timestep number kt')
+      iret = NF90_PUT_ATT(ntrajid, nstepid, 'units', 'count')
+      iret = NF90_PUT_ATT(ntrajid, nlonid, 'long_name', 'longitude')
+      iret = NF90_PUT_ATT(ntrajid, nlonid, 'units', 'degrees_E')
+      iret = NF90_PUT_ATT(ntrajid, nlatid, 'long_name', 'latitude')
+      iret = NF90_PUT_ATT(ntrajid, nlatid, 'units', 'degrees_N')
+      iret = NF90_PUT_ATT(ntrajid, nxid, 'long_name', 'x grid box position')
+      iret = NF90_PUT_ATT(ntrajid, nxid, 'units', 'fractional')
+      iret = NF90_PUT_ATT(ntrajid, nyid, 'long_name', 'y grid box position')
+      iret = NF90_PUT_ATT(ntrajid, nyid, 'units', 'fractional')
+      iret = NF90_PUT_ATT(ntrajid, nuvelid, 'long_name', 'zonal velocity')
+      iret = NF90_PUT_ATT(ntrajid, nuvelid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nvvelid, 'long_name', 'meridional velocity')
+      iret = NF90_PUT_ATT(ntrajid, nvvelid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nuoid, 'long_name', 'ocean u component')
+      iret = NF90_PUT_ATT(ntrajid, nuoid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nvoid, 'long_name', 'ocean v component')
+      iret = NF90_PUT_ATT(ntrajid, nvoid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nuaid, 'long_name', 'atmosphere u component')
+      iret = NF90_PUT_ATT(ntrajid, nuaid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nvaid, 'long_name', 'atmosphere v component')
+      iret = NF90_PUT_ATT(ntrajid, nvaid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nuiid, 'long_name', 'sea ice u component')
+      iret = NF90_PUT_ATT(ntrajid, nuiid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nviid, 'long_name', 'sea ice v component')
+      iret = NF90_PUT_ATT(ntrajid, nviid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(ntrajid, nsshxid, 'long_name', 'sea surface height gradient from x points')
+      iret = NF90_PUT_ATT(ntrajid, nsshxid, 'units', 'm/m')
+      iret = NF90_PUT_ATT(ntrajid, nsshyid, 'long_name', 'sea surface height gradient from y points')
+      iret = NF90_PUT_ATT(ntrajid, nsshyid, 'units', 'm/m')
+      iret = NF90_PUT_ATT(ntrajid, nsstid, 'long_name', 'sea surface temperature')
+      iret = NF90_PUT_ATT(ntrajid, nsstid, 'units', 'degC')
+      iret = NF90_PUT_ATT(ntrajid, ncntid, 'long_name', 'sea ice concentration')
+      iret = NF90_PUT_ATT(ntrajid, ncntid, 'units', 'degC')
+      iret = NF90_PUT_ATT(ntrajid, nthkid, 'long_name', 'sea ice thickness')
+      iret = NF90_PUT_ATT(ntrajid, nthkid, 'units', 'm')
+      iret = NF90_PUT_ATT(ntrajid, nmassid, 'long_name', 'mass')
+      iret = NF90_PUT_ATT(ntrajid, nmassid, 'units', 'kg')
+      iret = NF90_PUT_ATT(ntrajid, nthicknessid, 'long_name', 'thickness')
+      iret = NF90_PUT_ATT(ntrajid, nthicknessid, 'units', 'm')
+      iret = NF90_PUT_ATT(ntrajid, nwidthid, 'long_name', 'width')
+      iret = NF90_PUT_ATT(ntrajid, nwidthid, 'units', 'm')
+      iret = NF90_PUT_ATT(ntrajid, nlengthid, 'long_name', 'length')
+      iret = NF90_PUT_ATT(ntrajid, nlengthid, 'units', 'm')
+      iret = NF90_PUT_ATT(ntrajid, nyearid, 'long_name', 'calendar year')
+      iret = NF90_PUT_ATT(ntrajid, nyearid, 'units', 'years')
+      iret = NF90_PUT_ATT(ntrajid, ndayid, 'long_name', 'day of year')
+      iret = NF90_PUT_ATT(ntrajid, ndayid, 'units', 'days')
+      iret = NF90_PUT_ATT(ntrajid, nscaling_id, 'long_name', 'scaling factor for mass of berg')
+      iret = NF90_PUT_ATT(ntrajid, nscaling_id, 'units', 'none')
+      iret = NF90_PUT_ATT(ntrajid, nmass_of_bits_id, 'long_name', 'mass of bergy bits')
+      iret = NF90_PUT_ATT(ntrajid, nmass_of_bits_id, 'units', 'kg')
+      iret = NF90_PUT_ATT(ntrajid, nheat_density_id, 'long_name', 'heat density')
+      iret = NF90_PUT_ATT(ntrajid, nheat_density_id, 'units', 'J/kg')
+
+      ! End define mode
+      iret = NF90_ENDDEF(ntrajid)
+      !
+   END SUBROUTINE icb_trj_init
+
+
+   SUBROUTINE icb_trj_write( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_trj_write  ***
+      !!
+      !! ** Purpose :   write out iceberg trajectories
+      !!
+      !! ** Method  : - for the moment write out each snapshot of positions later
+      !!                can rewrite so that it is buffered and written out more efficiently
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt
+      !
+      INTEGER                ::   iret, jn
+      CHARACTER(len=80)      ::   cl_filename
+      TYPE(iceberg), POINTER ::   this
+      TYPE(point)  , POINTER ::   pt
+      !!----------------------------------------------------------------------
+
+      ! Write variables
+      ! sga - just write out the current point of the trajectory
+
+      this => first_berg
+      jn = num_traj
+      DO WHILE (ASSOCIATED(this))
+         pt => this%current_point
+         jn=jn+1
+
+         iret = NF90_PUT_VAR(ntrajid, numberid, this%number, (/1,jn/), (/nkounts,1/) )
+         iret = NF90_PUT_VAR(ntrajid, nstepid, kt, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nscaling_id, this%mass_scaling, (/ jn /) )
+
+         iret = NF90_PUT_VAR(ntrajid, nlonid, pt%lon, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nlatid, pt%lat, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nxid, pt%xi, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nyid, pt%yj, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nuvelid, pt%uvel, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nvvelid, pt%vvel, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nuoid, pt%uo, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nvoid, pt%vo, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nuaid, pt%ua, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nvaid, pt%va, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nuiid, pt%ui, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nviid, pt%vi, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nsshxid, pt%ssh_x, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nsshyid, pt%ssh_y, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nsstid, pt%sst, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, ncntid, pt%cn, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nthkid, pt%hi, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nmassid, pt%mass, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nthicknessid, pt%thickness, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nwidthid, pt%width, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nlengthid, pt%length, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nyearid, pt%year, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, ndayid, pt%day, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nmass_of_bits_id, pt%mass_of_bits, (/ jn /) )
+         iret = NF90_PUT_VAR(ntrajid, nheat_density_id, pt%heat_density, (/ jn /) )
+
+         this=>this%next
+      END DO
+      IF( lwp .and. nn_verbose_level > 0 ) WRITE(numout,*) 'trajectory write to frame ', jn
+      num_traj = jn
+      !
+   END SUBROUTINE icb_trj_write
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_trj_sync()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_trj_sync  ***
+      !!
+      !! ** Purpose :   
+      !!----------------------------------------------------------------------
+      INTEGER                               :: iret
+      !!----------------------------------------------------------------------
+      ! flush to file
+      iret = NF90_SYNC(ntrajid)
+      IF(iret /= NF90_NOERR)   CALL ctl_stop( 'icebergs, icb_trj_sync: nf_sync failed' )
+      !
+   END SUBROUTINE icb_trj_sync
+
+
+   SUBROUTINE icb_trj_end()
+      ! Local variables
+      INTEGER                               :: iret
+      !!----------------------------------------------------------------------
+      ! Finish up
+      iret = NF90_CLOSE(ntrajid)
+      IF (iret /= NF90_NOERR)   CALL ctl_stop( 'icebergs, icb_trj_end: nf_close failed' )
+      !
+   END SUBROUTINE icb_trj_end
+
+   !!-------------------------------------------------------------------------
+
+END MODULE icbtrj
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 3609)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 3609)
@@ -0,0 +1,754 @@
+MODULE icbutl
+
+   !!======================================================================
+   !!                       ***  MODULE  icbutl  ***
+   !! Icebergs:  various iceberg utility routines
+   !!======================================================================
+   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
+   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
+   !!            -    !                            Removal of mapping from another grid
+   !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_utl_interp   :
+   !!   icb_utl_bilin    :
+   !!   icb_utl_bilin_e  :
+   !!----------------------------------------------------------------------
+   USE par_oce                             ! ocean parameters
+   USE dom_oce                             ! ocean domain
+   USE in_out_manager                      ! IO parameters
+   USE lbclnk                              ! lateral boundary condition
+   USE lib_mpp                             ! MPI code and lk_mpp in particular
+   USE icb_oce                             ! define iceberg arrays
+   USE sbc_oce                             ! ocean surface boundary conditions
+#if defined key_lim2
+   USE ice_2,         ONLY: u_ice, v_ice   ! LIM-2 ice velocities  (CAUTION in C-grid do not use key_vp option)
+   USE ice_2,         ONLY: hicif          ! LIM-2 ice thickness
+#elif defined key_lim3
+   USE ice,           ONLY: u_ice, v_ice   ! LIM-3 variables  (always in C-grid)
+                                           ! gm  LIM3 case the mean ice thickness (i.e. averaged over categories)
+                                           ! gm            has to be computed somewhere in the ice and accessed here
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   icb_utl_copy          ! routine called in icbstp module
+   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
+   PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules
+   PUBLIC   icb_utl_bilin_x       ! routine called in icbdyn module
+   PUBLIC   icb_utl_add           ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
+   PUBLIC   icb_utl_delete        ! routine called in icblbc, icbthm modules
+   PUBLIC   icb_utl_destroy       ! routine called in icbstp module
+   PUBLIC   icb_utl_track         ! routine not currently used, retain just in case
+   PUBLIC   icb_utl_print_berg    ! routine called in icbthm module
+   PUBLIC   icb_utl_print         ! routine called in icbini, icbstp module
+   PUBLIC   icb_utl_count         ! routine called in icbdia, icbini, icblbc, icbrst modules
+   PUBLIC   icb_utl_incr          ! routine called in icbini, icbclv modules
+   PUBLIC   icb_utl_yearday       ! routine called in icbclv, icbstp module
+   PUBLIC   icb_utl_mass          ! routine called in icbdia module
+   PUBLIC   icb_utl_heat          ! routine called in icbdia module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!-------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE icb_utl_copy()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_utl_copy  ***
+      !!
+      !! ** Purpose :   iceberg initialization.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+
+      ! copy nemo forcing arrays into iceberg versions with extra halo
+      ! only necessary for variables not on T points
+      ! and ssh which is used to calculate gradients
+
+      uo_e(:,:) = 0._wp ;   uo_e(1:jpi, 1:jpj) = ssu_m(:,:)
+      vo_e(:,:) = 0._wp ;   vo_e(1:jpi, 1:jpj) = ssv_m(:,:)
+      ff_e(:,:) = 0._wp ;   ff_e(1:jpi, 1:jpj) = ff   (:,:)
+      ua_e(:,:) = 0._wp ;   ua_e(1:jpi, 1:jpj) = utau (:,:)
+      va_e(:,:) = 0._wp ;   va_e(1:jpi, 1:jpj) = vtau (:,:)
+
+      CALL lbc_lnk_e( uo_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_e( vo_e, 'V', -1._wp, 1, 1 )
+      CALL lbc_lnk_e( ff_e, 'F', +1._wp, 1, 1 )
+      CALL lbc_lnk_e( ua_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_e( va_e, 'V', -1._wp, 1, 1 )
+
+#if defined key_lim2 || defined key_lim3
+      ui_e(:,:) = 0._wp ;   ui_e(1:jpi, 1:jpj) = u_ice(:,:)
+      vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:)
+
+      CALL lbc_lnk_e( ui_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_e( vi_e, 'V', -1._wp, 1, 1 )
+#endif
+
+      !! special for ssh which is used to calculate slope
+      !! so fudge some numbers all the way around the boundary
+
+      ssh_e(:,:) = 0._wp ;   ssh_e(1:jpi, 1:jpj) = ssh_m(:,:)
+      ssh_e(0    ,    :) = ssh_e(1  ,  :)
+      ssh_e(jpi+1,    :) = ssh_e(jpi,  :)
+      ssh_e(:    ,    0) = ssh_e(:  ,  1)
+      ssh_e(:    ,jpj+1) = ssh_e(:  ,jpj)
+      ssh_e(0,0)         = ssh_e(1,1)
+      ssh_e(jpi+1,0)     = ssh_e(jpi,1)
+      ssh_e(0,jpj+1)     = ssh_e(1,jpj)
+      ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj)
+      CALL lbc_lnk_e( ssh_e, 'T', +1._wp, 1, 1 )
+      !
+   END SUBROUTINE icb_utl_copy
+
+
+   SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   &
+      &                       pj, pe2, pvo, pvi, pva, pssh_j,   &
+      &                       psst, pcn, phi, pff            )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_utl_interp  ***
+      !!
+      !! ** Purpose :   interpolation
+      !!
+      !! ** Method  : - interpolate from various ocean arrays onto iceberg position
+      !!
+      !!       !!gm  CAUTION here I do not care of the slip/no-slip conditions
+      !!             this can be done later (not that easy to do...)
+      !!             right now, U is 0 in land so that the coastal value of velocity parallel to the coast
+      !!             is half the off shore value, wile the normal-to-the-coast value is zero.
+      !!             This is OK as a starting point.
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential
+      REAL(wp), INTENT(  out) ::   pe1, pe2                       ! i- and j scale factors
+      REAL(wp), INTENT(  out) ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds
+      REAL(wp), INTENT(  out) ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
+      REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis
+      !
+      REAL(wp) ::   zcd, zmod       ! local scalars
+      !!----------------------------------------------------------------------
+
+      pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )         ! scale factors
+      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+      !
+      puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U' )             ! ocean velocities
+      pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V' )
+      psst = icb_utl_bilin( sst_m, pi, pj, 'T' )            ! SST
+      pcn  = icb_utl_bilin( fr_i , pi, pj, 'T' )            ! ice concentration
+      pff  = icb_utl_bilin_h( ff_e , pi, pj, 'F' )            ! Coriolis parameter
+      !
+      pua  = icb_utl_bilin_h( ua_e , pi, pj, 'U' )            ! 10m wind
+      pva  = icb_utl_bilin_h( va_e , pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed
+      zcd  = 1.22_wp * 1.5e-3_wp                                  ! air density * drag coefficient
+      zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  )
+      pua  = pua * zmod                                           ! note: stress module=0 necessarly implies ua=va=0
+      pva  = pva * zmod
+
+#if defined key_lim2 || defined key_lim3
+      pui = icb_utl_bilin_h( ui_e, pi, pj, 'U' )              ! sea-ice velocities
+      pvi = icb_utl_bilin_h( vi_e, pi, pj, 'V' )
+      phi = icb_utl_bilin(hicif, pi, pj, 'T' )              ! ice thickness
+#else
+      pui = 0._wp
+      pvi = 0._wp
+      phi = 0._wp
+#endif
+
+      ! Estimate SSH gradient in i- and j-direction (centred evaluation)
+      pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   &
+          &      icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 )
+      pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   &
+          &      icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 )
+      !
+   END SUBROUTINE icb_utl_interp
+
+
+   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
+      !!                this version deals with extra halo points
+      !!
+      !!       !!gm  CAUTION an optional argument should be added to handle
+      !!             the slip/no-slip conditions  ==>>> to be done later
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
+      CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
+      !
+      INTEGER  ::   ii, ij   ! local integer
+      REAL(wp) ::   zi, zj   ! local real
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( cd_type )
+         CASE ( 'T' )
+            ! note that here there is no +0.5 added
+            ! since we're looking for four T points containing quadrant we're in of 
+            ! current T cell
+            ii = INT( pi     )
+            ij = INT( pj      )    ! T-point
+            zi = pi - REAL(ii,wp)
+            zj = pj - REAL(ij,wp)
+         CASE ( 'U' )
+            ii = INT( pi-0.5 )
+            ij = INT( pj      )    ! U-point
+            zi = pi - 0.5 - REAL(ii,wp)
+            zj = pj - REAL(ij,wp)
+         CASE ( 'V' )
+            ii = INT( pi     )
+            ij = INT( pj -0.5 )    ! V-point
+            zi = pi - REAL(ii,wp)
+            zj = pj - 0.5 - REAL(ij,wp)
+         CASE ( 'F' )
+            ii = INT( pi-0.5 )
+            ij = INT( pj -0.5 )    ! F-point
+            zi = pi - 0.5 - REAL(ii,wp)
+            zj = pj - 0.5 - REAL(ij,wp)
+      END SELECT
+      !
+      ! find position in this processor
+      ii = mi1( ii )
+      ij = mj1( ij )
+      !
+      icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
+         &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
+      !
+   END FUNCTION icb_utl_bilin_h
+
+
+   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
+      !!
+      !!       !!gm  CAUTION an optional argument should be added to handle
+      !!             the slip/no-slip conditions  ==>>> to be done later
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
+      CHARACTER(len=1)            , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
+      !
+      INTEGER  ::   ii, ij   ! local integer
+      REAL(wp) ::   zi, zj   ! local real
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( cd_type )
+         CASE ( 'T' )
+            ! note that here there is no +0.5 added
+            ! since we're looking for four T points containing quadrant we're in of 
+            ! current T cell
+            ii = INT( pi     )
+            ij = INT( pj      )    ! T-point
+            zi = pi - REAL(ii,wp)
+            zj = pj - REAL(ij,wp)
+         CASE ( 'U' )
+            ii = INT( pi-0.5 )
+            ij = INT( pj      )    ! U-point
+            zi = pi - 0.5 - REAL(ii,wp)
+            zj = pj - REAL(ij,wp)
+         CASE ( 'V' )
+            ii = INT( pi     )
+            ij = INT( pj -0.5 )    ! V-point
+            zi = pi - REAL(ii,wp)
+            zj = pj - 0.5 - REAL(ij,wp)
+         CASE ( 'F' )
+            ii = INT( pi-0.5 )
+            ij = INT( pj -0.5 )    ! F-point
+            zi = pi - 0.5 - REAL(ii,wp)
+            zj = pj - 0.5 - REAL(ij,wp)
+      END SELECT
+      !
+      ! find position in this processor
+      ii = mi1( ii )
+      ij = mj1( ij )
+      !
+      icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
+         &          + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
+      !
+   END FUNCTION icb_utl_bilin
+
+
+   REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin_x  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
+      !!                Special case for interpolating longitude
+      !!
+      !!       !!gm  CAUTION an optional argument should be added to handle
+      !!             the slip/no-slip conditions  ==>>> to be done later
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
+      !
+      INTEGER                                  ::   ii, ij   ! local integer
+      REAL(wp)                                 ::   zi, zj   ! local real
+      REAL(wp)                                 ::   zret     ! local real
+      REAL(wp), DIMENSION(4)                   ::   z4
+      !!----------------------------------------------------------------------
+      !
+      ! note that here there is no +0.5 added
+      ! since we're looking for four T points containing quadrant we're in of 
+      ! current T cell
+      ii = INT( pi     )
+      ij = INT( pj      )    ! T-point
+      zi = pi - REAL(ii,wp)
+      zj = pj - REAL(ij,wp)
+      !
+      ! find position in this processor          !!gm use here mig, mjg arrays
+      ii = mi1( ii )
+      ij = mj1( ij )
+      z4(1) = pfld(ii  ,ij  )
+      z4(2) = pfld(ii+1,ij  )
+      z4(3) = pfld(ii  ,ij+1)
+      z4(4) = pfld(ii+1,ij+1)
+      IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
+         WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
+      ENDIF
+      !
+      zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
+      IF( zret > 180._wp ) zret = zret - 360._wp
+      icb_utl_bilin_x = zret
+      !
+   END FUNCTION icb_utl_bilin_x
+
+
+   REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION dom_init  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location of horizontal scale factor
+      !! ** Method  :   interpolation done using the 4 nearest grid points among
+      !!                t-, u-, v-, and f-points.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
+      REAL(wp)                , INTENT(in) ::   pi, pj               ! targeted coordinates in (i,j) referential
+      !
+      INTEGER  ::   ii, ij, icase   ! local integer
+      !
+      ! weights corresponding to corner points of a T cell quadrant
+      REAL(wp) ::   zi, zj          ! local real
+      !
+      ! values at corner points of a T cell quadrant
+      ! 00 = bottom left, 10 = bottom right, 01 = top left, 11 = top right
+      REAL(wp) ::   ze00, ze10, ze01, ze11
+      !!----------------------------------------------------------------------
+      !
+      ii = INT( pi )   ;   ij = INT( pj )    ! left bottom T-point (i,j) indices
+
+      ! fractional box spacing
+      ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell
+      ! 0.5 <= zi < 1  , 0   <= zj < 0.5   -->  NE quadrant
+      ! 0   <= zi < 0.5, 0.5 <= zj < 1     -->  SE quadrant
+      ! 0.5 <= zi < 1  , 0.5 <= zj < 1     -->  SW quadrant
+
+      zi = pi - REAL(ii,wp)          !!gm use here mig, mjg arrays
+      zj = pj - REAL(ij,wp)
+
+      ! find position in this processor
+      ii = mi1( ii )
+      ij = mj1( ij )
+
+      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
+         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NE quadrant
+            !                                                      !             i=I       i=I+1/2
+            ze01 = pev(ii  ,ij  )   ;   ze11 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
+            ze00 = pet(ii  ,ij  )   ;   ze10 = peu(ii  ,ij  )      !   j=J        T ------- U
+            zi = 2._wp * zi
+            zj = 2._wp * zj
+         ELSE                                                      !  SE quadrant
+            !                                                                    !             i=I       i=I+1/2
+            ze01 = pet(ii  ,ij+1)   ;   ze11 = peu(ii  ,ij+1)      !   j=J+1      T ------- U
+            ze00 = pev(ii  ,ij  )   ;   ze10 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
+            zi = 2._wp *  zi
+            zj = 2._wp * (zj-0.5_wp)
+         ENDIF
+      ELSE
+         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NW quadrant
+            !                                                                    !             i=I       i=I+1/2
+            ze01 = pef(ii  ,ij  )   ;   ze11 = pev(ii+1,ij)        !   j=J+1/2    F ------- V
+            ze00 = peu(ii  ,ij  )   ;   ze10 = pet(ii+1,ij)        !   j=J        U ------- T
+            zi = 2._wp * (zi-0.5_wp)
+            zj = 2._wp *  zj
+         ELSE                                                      !  SW quadrant
+            !                                                                    !             i=I+1/2   i=I+1
+            ze01 = peu(ii  ,ij+1)   ;   ze11 = pet(ii+1,ij+1)      !   j=J+1      U ------- T
+            ze00 = pef(ii  ,ij  )   ;   ze10 = pev(ii+1,ij  )      !   j=J+1/2    F ------- V
+            zi = 2._wp * (zi-0.5_wp)
+            zj = 2._wp * (zj-0.5_wp)
+         ENDIF
+      ENDIF
+      !
+      icb_utl_bilin_e = ( ze01 * (1.-zi) + ze11 * zi ) *     zj    &
+         &            + ( ze00 * (1.-zi) + ze10 * zi ) * (1.-zj)
+      !
+   END FUNCTION icb_utl_bilin_e
+
+
+   SUBROUTINE icb_utl_add( bergvals, ptvals )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_add           ***
+      !!
+      !! ** Purpose :   add a new berg to the iceberg list
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), INTENT(in)           ::   bergvals
+      TYPE(point)  , INTENT(in)           ::   ptvals
+      !
+      TYPE(iceberg), POINTER ::   new => NULL()
+      !!----------------------------------------------------------------------
+      !
+      new => NULL()
+      CALL icb_utl_create( new, bergvals, ptvals )
+      CALL icb_utl_insert( new )
+      new => NULL()     ! Clear new
+      !
+   END SUBROUTINE icb_utl_add         
+
+
+   SUBROUTINE icb_utl_create( berg, bergvals, ptvals )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_create  ***
+      !!
+      !! ** Purpose :   add a new berg to the iceberg list
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), INTENT(in) ::   bergvals
+      TYPE(point)  , INTENT(in) ::   ptvals
+      TYPE(iceberg), POINTER    ::   berg
+      !
+      TYPE(point)  , POINTER    ::   pt
+      INTEGER                   ::   istat
+      !!----------------------------------------------------------------------
+      !
+      IF( ASSOCIATED(berg) )   CALL ctl_stop( 'icebergs, icb_utl_create: berg already associated' )
+      ALLOCATE(berg, STAT=istat)
+      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
+      berg%number(:) = bergvals%number(:)
+      berg%mass_scaling = bergvals%mass_scaling
+      berg%prev => NULL()
+      berg%next => NULL()
+      !
+      ALLOCATE(pt, STAT=istat)
+      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
+      pt = ptvals
+      berg%current_point => pt
+      !
+   END SUBROUTINE icb_utl_create
+
+
+   SUBROUTINE icb_utl_insert( newberg )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_insert  ***
+      !!
+      !! ** Purpose :   add a new berg to the iceberg list
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER  ::   newberg
+      !
+      TYPE(iceberg), POINTER  ::   this, prev, last
+      !!----------------------------------------------------------------------
+      !
+      IF( ASSOCIATED( first_berg ) ) THEN
+         last => first_berg
+         DO WHILE (ASSOCIATED(last%next))
+            last => last%next
+         ENDDO
+         newberg%prev => last
+         last%next    => newberg
+         last         => newberg
+      ELSE                       ! list is empty so create it
+         first_berg => newberg
+      ENDIF
+      !
+   END SUBROUTINE icb_utl_insert
+
+
+   REAL(wp) FUNCTION icb_utl_yearday(kmon, kday, khr, kmin, ksec)
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION icb_utl_yearday  ***
+      !!
+      !! ** Purpose :   
+      !!
+      ! sga - improved but still only applies to 365 day year, need to do this properly
+      !
+      !!gm  all these info are already known in daymod, no???
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)     :: kmon, kday, khr, kmin, ksec
+      ! 
+      INTEGER, DIMENSION(12)  :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
+      !!----------------------------------------------------------------------
+      !
+      icb_utl_yearday = REAL( SUM( imonths(1:kmon) ), wp )
+      icb_utl_yearday = icb_utl_yearday + REAL(kday-1,wp) + (REAL(khr,wp) + (REAL(kmin,wp) + REAL(ksec,wp)/60.)/60.)/24.
+      !
+   END FUNCTION icb_utl_yearday
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_utl_delete( first, berg )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_delete  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER :: first, berg
+      !!----------------------------------------------------------------------
+      ! Connect neighbors to each other
+      IF ( ASSOCIATED(berg%prev) ) THEN
+        berg%prev%next => berg%next
+      ELSE
+        first => berg%next
+      ENDIF
+      IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
+      !
+      CALL icb_utl_destroy(berg)
+      !
+   END SUBROUTINE icb_utl_delete
+
+
+   SUBROUTINE icb_utl_destroy( berg )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_destroy  ***
+      !!
+      !! ** Purpose :   remove a single iceberg instance
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER :: berg
+      !!----------------------------------------------------------------------
+      !
+      ! Remove any points
+      IF( ASSOCIATED( berg%current_point ) )   DEALLOCATE( berg%current_point )
+      !
+      DEALLOCATE(berg)
+      !
+   END SUBROUTINE icb_utl_destroy
+
+
+   SUBROUTINE icb_utl_track( knum, cd_label, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_track  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(nkounts)    :: knum       ! iceberg number
+      CHARACTER(len=*)               :: cd_label   ! 
+      INTEGER                        :: kt         ! timestep number
+      ! 
+      TYPE(iceberg), POINTER         :: this
+      LOGICAL                        :: match
+      INTEGER                        :: k
+      !!----------------------------------------------------------------------
+      !
+      this => first_berg
+      DO WHILE( ASSOCIATED(this) )
+         match = .TRUE.
+         DO k = 1, nkounts
+            IF( this%number(k) /= knum(k) ) match = .FALSE.
+         END DO
+         IF( match )   CALL icb_utl_print_berg(this, kt)
+         this => this%next
+      END DO
+      !
+   END SUBROUTINE icb_utl_track
+
+
+   SUBROUTINE icb_utl_print_berg( berg, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_print_berg  ***
+      !!
+      !! ** Purpose :   print one
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER :: berg
+      TYPE(point)  , POINTER :: pt
+      INTEGER                :: kt      ! timestep number
+      !!----------------------------------------------------------------------
+      !
+      pt => berg%current_point
+      WRITE(numicb, 9200) kt, berg%number(1), &
+                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
+                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
+      CALL flush( numicb )
+ 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
+      !
+   END SUBROUTINE icb_utl_print_berg
+
+
+   SUBROUTINE icb_utl_print( cd_label, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_print  ***
+      !!
+      !! ** Purpose :   print many
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*)       :: cd_label
+      INTEGER                :: kt             ! timestep number
+      ! 
+      INTEGER                :: ibergs, inbergs
+      TYPE(iceberg), POINTER :: this
+      !!----------------------------------------------------------------------
+      !
+      this => first_berg
+      IF( ASSOCIATED(this) ) THEN
+         WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
+         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
+            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
+      ENDIF
+      DO WHILE( ASSOCIATED(this) )
+        CALL icb_utl_print_berg(this, kt)
+        this => this%next
+      END DO
+      ibergs = icb_utl_count()
+      inbergs = ibergs
+      IF( lk_mpp )   CALL mpp_sum(inbergs)
+      IF( ibergs > 0 )   WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)')   &
+         &                                  cd_label, ibergs, inbergs, narea
+      !
+   END SUBROUTINE icb_utl_print
+
+
+   SUBROUTINE icb_utl_incr()
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE icb_utl_incr  ***
+      !!
+      !! ** Purpose :   
+      !!
+      ! Small routine for coping with very large integer values labelling icebergs
+      ! num_bergs is a array of integers
+      ! the first member is incremented in steps of jpnij starting from narea
+      ! this means each iceberg is labelled with a unique number
+      ! when this gets to the maximum allowed integer the second and subsequent members are 
+      ! used to count how many times the member before cycles
+      !!----------------------------------------------------------------------
+      INTEGER ::   ii, ibig
+      !!----------------------------------------------------------------------
+
+      ibig = HUGE(num_bergs(1))
+      IF( ibig-jpnij < num_bergs(1) ) THEN
+         num_bergs(1) = narea
+         DO ii = 2,nkounts
+            IF( num_bergs(ii) == ibig ) THEN
+               num_bergs(ii) = 0
+               IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
+            ELSE
+               num_bergs(ii) = num_bergs(ii) + 1
+               EXIT
+            ENDIF
+         END DO
+      ELSE
+         num_bergs(1) = num_bergs(1) + jpnij
+      ENDIF
+      !
+   END SUBROUTINE icb_utl_incr
+
+
+   INTEGER FUNCTION icb_utl_count()
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION icb_utl_count  ***
+      !!
+      !! ** Purpose :   
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER :: this
+      !!----------------------------------------------------------------------
+      !
+      icb_utl_count = 0
+      this => first_berg
+      DO WHILE( ASSOCIATED(this) )
+         icb_utl_count = icb_utl_count+1
+         this => this%next
+      END DO
+      !
+   END FUNCTION icb_utl_count
+
+
+   REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION icb_utl_mass  ***
+      !!
+      !! ** Purpose :   compute the mass all iceberg, all berg bits or all bergs.
+      !!----------------------------------------------------------------------
+      TYPE(iceberg)      , POINTER  ::   first
+      TYPE(point)        , POINTER  ::   pt
+      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
+      !
+      TYPE(iceberg), POINTER ::   this
+      !!----------------------------------------------------------------------
+      icb_utl_mass = 0._wp
+      this => first
+      !
+      IF( PRESENT( justbergs  ) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_mass = icb_utl_mass + pt%mass         * this%mass_scaling
+            this => this%next
+         END DO
+      ELSEIF( PRESENT(justbits) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
+            this => this%next
+         END DO
+      ELSE
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
+            this => this%next
+         END DO
+      ENDIF
+      !
+   END FUNCTION icb_utl_mass
+
+
+   REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION icb_utl_heat  ***
+      !!
+      !! ** Purpose :   compute the heat in all iceberg, all bergies or all bergs.
+      !!----------------------------------------------------------------------
+      TYPE(iceberg)      , POINTER  ::   first
+      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
+      !
+      TYPE(iceberg)      , POINTER  ::   this
+      TYPE(point)        , POINTER  ::   pt
+      !!----------------------------------------------------------------------
+      icb_utl_heat = 0._wp
+      this => first
+      !
+      IF( PRESENT( justbergs  ) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_heat = icb_utl_heat + pt%mass         * this%mass_scaling * pt%heat_density
+            this => this%next
+         END DO
+      ELSEIF( PRESENT(justbits) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
+            this => this%next
+         END DO
+      ELSE
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
+            this => this%next
+         END DO
+      ENDIF
+      !
+   END FUNCTION icb_utl_heat
+
+   !!======================================================================
+END MODULE icbutl
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 3609)
@@ -55,5 +55,5 @@
 
    INTERFACE lbc_lnk_e
-      MODULE PROCEDURE lbc_lnk_2d
+      MODULE PROCEDURE lbc_lnk_2d_e
    END INTERFACE
 
@@ -270,4 +270,28 @@
    END SUBROUTINE lbc_lnk_2d
 
+   SUBROUTINE lbc_lnk_2d_e( pt2d, cd_type, psgn, jpri, jprj )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE lbc_lnk_2d  ***
+      !!
+      !! ** Purpose :   set lateral boundary conditions on a 2D array (non mpp case)
+      !!                special dummy routine to allow for use of halo indexing in mpp case
+      !!
+      !! ** Method  :   psign = -1 :    change the sign across the north fold
+      !!                      =  1 : no change of the sign across the north fold
+      !!                      =  0 : no change of the sign across the north fold and
+      !!                             strict positivity preserved: use inner row/column
+      !!                             for closed boundaries.
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1)            , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt2d      ! 2D array on which the lbc is applied
+      REAL(wp)                    , INTENT(in   )           ::   psgn      ! control of the sign 
+      INTEGER                     , INTENT(in   )           ::   jpri      ! size of extra halo (not needed in non-mpp)
+      INTEGER                     , INTENT(in   )           ::   jprj      ! size of extra halo (not needed in non-mpp)
+      !!----------------------------------------------------------------------
+
+      CALL lbc_lnk_2d( pt2d, cd_type, psgn )
+      !    
+   END SUBROUTINE lbc_lnk_2d_e
+
 #endif
 
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 3609)
@@ -67,4 +67,5 @@
    PUBLIC   mppobc, mpp_ini_ice, mpp_ini_znl
    PUBLIC   mppsize
+   PUBLIC   mppsend, mpprecv                          ! needed by ICB routines
    PUBLIC   lib_mpp_alloc   ! Called in nemogcm.F90
 
@@ -143,7 +144,7 @@
 
    ! Type of send : standard, buffered, immediate
-   CHARACTER(len=1) ::   cn_mpi_send = 'S'    ! type od mpi send/recieve (S=standard, B=bsend, I=isend)
-   LOGICAL          ::   l_isend = .FALSE.   ! isend use indicator (T if cn_mpi_send='I')
-   INTEGER          ::   nn_buffer = 0       ! size of the buffer in case of mpi_bsend 
+   CHARACTER(len=1)         ::   cn_mpi_send = 'S'    ! type od mpi send/recieve (S=standard, B=bsend, I=isend)
+   LOGICAL         , PUBLIC ::   l_isend = .FALSE.   ! isend use indicator (T if cn_mpi_send='I')
+   INTEGER                  ::   nn_buffer = 0       ! size of the buffer in case of mpi_bsend 
       
    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: tampon  ! buffer in case of bsend
@@ -159,6 +160,4 @@
    REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE, SAVE ::   t2ew, t2we   ! 2d for east-west & west-east
    REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE, SAVE ::   t2p1, t2p2   ! 2d for north fold
-   REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE, SAVE ::   tr2ns, tr2sn ! 2d for north-south & south-north + extra outer halo
-   REAL(wp), DIMENSION(:,:,:)    , ALLOCATABLE, SAVE ::   tr2ew, tr2we ! 2d for east-west   & west-east   + extra outer halo
 
    ! Arrays used in mpp_lbc_north_3d()
@@ -207,9 +206,4 @@
          &      t2ew(jpj,jpreci    ,2)   , t2we(jpj,jpreci    ,2)   ,                                            &
          &      t2p1(jpi,jprecj    ,2)   , t2p2(jpi,jprecj    ,2)   ,                                            &
-         !
-         &      tr2ns(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
-         &      tr2sn(1-jpr2di:jpi+jpr2di,jprecj+jpr2dj,2) ,                                                     &
-         &      tr2ew(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2) ,                                                     &
-         &      tr2we(1-jpr2dj:jpj+jpr2dj,jpreci+jpr2di,2) ,                                                     &
          !
          &      ztab(jpiglo,4,jpk) , znorthloc(jpi,4,jpk) , znorthgloio(jpi,4,jpk,jpni) ,                        &
@@ -947,5 +941,5 @@
 
 
-   SUBROUTINE mpp_lnk_2d_e( pt2d, cd_type, psgn )
+   SUBROUTINE mpp_lnk_2d_e( pt2d, cd_type, psgn, jpri, jprj )
       !!----------------------------------------------------------------------
       !!                  ***  routine mpp_lnk_2d_e  ***
@@ -958,6 +952,6 @@
       !!                    nlci   : first dimension of the local subdomain
       !!                    nlcj   : second dimension of the local subdomain
-      !!                    jpr2di : number of rows for extra outer halo
-      !!                    jpr2dj : number of columns for extra outer halo
+      !!                    jpri   : number of rows for extra outer halo
+      !!                    jprj   : number of columns for extra outer halo
       !!                    nbondi : mark for "east-west local boundary"
       !!                    nbondj : mark for "north-south local boundary"
@@ -968,9 +962,11 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
-      CHARACTER(len=1)                                            , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
-      !                                                                                         ! = T , U , V , F , W and I points
-      REAL(wp)                                                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the
-      !!                                                                                        ! north boundary, =  1. otherwise
+      INTEGER                                             , INTENT(in   ) ::   jpri
+      INTEGER                                             , INTENT(in   ) ::   jprj
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,1-jprj:jpj+jprj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)                                    , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      !                                                                                 ! = T , U , V , F , W and I points
+      REAL(wp)                                            , INTENT(in   ) ::   psgn     ! =-1 the sign change across the
+      !!                                                                                ! north boundary, =  1. otherwise
       INTEGER  ::   jl   ! dummy loop indices
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
@@ -978,8 +974,13 @@
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
-      !!----------------------------------------------------------------------
-
-      ipreci = jpreci + jpr2di      ! take into account outer extra 2D overlap area
-      iprecj = jprecj + jpr2dj
+      !!
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,jprecj+jprj,2) :: r2dns
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,jprecj+jprj,2) :: r2dsn
+      REAL(wp), DIMENSION(1-jprj:jpj+jprj,jpreci+jpri,2) :: r2dwe
+      REAL(wp), DIMENSION(1-jprj:jpj+jprj,jpreci+jpri,2) :: r2dew
+      !!----------------------------------------------------------------------
+
+      ipreci = jpreci + jpri      ! take into account outer extra 2D overlap area
+      iprecj = jprecj + jprj
 
 
@@ -989,16 +990,16 @@
       !
       !                                      !* North-South boundaries (always colsed)
-      IF( .NOT. cd_type == 'F' )   pt2d(:,  1-jpr2dj   :  jprecj  ) = 0.e0    ! south except at F-point
-                                   pt2d(:,nlcj-jprecj+1:jpj+jpr2dj) = 0.e0    ! north
+      IF( .NOT. cd_type == 'F' )   pt2d(:,  1-jprj   :  jprecj  ) = 0.e0    ! south except at F-point
+                                   pt2d(:,nlcj-jprecj+1:jpj+jprj) = 0.e0    ! north
                                 
       !                                      ! East-West boundaries
       !                                           !* Cyclic east-west
       IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
-         pt2d(1-jpr2di:     1    ,:) = pt2d(jpim1-jpr2di:  jpim1 ,:)       ! east
-         pt2d(   jpi  :jpi+jpr2di,:) = pt2d(     2      :2+jpr2di,:)       ! west
+         pt2d(1-jpri:     1    ,:) = pt2d(jpim1-jpri:  jpim1 ,:)       ! east
+         pt2d(   jpi  :jpi+jpri,:) = pt2d(     2      :2+jpri,:)       ! west
          !
       ELSE                                        !* closed
-         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpr2di   :jpreci    ,:) = 0.e0    ! south except at F-point
-                                      pt2d(nlci-jpreci+1:jpi+jpr2di,:) = 0.e0    ! north
+         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0.e0    ! south except at F-point
+                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0.e0    ! north
       ENDIF
       !
@@ -1009,5 +1010,5 @@
          !
          SELECT CASE ( jpni )
-         CASE ( 1 )     ;   CALL lbc_nfd        ( pt2d(1:jpi,1:jpj+jpr2dj), cd_type, psgn, pr2dj=jpr2dj )
+         CASE ( 1 )     ;   CALL lbc_nfd        ( pt2d(1:jpi,1:jpj+jprj), cd_type, psgn, pr2dj=jprj )
          CASE DEFAULT   ;   CALL mpp_lbc_north_e( pt2d                    , cd_type, psgn               )
          END SELECT 
@@ -1021,29 +1022,29 @@
       SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
       CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
-         iihom = nlci-nreci-jpr2di
+         iihom = nlci-nreci-jpri
          DO jl = 1, ipreci
-            tr2ew(:,jl,1) = pt2d(jpreci+jl,:)
-            tr2we(:,jl,1) = pt2d(iihom +jl,:)
+            r2dew(:,jl,1) = pt2d(jpreci+jl,:)
+            r2dwe(:,jl,1) = pt2d(iihom +jl,:)
          END DO
       END SELECT
       !
       !                           ! Migrations
-      imigr = ipreci * ( jpj + 2*jpr2dj)
+      imigr = ipreci * ( jpj + 2*jprj)
       !
       SELECT CASE ( nbondi )
       CASE ( -1 )
-         CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea )
+         CALL mppsend( 2, r2dwe(1-jprj,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, r2dew(1-jprj,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea )
-         CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe )
+         CALL mppsend( 1, r2dew(1-jprj,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, r2dwe(1-jprj,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, r2dew(1-jprj,1,2), imigr, noea )
+         CALL mpprecv( 2, r2dwe(1-jprj,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe )
+         CALL mppsend( 1, r2dew(1-jprj,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, r2dwe(1-jprj,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -1055,14 +1056,14 @@
       CASE ( -1 )
          DO jl = 1, ipreci
-            pt2d(iihom+jl,:) = tr2ew(:,jl,2)
+            pt2d(iihom+jl,:) = r2dew(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, ipreci
-            pt2d(jl-jpr2di,:) = tr2we(:,jl,2)
-            pt2d( iihom+jl,:) = tr2ew(:,jl,2)
+            pt2d(jl-jpri,:) = r2dwe(:,jl,2)
+            pt2d( iihom+jl,:) = r2dew(:,jl,2)
          END DO
       CASE ( 1 )
          DO jl = 1, ipreci
-            pt2d(jl-jpr2di,:) = tr2we(:,jl,2)
+            pt2d(jl-jpri,:) = r2dwe(:,jl,2)
          END DO
       END SELECT
@@ -1074,29 +1075,29 @@
       !
       IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
-         ijhom = nlcj-nrecj-jpr2dj
+         ijhom = nlcj-nrecj-jprj
          DO jl = 1, iprecj
-            tr2sn(:,jl,1) = pt2d(:,ijhom +jl)
-            tr2ns(:,jl,1) = pt2d(:,jprecj+jl)
+            r2dsn(:,jl,1) = pt2d(:,ijhom +jl)
+            r2dns(:,jl,1) = pt2d(:,jprecj+jl)
          END DO
       ENDIF
       !
       !                           ! Migrations
-      imigr = iprecj * ( jpi + 2*jpr2di )
+      imigr = iprecj * ( jpi + 2*jpri )
       !
       SELECT CASE ( nbondj )
       CASE ( -1 )
-         CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono )
+         CALL mppsend( 4, r2dsn(1-jpri,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, r2dns(1-jpri,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono )
-         CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso )
+         CALL mppsend( 3, r2dns(1-jpri,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, r2dsn(1-jpri,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, r2dns(1-jpri,1,2), imigr, nono )
+         CALL mpprecv( 4, r2dsn(1-jpri,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso )
+         CALL mppsend( 3, r2dns(1-jpri,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, r2dsn(1-jpri,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -1108,14 +1109,14 @@
       CASE ( -1 )
          DO jl = 1, iprecj
-            pt2d(:,ijhom+jl) = tr2ns(:,jl,2)
+            pt2d(:,ijhom+jl) = r2dns(:,jl,2)
          END DO
       CASE ( 0 )
          DO jl = 1, iprecj
-            pt2d(:,jl-jpr2dj) = tr2sn(:,jl,2)
-            pt2d(:,ijhom+jl ) = tr2ns(:,jl,2)
+            pt2d(:,jl-jprj) = r2dsn(:,jl,2)
+            pt2d(:,ijhom+jl ) = r2dns(:,jl,2)
          END DO
       CASE ( 1 ) 
          DO jl = 1, iprecj
-            pt2d(:,jl-jpr2dj) = tr2sn(:,jl,2)
+            pt2d(:,jl-jprj) = r2dsn(:,jl,2)
          END DO
       END SELECT
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 3609)
@@ -39,4 +39,5 @@
    LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS      
    LOGICAL , PUBLIC ::   ln_apr_dyn  = .FALSE.   !: Atmospheric pressure forcing used on dynamics (ocean & ice)
+   LOGICAL , PUBLIC ::   ln_icebergs = .FALSE.   !: Icebergs
    INTEGER , PUBLIC ::   nn_ice      = 0         !: flag on ice in the surface boundary condition (=0/1/2/3)
    INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget: 
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 3609)
@@ -43,4 +43,5 @@
    USE bdy_par          ! for lk_bdy
    USE bdyice_lim2      ! unstructured open boundary data  (bdy_ice_lim_2 routine)
+   USE icbstp           ! Icebergs!
 
    USE prtctl           ! Print control                    (prt_ctl routine)
@@ -290,4 +291,6 @@
       CASE(  4 )   ;       CALL sbc_ice_cice ( kt, nsbc )            ! CICE ice model
       END SELECT                                              
+
+      IF( ln_icebergs )    CALL icb_stp( kt )                   ! compute icebergs
 
       IF( ln_rnf       )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90	(revision 3609)
@@ -284,10 +284,10 @@
 
       IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN
-         CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1. )   ! lateral boundary conditions
-         CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1. )   ! lateral boundary conditions
-         CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1. )   ! lateral boundary conditions
-         CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1. )   ! lateral boundary conditions
-         CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions
-         CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions         
+         CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
+         CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
+         CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
+         CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
+         CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
+         CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions         
          IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
       END IF
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90	(revision 3609)
@@ -81,5 +81,5 @@
          !                                                    ! ==============
 
-         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
+         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
         
          ! Residus
@@ -104,5 +104,5 @@
          icount = icount + 1 
  
-         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
+         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
 
          ! Guess red update
@@ -167,5 +167,5 @@
       !  Output in gcx
       !  -------------
-      CALL lbc_lnk_e( gcx, c_solver_pt, 1. )    ! boundary conditions
+      CALL lbc_lnk_e( gcx, c_solver_pt, 1._wp, jpr2di, jpr2dj )    ! boundary conditions
       !
       CALL wrk_dealloc( jpi, jpj, ztab )
Index: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3608)
+++ branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3609)
@@ -61,4 +61,6 @@
    USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
    USE step            ! NEMO time-stepping                 (stp     routine)
+   USE icbini          ! handle bergs, initialisation
+   USE icbstp          ! handle bergs, calving, themodynamics and transport
 #if defined key_oasis3
    USE cpl_oasis3      ! OASIS3 coupling
@@ -161,5 +163,7 @@
 #endif
 
-      IF( lk_diaobs ) CALL dia_obs_wri
+      IF( lk_diaobs   )   CALL dia_obs_wri
+      !
+      IF( ln_icebergs )   CALL icb_end( nitend )
        
       !                            !------------------------!
@@ -360,4 +364,5 @@
       !                                     ! Misc. options
       IF( nn_cla == 1   )   CALL cla_init       ! Cross Land Advection
+                            CALL icb_init( rdt, nit000)   ! initialise icebergs instance
       
 #if defined key_top
