Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 3339)
@@ -0,0 +1,162 @@
+MODULE icb_oce
+
+   !!======================================================================
+   !!                       ***  MODULE  icb_oce  ***
+   !! Ocean physics:  declare 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)       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 fldread
+
+   IMPLICIT NONE
+   PUBLIC
+
+   ! Global constants
+   INTEGER, PUBLIC, PARAMETER :: nclasses = 10 ! Number of ice bergs classes
+   INTEGER, PUBLIC, PARAMETER :: nkounts  =  3 ! Number of integers combined for unique naming
+
+   !! icebergs_gridded type allows for interaction between the set of icebergs and the model grid
+   TYPE, PUBLIC :: icebergs_gridded
+      REAL(wp), DIMENSION(:,:)  , POINTER :: calving => NULL()     ! Calving mass rate [frozen runoff] (kg/s) (into stored ice)
+      REAL(wp), DIMENSION(:,:)  , POINTER :: calving_hflx=>NULL()  ! Calving heat flux [heat content of calving] (W/m2)
+      REAL(wp), DIMENSION(:,:)  , POINTER :: floating_melt=>NULL() ! Net melting rate to icebergs + bits (kg/s/m^2)
+      INTEGER , DIMENSION(:,:)  , POINTER :: maxclass=>NULL()      ! maximum class number at calving source point
+      REAL(wp), DIMENSION(:,:)  , POINTER :: tmp=>NULL()           ! Temporary work space
+      REAL(wp), DIMENSION(:,:,:), POINTER :: stored_ice=>NULL()    ! Accumulated ice mass flux at calving locations (kg)
+      REAL(wp), DIMENSION(:,:)  , POINTER :: stored_heat=>NULL()   ! Heat content of stored ice (J)
+   END TYPE icebergs_gridded
+
+   !! point type is just a convenience storage type for properties of an individual iceberg
+   TYPE, PUBLIC ::   point
+      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
+      REAL(wp) :: mass, thickness, width, length, uvel, vvel
+      REAL(wp) :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi
+      REAL(wp) :: mass_of_bits, heat_density
+   END TYPE point
+
+   !! iceberg type is really just a linked list
+   TYPE, PUBLIC :: iceberg
+      ! pointers to previous and next unique icebergs in linked list
+      TYPE(iceberg), POINTER      :: prev=>NULL(), next=>NULL()
+      ! variables which do not change for this iceberg
+      INTEGER, DIMENSION(nkounts) :: number
+      REAL(wp)                    :: mass_scaling
+      ! variables which change with time are held in a separate type
+      TYPE(point), POINTER        :: current_point => NULL()
+   END TYPE iceberg
+
+   !! master instances of gridded and linked list iceberg types
+   TYPE(icebergs_gridded), POINTER ::   berg_grid
+   TYPE(iceberg)         , POINTER ::   first_berg   => NULL()
+
+   !! parameters controlling iceberg characteristics and modelling
+   REAL(wp)                        ::   berg_dt                    ! Time-step between iceberg CALLs (should make adaptive?)
+   INTEGER                         ::   current_year
+   REAL(wp)                        ::   current_yearday            ! 1.00-365.99
+   REAL(wp), DIMENSION(:), POINTER ::   initial_width, initial_length
+   LOGICAL                         ::   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
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: vo_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ff_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ua_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: va_e
+#if defined key_lim2 || defined key_lim3
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ui_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: vi_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_e
+#endif
+
+   !!gm almost all those PARAM ARE defined in NEMO
+   REAL(wp), PUBLIC, PARAMETER :: rho_ice      = 916.7_wp   ! Density of fresh ice @ 0oC (kg/m^3)
+   REAL(wp), PUBLIC, PARAMETER :: rho_water    = 999.8_wp   ! Density of fresh water @ 0oC (kg/m^3)
+   REAL(wp), PUBLIC, PARAMETER :: rho_air      = 1.1_wp     ! Density of air @ 0oC (kg/m^3) ???
+   REAL(wp), PUBLIC, PARAMETER :: rho_seawater = 1025._wp   ! Approx. density of surface sea water @ 0oC (kg/m^3)
+   !!gm end
+   REAL(wp), PUBLIC, PARAMETER :: Cd_av = 1.3_wp      ! (Vertical) Drag coefficient between bergs and atmos (?)
+   REAL(wp), PUBLIC, PARAMETER :: Cd_ah = 0.0055_wp   ! (Horizontal) Drag coefficient between bergs and atmos (?)
+   REAL(wp), PUBLIC, PARAMETER :: Cd_wv = 0.9_wp      ! (Vertical) Drag coefficient between bergs and ocean (?)
+   REAL(wp), PUBLIC, PARAMETER :: Cd_wh = 0.0012_wp   ! (Horizontal) Drag coefficient between bergs and ocean (?)
+   REAL(wp), PUBLIC, PARAMETER :: Cd_iv = 0.9_wp      ! (Vertical) Drag coefficient between bergs and sea-ice (?)
+   !TOM> no horizontal drag for sea ice! real, PARAMETER :: Cd_ih=0.0012 ! (Horizontal) Drag coefficient 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 number
+   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.8e7, 4.1e8, 3.3e9, 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), ALLOCATABLE, PUBLIC, SAVE, DIMENSION(:,:)     ::   p_calving, p_calving_hflx
+   INTEGER              , PUBLIC, SAVE                     ::   numicb 
+   INTEGER              , PUBLIC, SAVE, DIMENSION(nkounts) ::   kount_bergs
+   INTEGER              , PUBLIC, SAVE                     ::   icbdi, icbei, icbdj, icbej
+   REAL(wp)             , PUBLIC, SAVE                     ::   icb_left, icb_right
+   INTEGER              , PUBLIC, SAVE                     ::   icbpack
+   INTEGER              , PUBLIC, SAVE                     ::   ktberg, knberg
+   INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:)       ::   icbfldpts
+   INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:)       ::   icbflddest
+   INTEGER, ALLOCATABLE , PUBLIC, SAVE, DIMENSION(:)       ::   icbfldproc
+
+   TYPE(FLD), ALLOCATABLE, PUBLIC     , DIMENSION(:)       ::   sf_icb       ! structure: file information, fields read
+
+END MODULE icb_oce
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 3339)
@@ -0,0 +1,182 @@
+MODULE icbclv
+
+   !!======================================================================
+   !!                       ***  MODULE  icbclv  ***
+   !! Ocean physics:  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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   accumulate_calving      :
+   !!   icb_gen       : generate test icebergs
+   !!   icb_nam       : read iceberg namelist
+   !!----------------------------------------------------------------------
+   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        ! define iceberg arrays
+   USE icbdia         ! iceberg utility routines
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   accumulate_calving  ! routine called in xxx.F90 module
+   PUBLIC   calve_icebergs      ! routine called in xxx.F90 module
+
+CONTAINS
+
+   SUBROUTINE accumulate_calving( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE accumulate_calving  ***
+      !!
+      !! ** Purpose :   ?
+      !!
+      !! ** input   : - ?
+      !!----------------------------------------------------------------------
+      !
+      INTEGER                         :: kt
+      REAL(wp)                        :: calving_used, rdist, ufact
+      INTEGER                         :: jn, ji, jj, nmax
+      LOGICAL, SAVE                   :: 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
+      ufact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 1000._wp
+      berg_grid%calving(:,:) = p_calving(:,:) * tmask_i(:,:) * ufact
+
+      ! Heat in units of W/m2, and mask (just in case)
+      berg_grid%calving_hflx(:,:) = p_calving_hflx(:,:) * tmask_i(:,:)
+
+      IF( first_call .AND. .NOT.restarted_bergs) THEN      ! This is a hack to simplify initialization
+         first_call = .FALSE.
+         !do jn=1, nclasses
+         !  where (berg_grid%calving==0.) berg_grid%stored_ice(:,:,jn)=0.
+         !enddo
+         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
+            nmax = berg_grid%maxclass(ji,jj)
+            rdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:nmax) )
+            DO jn = 1, nmax
+               berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) + &
+                                          berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * rdist
+            END DO
+         END DO
+      END DO
+
+      ! before changing the calving, save the amount we're about to use and do budget
+      calving_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 incoming_budget( kt,  calving_used, berg_grid%tmp )
+      !
+   END SUBROUTINE accumulate_calving
+
+! ##############################################################################
+
+   SUBROUTINE calve_icebergs()
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE calve_icebergs  ***
+      !!
+      !! ** Purpose :   This seems to be the routine that takes a stored ice field and calves to the ocean,
+      !!                so I assume that 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 (largest first?)
+      !!                Note that only the non-overlapping part of the processor where icebergs are allowed
+      !!                is considered
+      !!----------------------------------------------------------------------
+      !
+      INTEGER                         ::   icnt,icntmax
+      TYPE(iceberg)                   ::   newberg
+      TYPE(point)                     ::   newpt
+      LOGICAL                         ::   lret
+      INTEGER                         ::   ji, jj, jn   ! dummy loop indices
+      REAL(wp)                        ::   xi, yj, ddt, calved_to_berg, heat_to_berg
+      !!----------------------------------------------------------------------
+
+      icntmax          = 0
+
+      DO jn = 1, nclasses
+         DO jj = icbdj, icbej
+            DO ji = icbdi, icbei
+               !
+               ddt  = 0._wp
+               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( nimpp+ji-1, wp )
+                  newpt%yj  = REAL( njmpp+jj-1, 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        = initial_width    (jn)
+                  newpt%length       = initial_length   (jn)
+                  newberg%mass_scaling = rn_mass_scaling     (jn)
+                  newpt%mass_of_bits = 0._wp                          ! no bergy
+                  !
+                  newpt%year   = current_year
+                  newpt%day    = current_yearday + ddt / rday
+                  newpt%heat_density = berg_grid%stored_heat(ji,jj) / berg_grid%stored_ice(ji,jj,jn)   ! This is in J/kg
+                  !
+                  CALL increment_kounter()
+                  newberg%number(:) = kount_bergs(:)
+                  !
+                  CALL add_new_berg_to_list( newberg, newpt )
+                  !
+                  calved_to_berg = rn_initial_mass(jn) * rn_mass_scaling(jn)           ! Units of kg
+                  !                                ! Heat content
+                  heat_to_berg           = calved_to_berg * newpt%heat_density             ! Units of J
+                  berg_grid%stored_heat(ji,jj) = berg_grid%stored_heat(ji,jj) - heat_to_berg
+                  !                                ! Stored mass
+                  berg_grid%stored_ice(ji,jj,jn)   = berg_grid%stored_ice(ji,jj,jn) - calved_to_berg
+                  !
+                  ddt  = ddt + berg_dt * 2._wp / 17._wp    ! Minor offset to start day    !!gm why???
+                  icnt = icnt + 1
+                  !
+                  CALL calving_budget(ji, jj, jn,  calved_to_berg, heat_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 )
+      ENDDO
+      CALL lbc_lnk( berg_grid%stored_heat, 'T', 1._wp )
+      !
+      IF( nn_verbose_level > 0 .AND. icntmax > 1 )   WRITE(numicb,*) 'calve_icebergs: icnt=', icnt,' on', narea
+      !
+   END SUBROUTINE  calve_icebergs
+
+END MODULE icbclv
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 3339)
@@ -0,0 +1,568 @@
+MODULE icbdia
+
+   !!======================================================================
+   !!                       ***  MODULE  icbdia  ***
+   !! Ocean physics:  initialise variables for iceberg budgets and diagnostics
+   !!======================================================================
+   !! 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 are now all here with lots
+   !!            -    !                            of silly routines to call to get values in
+   !!            -    !                            from the right points in the code
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !! icb_budget_end  : end        iceberg budgeting
+   !! icb_budget_init : initialise iceberg budgeting
+   !!----------------------------------------------------------------------
+   USE par_oce        ! nemo parameters
+   USE dom_oce        ! ocean domain
+   USE in_out_manager ! nemo IO
+   USE lib_mpp
+   USE iom
+
+   USE icb_oce        ! define iceberg arrays
+   USE icbutl         ! iceberg utility routines
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: berg_melt=>NULL()    ! Melting+erosion rate of icebergs (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: melt_buoy=>NULL()    ! Buoyancy component of melting rate (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: melt_eros=>NULL()    ! Erosion component of melting rate (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: melt_conv=>NULL()    ! Convective component of melting rate (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: bits_src=>NULL()    ! Mass flux from berg erosion into bergy bits (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: bits_melt=>NULL()   ! Melting rate of bergy bits (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: bits_mass=>NULL()   ! Mass distribution of bergy bits (kg/s/m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: virtual_area=>NULL() ! Virtual surface coverage by icebergs (m^2)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PUBLIC  :: berg_mass=>NULL()         ! Mass distribution (kg/m^2)
+   REAL(wp), DIMENSION(:,:,:), POINTER, PUBLIC  :: real_calving=>NULL() ! Calving rate into iceberg class at calving locations (kg/s)
+   REAL(wp), DIMENSION(:,:)  , POINTER, PRIVATE :: tmpc=>NULL()         ! Temporary work space
+   REAL(wp), DIMENSION(:)    , POINTER, PRIVATE :: zsumbuf=>NULL()         ! Temporary work space
+   INTEGER , DIMENSION(:)    , POINTER, PRIVATE :: isumbuf=>NULL()         ! Temporary work space
+
+   REAL(wp)                           , PRIVATE ::  net_berg_melt
+   REAL(wp)                           , PRIVATE ::  net_bits_src
+   REAL(wp)                           , PRIVATE ::  net_bits_melt
+   REAL(wp)                           , PRIVATE ::  bits_mass_start, bits_mass_end
+   REAL(wp)                           , PRIVATE ::  floating_heat_start, floating_heat_end
+   REAL(wp)                           , PRIVATE ::  floating_mass_start, floating_mass_end
+   REAL(wp)                           , PRIVATE ::  icebergs_mass_start, icebergs_mass_end
+   REAL(wp)                           , PRIVATE ::  stored_start, stored_heat_start
+   REAL(wp)                           , PRIVATE ::  stored_end  , stored_heat_end
+   REAL(wp)                           , PRIVATE ::  net_incoming_calving, net_outgoing_calving
+   REAL(wp)                           , PRIVATE ::  net_incoming_calving_heat, net_outgoing_calving_heat
+   REAL(wp)                           , PRIVATE ::  net_incoming_calving_heat_used
+   REAL(wp)                           , PRIVATE ::  net_calving_received, net_calving_returned, net_calving_used
+   REAL(wp)                           , PRIVATE ::  net_heat_to_bergs, net_heat_to_ocean, net_melt
+   REAL(wp)                           , PRIVATE ::  net_calving_to_bergs
+
+   INTEGER                            , PRIVATE ::  nbergs_start, nbergs_end, nbergs_calved
+   INTEGER                            , PRIVATE ::  nbergs_melted
+   INTEGER                            , PRIVATE ::  nspeeding_tickets
+   INTEGER , DIMENSION(nclasses)      , PRIVATE ::  nbergs_calved_by_class
+
+   PUBLIC   icb_budget_end  ! routine called in xxx.F90 module
+   PUBLIC   icb_budget_init ! routine called in xxx.F90 module
+   PUBLIC   icb_budget      ! routine called in xxx.F90 module
+   PUBLIC   icb_budget_step ! routine called in xxx.F90 module
+   PUBLIC   icb_budget_put  ! routine called in xxx.F90 module
+   PUBLIC   melt_budget     ! routine called in xxx.F90 module
+   PUBLIC   size_budget     ! routine called in xxx.F90 module
+   PUBLIC   speed_budget     ! routine called in xxx.F90 module
+   PUBLIC   calving_budget  ! routine called in xxx.F90 module
+   PUBLIC   incoming_budget   ! routine called in xxx.F90 module
+
+CONTAINS
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_budget_end
+
+      IF( .NOT. ln_bergdia ) RETURN
+      DEALLOCATE( berg_melt )
+      DEALLOCATE( melt_buoy )
+      DEALLOCATE( melt_eros )
+      DEALLOCATE( melt_conv )
+      DEALLOCATE( bits_src )
+      DEALLOCATE( bits_melt )
+      DEALLOCATE( bits_mass )
+      DEALLOCATE( virtual_area )
+      DEALLOCATE( berg_mass )
+      DEALLOCATE( real_calving )
+      DEALLOCATE( tmpc )
+      IF( lk_mpp ) THEN
+         DEALLOCATE( zsumbuf )
+         DEALLOCATE( isumbuf )
+      ENDIF
+
+   END SUBROUTINE icb_budget_end
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_budget_init( )
+
+      IF( .NOT. ln_bergdia ) RETURN
+      ALLOCATE( berg_melt    (jpi,jpj)   )           ;   berg_melt    (:,:) = 0._wp
+      ALLOCATE( melt_buoy    (jpi,jpj)   )           ;   melt_buoy    (:,:) = 0._wp
+      ALLOCATE( melt_eros    (jpi,jpj)   )           ;   melt_eros    (:,:) = 0._wp
+      ALLOCATE( melt_conv    (jpi,jpj)   )           ;   melt_conv    (:,:) = 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.
+      ALLOCATE( tmpc(jpi,jpj) )                      ;   tmpc(:,:)=0.
+
+      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
+      icebergs_mass_end              = 0._wp
+      bits_mass_end                 = 0._wp
+      stored_heat_start              = 0._wp
+      floating_heat_start            = 0._wp
+      floating_mass_start            = 0._wp
+      icebergs_mass_start            = 0._wp
+      bits_mass_start               = 0._wp
+      bits_mass_end                 = 0._wp
+      net_calving_used               = 0._wp
+      net_calving_to_bergs           = 0._wp
+      net_heat_to_bergs              = 0._wp
+      net_heat_to_ocean              = 0._wp
+      net_calving_received           = 0._wp
+      net_calving_returned           = 0._wp
+      net_incoming_calving           = 0._wp
+      net_outgoing_calving           = 0._wp
+      net_incoming_calving_heat      = 0._wp
+      net_incoming_calving_heat_used = 0._wp
+      net_outgoing_calving_heat      = 0._wp
+      net_melt                       = 0._wp
+      net_berg_melt                  = 0._wp
+      net_bits_melt                 = 0._wp
+      net_bits_src                  = 0._wp
+
+      floating_mass_start            = sum_mass( first_berg )
+      icebergs_mass_start            = sum_mass( first_berg, justbergs=.true. )
+      bits_mass_start               = sum_mass( first_berg, justbits=.true. )
+      IF( lk_mpp ) THEN
+         ALLOCATE( zsumbuf(23) )          ; zsumbuf(:) = 0._wp
+         ALLOCATE( isumbuf(4+nclasses) )  ; isumbuf(:) = 0
+         zsumbuf(1) = floating_mass_start
+         zsumbuf(2) = icebergs_mass_start
+         zsumbuf(3) = bits_mass_start
+         CALL mpp_sum( zsumbuf(1:3), 3 )
+         floating_mass_start = zsumbuf(1)
+         icebergs_mass_start = zsumbuf(2)
+         bits_mass_start = zsumbuf(3)
+      ENDIF
+
+   END SUBROUTINE icb_budget_init
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_budget( lbudge )
+      ! Arguments
+      LOGICAL                         ::   lbudge
+      ! Local variables
+      INTEGER                         ::   k
+      REAL(wp)                        ::   unused_calving, tmpsum, grdd_berg_mass, grdd_bits_mass
+
+      IF( .NOT. ln_bergdia ) RETURN
+
+      unused_calving            = SUM( berg_grid%calving(:,:) )
+      tmpsum                    = SUM( berg_grid%floating_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      net_melt                  = net_melt + tmpsum * berg_dt
+      net_outgoing_calving      = net_outgoing_calving + ( unused_calving + tmpsum ) * berg_dt
+      tmpsum                    = SUM( berg_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      net_berg_melt             = net_berg_melt + tmpsum * berg_dt
+      tmpsum                    = SUM( bits_src(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      net_bits_src             = net_bits_src + tmpsum * berg_dt
+      tmpsum                    = SUM( bits_melt(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      net_bits_melt            = net_bits_melt + tmpsum * berg_dt
+      tmpsum                    = SUM( p_calving(:,:) * tmask_i(:,:) )
+      net_calving_returned      = net_calving_returned+tmpsum * berg_dt
+      tmpsum                    = SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) )
+      net_outgoing_calving_heat = net_outgoing_calving_heat + tmpsum * berg_dt   ! Units of J
+
+      IF( lbudge ) THEN
+         stored_end             = SUM( berg_grid%stored_ice(:,:,:) )
+         stored_heat_end        = SUM( berg_grid%stored_heat(:,:) )
+         floating_mass_end      = sum_mass( first_berg )
+         icebergs_mass_end      = sum_mass( first_berg,justbergs=.true. )
+         bits_mass_end         = sum_mass( first_berg,justbits=.true. )
+         floating_heat_end      = sum_heat( first_berg )
+
+         nbergs_end             = count_bergs()
+         grdd_berg_mass         = SUM( berg_mass (:,:)*e1e2t(:,:)*tmask_i(:,:) )
+         grdd_bits_mass        = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) )
+
+         IF( lk_mpp ) THEN
+            zsumbuf( 1) = stored_end
+            zsumbuf( 2) = stored_heat_end
+            zsumbuf( 3) = floating_mass_end
+            zsumbuf( 4) = icebergs_mass_end
+            zsumbuf( 5) = bits_mass_end
+            zsumbuf( 6) = floating_heat_end
+            zsumbuf( 7) = net_calving_returned
+            zsumbuf( 8) = net_outgoing_calving
+            zsumbuf( 9) = net_calving_received
+            zsumbuf(10) = net_incoming_calving
+            zsumbuf(11) = net_incoming_calving_heat
+            zsumbuf(12) = net_incoming_calving_heat_used
+            zsumbuf(13) = net_outgoing_calving_heat
+            zsumbuf(14) = net_calving_used
+            zsumbuf(15) = net_calving_to_bergs
+            zsumbuf(16) = net_heat_to_bergs
+            zsumbuf(17) = net_heat_to_ocean
+            zsumbuf(18) = net_melt
+            zsumbuf(19) = net_berg_melt
+            zsumbuf(20) = net_bits_src
+            zsumbuf(21) = net_bits_melt
+            zsumbuf(22) = grdd_berg_mass
+            zsumbuf(23) = grdd_bits_mass
+
+            CALL mpp_sum( zsumbuf(1:23), 23)
+
+            stored_end                     = zsumbuf( 1)
+            stored_heat_end                = zsumbuf( 2)
+            floating_mass_end              = zsumbuf( 3)
+            icebergs_mass_end              = zsumbuf( 4)
+            bits_mass_end                 = zsumbuf( 5)
+            floating_heat_end              = zsumbuf( 6)
+            net_calving_returned           = zsumbuf( 7)
+            net_outgoing_calving           = zsumbuf( 8)
+            net_calving_received           = zsumbuf( 9)
+            net_incoming_calving           = zsumbuf(10)
+            net_incoming_calving_heat      = zsumbuf(11)
+            net_incoming_calving_heat_used = zsumbuf(12)
+            net_outgoing_calving_heat      = zsumbuf(13)
+            net_calving_used               = zsumbuf(14)
+            net_calving_to_bergs           = zsumbuf(15)
+            net_heat_to_bergs              = zsumbuf(16)
+            net_heat_to_ocean              = zsumbuf(17)
+            net_melt                       = zsumbuf(18)
+            net_berg_melt                  = zsumbuf(19)
+            net_bits_src                  = zsumbuf(20)
+            net_bits_melt                 = zsumbuf(21)
+            grdd_berg_mass                 = zsumbuf(22)
+            grdd_bits_mass                = zsumbuf(23)
+
+            isumbuf(1) = nbergs_end
+            isumbuf(2) = nbergs_calved
+            isumbuf(3) = nbergs_melted
+            isumbuf(4) = nspeeding_tickets
+            DO k = 1,nclasses
+               isumbuf(4+k) = nbergs_calved_by_class(k)
+            ENDDO
+
+            CALL mpp_sum( isumbuf(1:nclasses+4), nclasses+4 )
+
+            nbergs_end        = isumbuf(1)
+            nbergs_calved     = isumbuf(2)
+            nbergs_melted     = isumbuf(3)
+            nspeeding_tickets = isumbuf(4)
+            DO k = 1,nclasses
+               nbergs_calved_by_class(k)= isumbuf(4+k)
+            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','',icebergs_mass_start,'',icebergs_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',net_calving_used, &
+                                               'bergs',net_calving_to_bergs, &
+                                               'stored mass',stored_start,stored_end)
+         CALL report_budget('floating mass','kg','calving used',net_calving_to_bergs, &
+                                                 'bergs',net_melt, &
+                                                 'stored mass',floating_mass_start,floating_mass_end)
+         CALL report_budget('berg mass','kg','calving',net_calving_to_bergs, &
+                                             'melt+eros',net_berg_melt, &
+                                             'berg mass',icebergs_mass_start,icebergs_mass_end)
+         CALL report_budget('bits mass','kg','eros used',net_bits_src, &
+                                             'bergs',net_bits_melt, &
+                                             'stored mass',bits_mass_start,bits_mass_end)
+         CALL report_budget('net mass','kg','recvd',net_calving_received, &
+                                            'rtrnd',net_calving_returned, &
+                                            'net mass',stored_start+floating_mass_start, &
+                                                       stored_end+floating_mass_end)
+         CALL report_consistant('iceberg mass','kg','gridded',grdd_berg_mass,'bergs',icebergs_mass_end)
+         CALL report_consistant('bits mass','kg','gridded',grdd_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',net_incoming_calving_heat, &
+                                           'net heat',net_outgoing_calving_heat, &
+                                           'net heat',stored_heat_start+floating_heat_start, &
+                                                      stored_heat_end+floating_heat_end)
+         CALL report_budget('stored heat','J','calving used',net_incoming_calving_heat_used, &
+                                              'bergs',net_heat_to_bergs, &
+                                              'net heat',stored_heat_start,stored_heat_end)
+         CALL report_budget('flting heat','J','calved',net_heat_to_bergs, &
+                                              'melt',net_heat_to_ocean, &
+                                              'net heat',floating_heat_start,floating_heat_end)
+         IF (nn_verbose_level >= 1) THEN
+            CALL report_consistant('top interface','kg','from SIS',net_incoming_calving, &
+                                   'received',net_calving_received)
+            CALL report_consistant('bot interface','kg','sent',net_outgoing_calving, &
+                                   'returned',net_calving_returned)
+         ENDIF
+         WRITE(numicb,'("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(k),k=1,nclasses)
+         IF ( nspeeding_tickets .GT. 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
+         icebergs_mass_start            = icebergs_mass_end
+         bits_mass_start               = bits_mass_end
+         net_calving_used               = 0._wp
+         net_calving_to_bergs           = 0._wp
+         net_heat_to_bergs              = 0._wp
+         net_heat_to_ocean              = 0._wp
+         net_calving_received           = 0._wp
+         net_calving_returned           = 0._wp
+         net_incoming_calving           = 0._wp
+         net_outgoing_calving           = 0._wp
+         net_incoming_calving_heat      = 0._wp
+         net_incoming_calving_heat_used = 0._wp
+         net_outgoing_calving_heat      = 0._wp
+         net_melt                       = 0._wp
+         net_berg_melt                  = 0._wp
+         net_bits_melt                 = 0._wp
+         net_bits_src                  = 0._wp
+      ENDIF
+
+   END SUBROUTINE icb_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_budget_step
+      !! things to reset at the beginning of each timestep
+      !! this probably screws up fields going to diawri, so needs to be looked at - sga
+
+      IF( .NOT. ln_bergdia ) RETURN
+      berg_melt    (:,:)   = 0._wp
+      melt_buoy    (:,:)   = 0._wp
+      melt_eros    (:,:)   = 0._wp
+      melt_conv    (:,:)   = 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_budget_step
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_budget_put
+
+      IF( .NOT. ln_bergdia ) RETURN
+      CALL iom_put( "berg_melt"         , berg_melt    (:,:)   )  ! 'Melt rate of icebergs'                    , 'kg/m2/s'
+      CALL iom_put( "berg_melt_buoy"    , melt_buoy    (:,:)   )  ! 'Buoyancy component of iceberg melt rate'  , 'kg/m2/s'
+      CALL iom_put( "berg_melt_eros"    , melt_eros    (:,:)   )  ! 'Erosion component of iceberg melt rate'   , 'kg/m2/s'
+      CALL iom_put( "berg_melt_conv"    , melt_conv    (:,:)   )  ! '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_budget_put
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE calving_budget( ji, jj, jn, calved, heated )
+      INTEGER   :: ji, jj, jn
+      REAL(wp)  :: calved
+      REAL(wp)  :: heated
+
+      IF( .NOT. ln_bergdia ) RETURN
+      real_calving(ji,jj,jn)     = real_calving(ji,jj,jn) + calved / berg_dt
+      nbergs_calved              = nbergs_calved              + 1
+      nbergs_calved_by_class(jn) = nbergs_calved_by_class(jn) + 1
+      net_calving_to_bergs       = net_calving_to_bergs + calved
+      net_heat_to_bergs          = net_heat_to_bergs    + heated
+
+   END SUBROUTINE calving_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE incoming_budget( kt,  calving_used, heat_used )
+      INTEGER ,                 INTENT(in)  :: kt
+      REAL(wp),                 INTENT(in)  :: calving_used
+      REAL(wp), DIMENSION(:,:), INTENT(in)  :: heat_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)')   'accumulate_calving: 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)')    'accumulate_calving: initial stored heat=',stored_heat_start,' J'
+      ENDIF
+
+      net_calving_received = net_calving_received + SUM( berg_grid%calving(:,:) ) * berg_dt
+      net_incoming_calving = net_calving_received
+      net_incoming_calving_heat = net_incoming_calving_heat +  &
+                                  SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) ) * berg_dt   ! Units of J
+      net_calving_used = net_calving_used + calving_used * berg_dt
+      net_incoming_calving_heat_used = net_incoming_calving_heat_used + SUM( heat_used(:,:) )
+
+   END SUBROUTINE incoming_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE size_budget(ji, jj, Wn, Ln, Abits, mass_scale, Mnew, nMbits, z1_e1e2)
+      INTEGER             :: ji, jj
+      REAL(wp)            :: Wn, Ln, Abits, mass_scale, Mnew, nMbits, z1_e1e2
+
+      IF( .NOT. ln_bergdia ) RETURN
+      virtual_area(ji,jj) = virtual_area(ji,jj)+(Wn*Ln+Abits)*mass_scale        ! m^2
+      berg_mass(ji,jj)    = berg_mass(ji,jj) + Mnew * z1_e1e2                   ! kg/m2
+      bits_mass(ji,jj)   = bits_mass(ji,jj) + nMbits * z1_e1e2                ! kg/m2
+
+   END SUBROUTINE size_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE speed_budget()
+
+      IF( .NOT. ln_bergdia ) RETURN
+      nspeeding_tickets = nspeeding_tickets + 1
+
+   END SUBROUTINE speed_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE melt_budget(ji, jj, mnew, heat, mass_scale, dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2)
+
+      INTEGER               ::  ji, jj
+      REAL(wp), INTENT(in)  ::  mnew, heat, mass_scale
+      REAL(wp), INTENT(in)  ::  dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2
+
+      IF( .NOT. ln_bergdia ) RETURN
+
+      berg_melt    (ji,jj) = berg_melt    (ji,jj) + dM      * z1_dt_e1e2   ! kg/m2/s
+      bits_src    (ji,jj) = bits_src    (ji,jj) + dMbitsE * z1_dt_e1e2   ! mass flux into bergy bitskg/m2/s
+      bits_melt   (ji,jj) = bits_melt   (ji,jj) + dMbitsM * z1_dt_e1e2   ! melt rate of bergy bits kg/m2/s
+      melt_buoy    (ji,jj) = melt_buoy    (ji,jj) + dMb     * z1_dt_e1e2   ! kg/m2/s
+      melt_eros    (ji,jj) = melt_eros    (ji,jj) + dMe     * z1_dt_e1e2   ! erosion rate kg/m2/s
+      melt_conv    (ji,jj) = melt_conv    (ji,jj) + dMv     * z1_dt_e1e2   ! kg/m2/s
+      net_heat_to_ocean = net_heat_to_ocean + heat * mass_scale * berg_dt         ! J
+      IF( mnew <= 0._wp ) nbergs_melted = nbergs_melted + 1                        ! Delete the berg if completely melted
+
+
+   END SUBROUTINE melt_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE report_state(budgetstr,budgetunits,startstr,startval,endstr,endval,delstr,nbergs)
+      ! Arguments
+      CHARACTER*(*), INTENT(in)           :: budgetstr, budgetunits, startstr, endstr, delstr
+      REAL(wp),      INTENT(in)           :: startval, endval
+      INTEGER,       INTENT(in), OPTIONAL :: nbergs
+
+      if (present(nbergs)) then
+        WRITE(numicb,100) budgetstr//' state:', &
+                            startstr//' start',startval,budgetunits, &
+                            endstr//' end',endval,budgetunits, &
+                            'Delta '//delstr,endval-startval,budgetunits, &
+                            '# of bergs',nbergs
+      else
+        WRITE(numicb,100) budgetstr//' state:', &
+                            startstr//' start',startval,budgetunits, &
+                            endstr//' end',endval,budgetunits, &
+                            delstr//'Delta',endval-startval,budgetunits
+      endif
+      100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8)
+   END SUBROUTINE report_state
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE report_consistant(budgetstr,budgetunits,startstr,startval,endstr,endval)
+      ! Arguments
+      CHARACTER*(*), INTENT(in) :: budgetstr, budgetunits, startstr, endstr
+      REAL(wp),      INTENT(in) :: startval, endval
+
+      WRITE(numicb,200) budgetstr//' check:', &
+                        startstr,startval,budgetunits, &
+                        endstr,endval,budgetunits, &
+                        'error',(endval-startval)/((endval+startval)+1e-30),'nd'
+      200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,","))
+   END SUBROUTINE report_consistant
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE report_budget(budgetstr,budgetunits,instr,inval,outstr,outval,delstr,startval,endval)
+      ! Arguments
+      CHARACTER*(*), INTENT(in) :: budgetstr, budgetunits, instr, outstr, delstr
+      REAL(wp),      INTENT(in) :: inval, outval, startval, endval
+
+      WRITE(numicb,200) budgetstr//' budget:', &
+                        instr//' in',inval,budgetunits, &
+                        outstr//' out',outval,budgetunits, &
+                        'Delta '//delstr,inval-outval,budgetunits, &
+                        'error',((endval-startval)-(inval-outval))/max(1.e-30,max(abs(endval-startval),abs(inval-outval))),'nd'
+      200 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a8,"=",es10.3,x,a2)
+   END SUBROUTINE report_budget
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE report_istate(budgetstr,startstr,startval,endstr,endval,delstr)
+      ! Arguments
+      CHARACTER*(*), INTENT(in) :: budgetstr, startstr, endstr, delstr
+      INTEGER,       INTENT(in) :: startval, endval
+
+      WRITE(numicb,100) budgetstr//' state:', &
+                        startstr//' start',startval, &
+                        endstr//' end',endval, &
+                        delstr//'Delta',endval-startval
+      100 FORMAT(a19,3(a18,"=",i14,x,:,","))
+   END SUBROUTINE report_istate
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE report_ibudget(budgetstr,instr,inval,outstr,outval,delstr,startval,endval)
+      ! Arguments
+      CHARACTER*(*), INTENT(in) :: budgetstr, instr, outstr, delstr
+      INTEGER,       INTENT(in) :: inval, outval, startval, endval
+
+      WRITE(numicb,200) budgetstr//' budget:', &
+                          instr//' in',inval, &
+                          outstr//' out',outval, &
+                          'Delta '//delstr,inval-outval, &
+                          'error',((endval-startval)-(inval-outval))
+      200 FORMAT(a19,10(a18,"=",i14,x,:,","))
+   END SUBROUTINE report_ibudget
+   !!-------------------------------------------------------------------------
+
+END MODULE icbdia
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 3339)
@@ -0,0 +1,378 @@
+MODULE icbdyn
+
+   !!======================================================================
+   !!                       ***  MODULE  icbdyn  ***
+   !! Ocean physics:  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)
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_init      : initialise
+   !!   icb_gen       : generate test icebergs
+   !!   icb_nam       : read iceberg namelist
+   !!----------------------------------------------------------------------
+   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   evolve_icebergs  ! routine called in xxx.F90 module
+
+CONTAINS
+
+   SUBROUTINE evolve_icebergs()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE evolve_icebergs  ***
+      !!
+      !! ** Purpose :   iceberg evolution.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      !
+      REAL(wp)                        ::   uvel1 , vvel1 , u1, v1, ax1, ay1, xi1 , yj1
+      REAL(wp)                        ::   uvel2 , vvel2 , u2, v2, ax2, ay2, xi2 , yj2
+      REAL(wp)                        ::   uvel3 , vvel3 , u3, v3, ax3, ay3, xi3 , yj3
+      REAL(wp)                        ::   uvel4 , vvel4 , u4, v4, ax4, ay4, xi4 , yj4
+      REAL(wp)                        ::   uvel_n, vvel_n                  , xi_n, yj_n
+      REAL(wp)                        ::   zdt, zdt_2, zdt_6, e1, e2
+      LOGICAL                         ::   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
+
+         bounced = .false.
+
+
+         ! STEP 1 !
+         ! ====== !
+         xi1 = pt%xi   ;   uvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
+         yj1 = pt%yj   ;   vvel1 = pt%vvel
+
+
+         !                                         !**   A1 = A(X1,V1)
+         CALL accel(        xi1, e1, uvel1, uvel1, ax1,     &
+            &        berg , yj1, e2, vvel1, vvel1, ay1, zdt_2 )
+         !
+         u1 = uvel1 / e1                           !**   V1 in d(i,j)/dt
+         v1 = vvel1 / e2
+
+         ! STEP 2 !
+         ! ====== !
+         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
+         ! position using di/dt & djdt   !   V2  in m/s
+         xi2 = xi1 + zdt_2 * u1          ;   uvel2 = uvel1 + zdt_2 * ax1
+         yj2 = yj1 + zdt_2 * v1          ;   vvel2 = vvel1 + zdt_2 * ay1
+         !
+         CALL adjust_to_ground( xi2, xi1, u1, yj2, yj1, v1, bounced )
+
+         !                                         !**   A2 = A(X2,V2)
+         CALL accel(        xi2, e1, uvel2, uvel1, ax2,    &
+            &        berg , yj2, e2, vvel2, vvel1, ay2, zdt_2 )
+         !
+         u2 = uvel2 / e1                           !**   V2 in d(i,j)/dt
+         v2 = vvel2 / e2
+         !
+         ! STEP 3 !
+         ! ====== !
+         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
+         xi3  = xi1  + zdt_2 * u2   ;   uvel3 = uvel1 + zdt_2 * ax2
+         yj3  = yj1  + zdt_2 * v2   ;   vvel3 = vvel1 + zdt_2 * ay2
+         !
+         CALL adjust_to_ground( xi3, xi1, u3, yj3, yj1, v3, bounced )
+
+         !                                         !**   A3 = A(X3,V3)
+         CALL accel(        xi3, e1, uvel3, uvel1, ax3,    &
+            &        berg , yj3, e2, vvel3, vvel1, ay3, zdt )
+         !
+         u3 = uvel3 / e1                           !**   V3 in d(i,j)/dt
+         v3 = vvel3 / e2
+
+         ! STEP 4 !
+         ! ====== !
+         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
+         xi4 = xi1 + zdt * u3   ;   uvel4 = uvel1 + zdt * ax3
+         yj4 = yj1 + zdt * v3   ;   vvel4 = vvel1 + zdt * ay3
+
+         CALL adjust_to_ground( xi4, xi1, u4, yj4, yj1, v4, bounced )
+
+         !                                         !**   A4 = A(X4,V4)
+         CALL accel(        xi4, e1, uvel4, uvel1, ax4,    &
+            &        berg , yj4, e2, vvel4, vvel1, ay4, zdt )
+
+         u4 = uvel4 / e1                           !**   V4 in d(i,j)/dt
+         v4 = vvel4 / e2
+
+         ! FINAL STEP !
+         ! ========== !
+         !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
+         !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
+         xi_n   = pt%xi   + zdt_6 * (  u1  + 2.*(u2  + u3 ) + u4  )
+         yj_n   = pt%yj   + zdt_6 * (  v1  + 2.*(v2  + v3 ) + v4  )
+         uvel_n = pt%uvel + zdt_6 * (  ax1 + 2.*(ax2 + ax3) + ax4 )
+         vvel_n = pt%vvel + zdt_6 * (  ay1 + 2.*(ay2 + ay3) + ay4 )
+
+         CALL adjust_to_ground( xi_n, xi1, uvel_n, yj_n, yj1, vvel_n, bounced )
+
+         pt%uvel = uvel_n                        !** save in berg structure
+         pt%vvel = vvel_n
+         pt%xi   = xi_n
+         pt%yj   = yj_n
+
+         ! sga - update actual position
+         pt%lon  = bilin_x(glamt, pt%xi, pt%yj )
+         pt%lat  = bilin(gphit, pt%xi, pt%yj, 'T', 0, 0 )
+
+         berg => berg%next                         ! switch to the next berg
+         !
+      END DO                                  !==  end loop over all bergs  ==!
+      !
+   END SUBROUTINE evolve_icebergs
+
+
+   SUBROUTINE adjust_to_ground( pi, pi0, pu, pj, pj0, pv, bounced )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE adjust_to_ground  ***
+      !!
+      !! ** Purpose :   iceberg grounding.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      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)         ::   bounced    ! bounced indicator
+      !
+      REAL(wp), PARAMETER ::   posn_eps = 0.05_wp     ! bouncing distance in (i,j) referential
+      !
+      INTEGER  ::   ii, ii0
+      INTEGER  ::   ij, ij0
+      INTEGER  ::   bounce_method
+      !!----------------------------------------------------------------------
+      !
+      bounced = .false.
+      bounce_method = 2
+      !
+      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 = ii0 - nimpp + 1
+      ij0 = ij0 - njmpp + 1
+      ii  = ii  - nimpp + 1
+      ij  = ij  - njmpp + 1
+      !
+      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
+      ! -------------------------------
+      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
+
+      SELECT CASE ( bounce_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 adjust_to_ground
+
+
+   SUBROUTINE accel(       xi, e1, uvel, uvel0, ax,                &
+      &             berg , yj, e2, vvel, vvel0, ay, dt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE accel  ***
+      !!
+      !! ** Purpose :   compute the iceberg acceleration.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      TYPE(iceberg ), POINTER           ::   berg
+      REAL(wp), INTENT(in   )           ::   xi   , yj      ! berg position in (i,j) referential
+      REAL(wp), INTENT(in   )           ::   uvel , vvel    ! berg velocity [m/s]
+      REAL(wp), INTENT(in   )           ::   uvel0, vvel0   ! initial berg velocity [m/s]
+      REAL(wp), INTENT(  out)           ::   e1, e2         ! horizontal scale factor at (xi,yj)
+      REAL(wp), INTENT(inout)           ::   ax, ay         ! berg acceleration
+      REAL(wp), INTENT(in   )           ::   dt             ! berg time step
+      !
+      REAL(wp), PARAMETER ::   alpha     = 0._wp      !
+      REAL(wp), PARAMETER ::   beta      = 1._wp      !
+      REAL(wp), PARAMETER ::   vel_lim   =15._wp      ! max allowed berg speed
+      REAL(wp), PARAMETER ::   accel_lim = 1.e-2_wp   ! max allowed berg acceleration
+      REAL(wp), PARAMETER ::   Cr0       = 0.06_wp    !
+      !
+      INTEGER  ::   itloop
+      REAL(wp) ::   uo, vo, ui, vi, ua, va, uwave, vwave, ssh_x, ssh_y, sst, cn, hi
+      REAL(wp) ::   zff, T, D, W, L, M, F
+      REAL(wp) ::   drag_ocn, drag_atm, drag_ice, wave_rad
+      REAL(wp) ::   c_ocn, c_atm, c_ice
+      REAL(wp) ::   ampl, wmod, Cr, Lwavelength, Lcutoff, Ltop
+      REAL(wp) ::   lambda, detA, A11, A12, axe, aye, D_hi
+      REAL(wp) ::   uveln, vveln, us, vs, speed, loc_dx, new_speed
+      !!----------------------------------------------------------------------
+
+      ! Interpolate gridded fields to berg
+      knberg = berg%number(1)
+      CALL interp_flds( xi, e1, uo, ui, ua, ssh_x,                     &
+         &              yj, e2, vo, vi, va, ssh_y, sst, cn, hi, zff )
+
+      M = berg%current_point%mass
+      T = berg%current_point%thickness                           ! total thickness
+      D = ( rn_rho_bergs / rho_seawater ) * T   ! draught (keel depth)
+      F = T - D                                    ! freeboard
+      W = berg%current_point%width
+      L = berg%current_point%length
+
+      hi   = MIN( hi   , D    )
+      D_hi = MAX( 0._wp, D-hi )
+
+      ! Wave radiation
+      uwave = ua - uo   ;   vwave = va - vo     ! Use wind speed rel. to ocean for wave model
+      wmod  = uwave*uwave + vwave*vwave         ! The wave amplitude and length depend on the  current;
+      !                                         ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
+      ampl        = 0.5 * 0.02025 * wmod        ! This is "a", the wave amplitude
+      Lwavelength =       0.32    * wmod        ! Surface wave length fitted to data in table at
+      !                                         ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
+      Lcutoff     = 0.125 * Lwavelength
+      Ltop        = 0.25  * Lwavelength
+      Cr          = Cr0 * MIN(  MAX( 0., (L-Lcutoff) / ((Ltop-Lcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
+      !                                         ! fitted to graph from Carrieres et al.,  POAC Drift Model.
+      wave_rad    = 0.5 * rho_seawater / M * Cr * grav * ampl * MIN( ampl,F ) * (2.*W*L) / (W+L)
+      wmod        = SQRT( ua*ua + va*va )       ! Wind speed
+      IF( wmod /= 0._wp ) THEN
+         uwave = ua/wmod ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
+         vwave = va/wmod
+      ELSE
+         uwave = 0.   ;    vwave=0.   ;    wave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
+      ENDIF
+
+      ! Weighted drag coefficients
+      c_ocn = rho_seawater / M * (0.5*Cd_wv*W*(D_hi)+Cd_wh*W*L)
+      c_atm = rho_air      / M * (0.5*Cd_av*W*F     +Cd_ah*W*L)
+      c_ice = rho_ice      / M * (0.5*Cd_iv*W*hi              )
+      IF( abs(ui) + abs(vi) == 0._wp )   c_ice = 0._wp
+
+      uveln = uvel   ;   vveln = vvel ! Copy starting uvel, vvel
+      !
+      DO itloop = 1, 2  ! Iterate on drag coefficients
+         !
+         us = 0.5 * ( uveln + uvel )
+         vs = 0.5 * ( vveln + vvel )
+         drag_ocn = c_ocn * SQRT( (us-uo)*(us-uo) + (vs-vo)*(vs-vo) )
+         drag_atm = c_atm * SQRT( (us-ua)*(us-ua) + (vs-va)*(vs-va) )
+         drag_ice = c_ice * SQRT( (us-ui)*(us-ui) + (vs-vi)*(vs-vi) )
+         !
+         ! Explicit accelerations
+         !axe= zff*vvel -grav*ssh_x +wave_rad*uwave &
+         !    -drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)
+         !aye=-zff*uvel -grav*ssh_y +wave_rad*vwave &
+         !    -drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)
+         axe = -grav * ssh_x + wave_rad * uwave
+         aye = -grav * ssh_y + wave_rad * vwave
+         IF( alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
+            axe = axe + zff*vvel0
+            aye = aye - zff*uvel0
+         else
+            axe=axe+zff*vvel
+            aye=aye-zff*uvel
+         endif
+         if (beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest
+            axe=axe-drag_ocn*(uvel0-uo) -drag_atm*(uvel0-ua) -drag_ice*(uvel0-ui)
+            aye=aye-drag_ocn*(vvel0-vo) -drag_atm*(vvel0-va) -drag_ice*(vvel0-vi)
+         else
+            axe=axe-drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)
+            aye=aye-drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)
+         endif
+
+         ! Solve for implicit accelerations
+         IF( alpha + beta > 0._wp ) THEN
+            lambda = drag_ocn + drag_atm + drag_ice
+            A11 = 1.+beta*dt*lambda
+            A12 = alpha*dt*zff
+            detA = 1._wp / ( A11*A11 + A12*A12 )
+            ax = detA * ( A11*axe+A12*aye )
+            ay = detA * ( A11*aye-A12*axe )
+         else
+            ax=axe   ;   ay=aye
+         endif
+
+         uveln = uvel0 + dt*ax
+         vveln = vvel0 + dt*ay
+         !
+      END DO      ! itloop
+
+      IF( rn_speed_limit > 0.) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
+         speed = SQRT( uveln*uveln + vveln*vveln )    ! Speed of berg
+         IF( speed > 0.) THEN
+            loc_dx = MIN( e1, e2 )                          ! minimum grid spacing
+            new_speed = loc_dx / dt * rn_speed_limit     ! Speed limit as a factor of dx / dt
+            IF( new_speed < speed ) THEN
+               uveln = uveln * ( new_speed / speed )        ! Scale velocity to reduce speed
+               vveln = vveln * ( new_speed / speed )        ! without changing the direction
+               CALL speed_budget()
+            ENDIF
+         ENDIF
+      ENDIF
+      !                                      ! check the speed and acceleration limits
+      IF( ABS( uveln ) > vel_lim   .OR. ABS( vveln ) > vel_lim   )   &
+         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
+      IF( ABS( ax    ) > accel_lim .OR. ABS( ay    ) > accel_lim )   &
+         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
+      !
+   END SUBROUTINE accel
+
+END MODULE icbdyn
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 3339)
@@ -0,0 +1,490 @@
+MODULE icbini
+
+   !!======================================================================
+   !!                       ***  MODULE  icbini  ***
+   !! Ocean physics:  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
+   !!            -    !                            Restore restart routines
+   !!            -    !  2011-05  (Alderson)       generate_test_icebergs restored
+   !!            -    !  2011-05  (Alderson)       new forcing arrays with extra halo
+   !!            -    !  2011-05  (Alderson)       north fold exchange arrays added
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icb_init      : initialise
+   !!   icb_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
+   USE iom
+   USE fldread
+   USE lbclnk
+
+   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
+
+   CHARACTER(len=100), PRIVATE ::   cn_dir       = './'    ! Root directory for location of icb files
+   TYPE(FLD_N), PRIVATE        ::   sn_icb                 ! information about the calving file to be read
+
+   PRIVATE  icb_alloc      ! routine called in xxx.F90 module
+   PUBLIC   icb_init       ! routine called in xxx.F90 module
+   PUBLIC   icb_gen        ! routine called in xxx.F90 module
+   PUBLIC   icb_nam        ! routine called in xxx.F90 module
+
+CONTAINS
+
+   INTEGER FUNCTION icb_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_alloc  ***
+      !!----------------------------------------------------------------------
+      !
+      INTEGER                              ::   ill
+      !
+      ! 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 )
+
+      icb_alloc = 0
+      ALLOCATE(berg_grid, STAT=ill)
+      icb_alloc = icb_alloc + ill
+      !
+      ALLOCATE( berg_grid%calving      (jpi,jpj)   , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( berg_grid%calving_hflx (jpi,jpj)   , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( berg_grid%stored_heat  (jpi,jpj)   , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( berg_grid%floating_melt(jpi,jpj)   , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( berg_grid%maxclass(jpi,jpj)        , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      !
+      ALLOCATE( berg_grid%stored_ice   (jpi,jpj,nclasses) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      !
+      ALLOCATE( 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) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( vo_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( ff_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( ua_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( va_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+#if defined key_lim2 || defined key_lim3
+      ALLOCATE( ui_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( vi_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+#endif
+      ALLOCATE( ssh_e(0:jpi+1,0:jpj+1) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( initial_width    (nclasses) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( initial_length   (nclasses) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( p_calving(jpi,jpj) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( p_calving_hflx(jpi,jpj) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+
+      ALLOCATE( icbfldpts(jpi) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( icbflddest(jpi) , STAT=ill)
+      icb_alloc = icb_alloc + ill
+      ALLOCATE( icbfldproc(jpni) , 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
+
+   SUBROUTINE icb_init( dt, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_init  ***
+      !!
+      !! ** Purpose :   iceberg initialization.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      REAL(wp)                , INTENT(in) ::   dt      ! iceberg time-step (rdt*nn_fsbc)
+      INTEGER                 , INTENT(in) ::   kt      ! time step number
+      !
+      INTEGER                              ::   ji, jj, jn, i1, j1
+      INTEGER                              ::   ii, inum, ivar, k1
+      INTEGER                              ::   istat1, istat2, istat3
+      CHARACTER(len=300)                   ::   sdist
+      !!----------------------------------------------------------------------
+
+      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' )
+
+      ! set parameters (mostly from namelist)
+      !
+      berg_dt                  = dt
+      initial_width (:) = SQRT(  rn_initial_mass(:) / ( rn_LoW_ratio*rn_rho_bergs*rn_initial_thickness(:) )  )
+      initial_length(:) = rn_LoW_ratio * initial_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
+      p_calving        (:,:)   = 0._wp
+      p_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')
+
+      ! 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 p_calving arrays for this
+      !
+      ! pack i and j together using a scaling of a power of 10
+      icbpack = 10000
+      IF( jpiglo >= icbpack ) CALL ctl_stop('icbini: processor index packing failure')
+      icbfldproc(:) = -1
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            p_calving_hflx(ji,jj) = narea
+            p_calving(ji,jj)      = icbpack*(njmpp+jj-1) + nimpp+ji-1
+         ENDDO
+      ENDDO
+      CALL lbc_lnk( p_calving_hflx, 'T', 1._wp )
+      CALL lbc_lnk( p_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
+      icbdi = -1
+      icbei = -1
+      DO ji = 1,jpi
+         k1 = INT( p_calving(ji,jj) )
+         j1 = INT( k1/icbpack )
+         i1 = k1 - j1*icbpack
+         IF( i1 == nimpp+ji-1 ) THEN
+            IF( icbdi < 0 ) THEN
+               icbdi = ji
+            ELSE
+               icbei = ji
+            ENDIF
+         ENDIF
+      END DO
+     
+      ! repeat for j direction
+      ji = jpi/2
+      icbdj = -1
+      icbej = -1
+      DO jj = 1,jpj
+         k1 = INT( p_calving(ji,jj) )
+         j1 = INT( k1/icbpack )
+         i1 = k1 - j1*icbpack
+         IF( j1 == njmpp+jj-1 ) THEN
+            IF( icbdj < 0 ) THEN
+               icbdj = jj
+            ELSE
+               icbej = jj
+            ENDIF
+         ENDIF
+      END DO
+         
+      ! special for east-west boundary exchange we save the destination index
+      i1 = MAX( icbdi-1, 1)
+      k1 = INT( p_calving(i1,jpj/2) )
+      jj = INT( k1/icbpack )
+      icb_left = REAL( k1 - icbpack*jj, wp )
+      i1 = MIN( icbei+1, jpi )
+      k1 = INT( p_calving(i1,jpj/2) )
+      jj = INT( k1/icbpack )
+      icb_right = REAL( k1 - icbpack*jj, wp )
+      
+      ! north fold
+      IF( npolj > 0 ) THEN
+         !
+         ! icebergs in row icbej+1 get passed across fold
+         icbfldpts(:)  = INT( p_calving(:,icbej+1) )
+         icbflddest(:) = INT( p_calving_hflx(:,icbej+1) )
+         !
+         ! work out list of unique processors to talk to
+         DO ji = icbdi, icbei
+            ii = icbflddest(ji)
+            DO jn = 1, jpni
+               IF( icbfldproc(jn) == -1 ) THEN
+                  icbfldproc(jn) = ii
+                  EXIT
+               ENDIF
+               IF( icbfldproc(jn) == ii ) EXIT
+            ENDDO
+         ENDDO
+      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 ', icbdi, icbei
+         WRITE(numicb,*) 'berg j interior ', icbdj, icbej
+         WRITE(numicb,*) 'berg left       ', icb_left
+         WRITE(numicb,*) 'berg right      ', icb_right
+         jj = jpj/2
+         WRITE(numicb,*) "central j line:"
+         WRITE(numicb,*) "i processor"
+         WRITE(numicb,*) (INT(p_calving_hflx(ji,jj)), ji=1,jpi)
+         WRITE(numicb,*) "i point"
+         WRITE(numicb,*) (INT(p_calving(ji,jj)), ji=1,jpi)
+         ji = jpi/2
+         WRITE(numicb,*) "central i line:"
+         WRITE(numicb,*) "j processor"
+         WRITE(numicb,*) (INT(p_calving_hflx(ji,jj)), jj=1,jpj)
+         WRITE(numicb,*) "j point"
+         WRITE(numicb,*) (INT(p_calving(ji,jj)), jj=1,jpj)
+         IF( npolj > 0 ) THEN
+            WRITE(numicb,*) 'north fold destination points '
+            WRITE(numicb,*) icbfldpts
+            WRITE(numicb,*) 'north fold destination procs  '
+            WRITE(numicb,*) icbflddest
+         ENDIF
+         CALL flush(numicb)
+      ENDIF
+      
+      p_calving(:,:)      = 0._wp
+      p_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
+      kount_bergs(:) = 0
+      kount_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
+         sdist = TRIM( cn_dir )//TRIM( sn_icb%clname )
+         CALL iom_open ( sdist, inum )                              ! open file
+         ivar = iom_varid( inum, 'maxclass', ldstop=.FALSE. )
+         IF( ivar > 0 ) THEN
+            CALL iom_get  ( inum, jpdom_data, 'maxclass', p_calving )   ! read the max distribution array
+            berg_grid%maxclass(:,:) = INT( p_calving )
+            p_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_gen()
+      ELSE
+         IF( nn_test_icebergs > 0 ) THEN
+            CALL icb_gen()
+         ELSE
+            !
+            CALL icebergs_read_restart()
+            restarted_bergs = .TRUE.
+         ENDIF
+      ENDIF
+      !
+      IF( nn_sample_rate .GT. 0 ) CALL traj_init( nitend )
+      !
+      CALL icb_budget_init()
+      !
+      IF( nn_verbose_level >= 2 )   CALL print_bergs('icb_init, initial status', nit000-1)
+      !
+   END SUBROUTINE icb_init
+
+   SUBROUTINE icb_gen()
+
+      ! Local variables
+      INTEGER                         ::   i, j, nbergs
+      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 j = icbdj,icbej
+         DO i = icbdi,icbei
+            IF (tmask(i,j,1) .GT. 0._wp .AND. &
+                gphit(i,j) .GT. rn_test_box(3) .AND. gphit(i,j) .LT. rn_test_box(4) .AND. &
+                glamt(i,j) .GT. rn_test_box(1) .AND. glamt(i,j) .LT. rn_test_box(2)) THEN
+               localberg%mass_scaling = rn_mass_scaling(iberg)
+               localpt%xi = REAL( nimpp+i-1, wp )
+               localpt%yj = REAL( njmpp+j-1, wp )
+               localpt%lon = bilin(glamt, localpt%xi, localpt%yj, 'T', 0, 0 )
+               localpt%lat = bilin(gphit, localpt%xi, localpt%yj, 'T', 0, 0 )
+               localpt%mass = rn_initial_mass(iberg)
+               localpt%thickness = rn_initial_thickness(iberg)
+               localpt%width = initial_width(iberg)
+               localpt%length = initial_length(iberg)
+               localpt%year = iyr
+               localpt%day = FLOAT(iday)+(FLOAT(ihr)+FLOAT(imin)/60._wp)/24._wp
+               localpt%mass_of_bits = 0._wp
+               localpt%heat_density = 0._wp
+               localpt%uvel = 0._wp
+               localpt%vvel = 0._wp
+               CALL increment_kounter()
+               localberg%number(:) = kount_bergs(:)
+               call add_new_berg_to_list(localberg, localpt)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      nbergs = count_bergs()
+      IF( lk_mpp ) CALL mpp_sum(nbergs)
+      WRITE(numicb,'(a,i6,a)') 'diamonds, icb_gen: ',nbergs,' were generated'
+
+   END SUBROUTINE icb_gen
+
+   SUBROUTINE icb_nam
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE icb_nam  ***
+      !!
+      !! ** Purpose :   read domaine namelists and print the variables.
+      !!
+      !! ** input   : - namberg namelist
+      !!----------------------------------------------------------------------
+
+      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
+      INTEGER   ::   k  ! loop counter
+      REAL(wp)  ::   dfact
+      !!----------------------------------------------------------------------
+
+      ! (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 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
+         WRITE(numout,*) '~~~~~~~~~~~~~ '
+         WRITE(numout,*) 'NO icebergs used'
+         RETURN
+      ENDIF
+
+      IF( nn_test_icebergs > nclasses ) THEN
+          IF( lwp ) WRITE(numout,*) 'Resetting nn_test_icebergs to ',nclasses
+          nn_test_icebergs = nclasses
+      ENDIF
+
+      dfact = SUM( rn_distribution )
+      IF( dfact < 1._wp ) THEN
+         IF( dfact <= 0._wp ) THEN
+            CALL ctl_stop('icb_init: sum of berg distribution equal to zero')
+         ELSE
+            rn_distribution(:) = rn_distribution(:) / dfact
+            CALL ctl_warn('icb_init: 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 k=1,nclasses
+            WRITE(numout,'(a,f15.2)') '                                                                ',rn_initial_mass(k)
+         ENDDO
+         WRITE(numout,*) '   Fraction of calving to apply to this class (non-dim)         rn_distribution     ='
+         DO k=1,nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_distribution(k)
+         ENDDO
+         WRITE(numout,*) '   Ratio between effective and real iceberg mass (non-dim)      rn_mass_scaling     = '
+         DO k=1,nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_mass_scaling(k)
+         ENDDO
+         WRITE(numout,*) '   Total thickness of newly calved bergs (m)                   rn_initial_thickness = '
+         DO k=1,nclasses
+            WRITE(numout,'(a,f10.2)') '                                                                ',rn_initial_thickness(k)
+         ENDDO
+         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)  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_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 3339)
@@ -0,0 +1,868 @@
+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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   mpp_send_bergs   :  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 (icbdi-icbei, icbdj-icbej)
+   !!                       so that iceberg does not exist in more than one processor
+   !!                       North fold exchanges controlled by three arrays:
+   !!                          icbflddest - unique processor numbers that current one exchanges with
+   !!                          icbfldproc - processor number that current grid point exchanges with
+   !!                          icbfldpts  - 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, PUBLIC      ::   delta_buf = 25             ! Size by which to increment buffers
+   INTEGER, PARAMETER, PUBLIC      ::   buffer_width = 15+nkounts  ! items to store for each berg
+
+#endif
+
+   PUBLIC   lbc_send_bergs
+   PRIVATE  lbc_nfld_bergs
+   PUBLIC   mpp_send_bergs
+   PUBLIC   dealloc_buffers
+
+#if defined key_mpp_mpi
+   PRIVATE  mpp_nfld_bergs
+   PRIVATE  dealloc_buffer
+   PRIVATE  pack_berg_into_buffer
+   PRIVATE  unpack_berg_from_buffer
+   PRIVATE  increase_buffer
+   PRIVATE  increase_ibuffer
+#endif
+
+   !!-------------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lbc_send_bergs()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE lbc_send_bergs  ***
+      !!
+      !! ** Purpose :   in non-mpp case need to deal with cyclic conditions
+      !!                including north-fold
+      !!----------------------------------------------------------------------
+      ! Local variables
+      TYPE(iceberg), POINTER :: this
+      TYPE(point)  , POINTER :: pt
+      INTEGER                :: ine
+
+      !! 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
+            ine = INT( pt%xi + 0.5 )
+            IF( ine .GT. nimpp+icbei-1 ) THEN
+               pt%xi = icb_right + MOD(pt%xi, 1._wp ) - 1._wp
+            ELSE IF( ine .LT. nimpp+icbdi-1 ) THEN
+               pt%xi = icb_left + MOD(pt%xi, 1._wp )
+            ENDIF
+            this => this%next
+         ENDDO
+
+      ENDIF
+
+      !! north/south boundaries
+      !! ======================
+
+      !! south symmetric
+
+      IF( nperio == 2 ) THEN
+         CALL ctl_stop(' south symmetric condition not implemented for icebergs')
+      ENDIF
+
+      !! north fold
+
+      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) CALL lbc_nfld_bergs()
+
+   END SUBROUTINE lbc_send_bergs
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE lbc_nfld_bergs()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE lbc_nfld_bergs  ***
+      !!
+      !! ** Purpose :   single processor north fold exchange
+      !!----------------------------------------------------------------------
+      ! 
+      ! Local variables
+      TYPE(iceberg), POINTER :: this
+      TYPE(point)  , POINTER :: pt
+      INTEGER                :: ine, jne, kpts
+      INTEGER                :: jiglo, jjglo
+
+      this => first_berg
+      DO WHILE( ASSOCIATED(this) )
+         pt => this%current_point
+         jne = INT( pt%yj + 0.5 )
+         IF( jne .GT. njmpp+icbej-1 ) THEN
+            !
+            ine = INT( pt%xi + 0.5 )
+            kpts  = icbfldpts (ine-nimpp+1)
+            !
+            ! moving across the cut line means both position and
+            ! velocity must change
+            jjglo = INT( kpts/icbpack )
+            jiglo = kpts - icbpack*jjglo
+            pt%xi = jiglo - ( pt%xi - REAL(ine,wp) )
+            pt%yj = jjglo - ( pt%yj - REAL(jne,wp) )
+            pt%uvel = -1._wp * pt%uvel
+            pt%vvel = -1._wp * pt%vvel
+         ENDIF
+         this => this%next
+      ENDDO
+
+   END SUBROUTINE lbc_nfld_bergs
+
+   !!-------------------------------------------------------------------------
+
+#if defined key_mpp_mpi
+
+   SUBROUTINE mpp_send_bergs()
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE mpp_send_bergs  ***
+      !!
+      !! ** 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
+      !!----------------------------------------------------------------------
+
+      ! Local variables
+      TYPE(iceberg)         , POINTER     :: tmpberg, this
+      TYPE(point)           , POINTER     :: pt
+      INTEGER                             :: nbergs_to_send_e, nbergs_to_send_w
+      INTEGER                             :: nbergs_to_send_n, nbergs_to_send_s
+      INTEGER                             :: nbergs_rcvd_from_e, nbergs_rcvd_from_w
+      INTEGER                             :: nbergs_rcvd_from_n, nbergs_rcvd_from_s
+      INTEGER                             :: i, nbergs_start, nbergs_end
+      INTEGER                             :: ine, jne
+      INTEGER                             :: pe_N, pe_S, pe_W, pe_E
+      REAL(wp), DIMENSION(2)              :: ewbergs, webergs, nsbergs, snbergs
+      INTEGER                             :: ml_req1, ml_req2, ml_req3, ml_req4
+      INTEGER                             :: ml_req5, ml_req6, ml_req7, ml_req8, ml_err
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat
+
+      ! set up indices of neighbouring processors
+      pe_N = -1
+      pe_S = -1
+      pe_W = -1
+      pe_E = -1
+      IF( nbondi .EQ.  0 .OR. nbondi .EQ. 1) pe_W = nowe
+      IF( nbondi .EQ. -1 .OR. nbondi .EQ. 0) pe_E = noea
+      IF( nbondj .EQ.  0 .OR. nbondj .EQ. 1) pe_S = noso
+      IF( nbondj .EQ. -1 .OR. nbondj .EQ. 0) pe_N = nono
+      !
+      ! at northern line of processors with north fold handle bergs differently
+      IF( npolj > 0 ) pe_N = -1
+
+      ! if there's only one processor in x direction then don't let mpp try to handle periodicity
+      IF( jpni == 1 ) THEN
+         pe_E = -1
+         pe_W = -1
+      ENDIF
+
+      IF( nn_verbose_level >= 2 ) THEN
+         WRITE(numicb,*) 'processor west  : ', pe_W
+         WRITE(numicb,*) 'processor east  : ', pe_E
+         WRITE(numicb,*) 'processor north : ', pe_N
+         WRITE(numicb,*) 'processor south : ', pe_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 lbc_send_bergs and called here
+
+      IF( jpni == 1 ) CALL lbc_send_bergs()
+
+      ! 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
+         nbergs_start = count_bergs()
+      ENDIF
+
+      nbergs_to_send_e = 0
+      nbergs_to_send_w = 0
+      nbergs_to_send_n = 0
+      nbergs_to_send_s = 0
+      nbergs_rcvd_from_e = 0
+      nbergs_rcvd_from_w = 0
+      nbergs_rcvd_from_n = 0
+      nbergs_rcvd_from_s = 0
+
+      ! Find number of bergs that headed east/west
+      IF( ASSOCIATED(first_berg) ) THEN
+         this => first_berg
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            ine = INT( pt%xi + 0.5 )
+            IF( pe_E >= 0 .AND. ine .GT. nimpp+icbei-1 ) THEN
+               tmpberg => this
+               this => this%next
+               nbergs_to_send_e = nbergs_to_send_e + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',ktberg,' packing berg ',tmpberg%number(:),' for transfer to east'
+                  CALL flush( numicb )
+               ENDIF
+               ! deal with periodic case
+               tmpberg%current_point%xi = icb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp
+               ! now pack it into buffer and delete from list
+               CALL pack_berg_into_buffer( tmpberg, obuffer_e, nbergs_to_send_e)
+               CALL delete_iceberg_from_list(first_berg, tmpberg)
+            ELSE IF( pe_W >= 0 .AND. ine .LT. nimpp+icbdi-1 ) THEN
+               tmpberg => this
+               this => this%next
+               nbergs_to_send_w = nbergs_to_send_w + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',ktberg,' packing berg ',tmpberg%number(:),' for transfer to west'
+                  CALL flush( numicb )
+               ENDIF
+               ! deal with periodic case
+               tmpberg%current_point%xi = icb_left + MOD(tmpberg%current_point%xi, 1._wp )
+               ! now pack it into buffer and delete from list
+               CALL pack_berg_into_buffer( tmpberg, obuffer_w, nbergs_to_send_w)
+               CALL delete_iceberg_from_list(first_berg, tmpberg)
+            ELSE
+               this => this%next
+            ENDIF
+         ENDDO
+      ENDIF
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',ktberg,' send ew: ', nbergs_to_send_e, nbergs_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 )
+         webergs(1) = nbergs_to_send_e
+         CALL mppsend( 12, webergs(1), 1, pe_E, ml_req1)
+         CALL mpprecv( 11, ewbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err )
+         nbergs_rcvd_from_e = INT( ewbergs(2) )
+      CASE(  0 )
+         ewbergs(1) = nbergs_to_send_w
+         webergs(1) = nbergs_to_send_e
+         CALL mppsend( 11, ewbergs(1), 1, pe_W, ml_req2)
+         CALL mppsend( 12, webergs(1), 1, pe_E, ml_req3)
+         CALL mpprecv( 11, ewbergs(2), 1 )
+         CALL mpprecv( 12, webergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err )
+         IF( l_isend ) CALL mpi_wait( ml_req3, ml_stat, ml_err )
+         nbergs_rcvd_from_e = INT( ewbergs(2) )
+         nbergs_rcvd_from_w = INT( webergs(2) )
+      CASE(  1 )
+         ewbergs(1) = nbergs_to_send_w
+         CALL mppsend( 11, ewbergs(1), 1, pe_W, ml_req4)
+         CALL mpprecv( 12, webergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req4, ml_stat, ml_err )
+         nbergs_rcvd_from_w = INT( webergs(2) )
+      END SELECT
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',ktberg,' recv ew: ', nbergs_rcvd_from_w, nbergs_rcvd_from_e
+         call flush(numicb)
+      endif
+
+      SELECT CASE ( nbondi )
+      CASE( -1 )
+         IF( nbergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, nbergs_to_send_e*buffer_width, pe_E, ml_req1 )
+         IF( nbergs_rcvd_from_e > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_e, nbergs_rcvd_from_e)
+            CALL mpprecv( 13, ibuffer_e%data, nbergs_rcvd_from_e*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_e
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_e, i)
+         ENDDO
+      CASE(  0 )
+         IF( nbergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, nbergs_to_send_w*buffer_width, pe_W, ml_req2 )
+         IF( nbergs_to_send_e > 0 ) CALL mppsend( 14, obuffer_e%data, nbergs_to_send_e*buffer_width, pe_E, ml_req3 )
+         IF( nbergs_rcvd_from_e > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_e, nbergs_rcvd_from_e)
+            CALL mpprecv( 13, ibuffer_e%data, nbergs_rcvd_from_e*buffer_width )
+         ENDIF
+         IF( nbergs_rcvd_from_w > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_w, nbergs_rcvd_from_w)
+            CALL mpprecv( 14, ibuffer_w%data, nbergs_rcvd_from_w*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err )
+         IF( nbergs_to_send_e > 0 .AND. l_isend ) CALL mpi_wait( ml_req3, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_e
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_e%data(16,i)),' from east'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_e, i)
+         ENDDO
+         DO i = 1, nbergs_rcvd_from_w
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_w, i)
+         ENDDO
+      CASE(  1 )
+         IF( nbergs_to_send_w > 0 ) CALL mppsend( 13, obuffer_w%data, nbergs_to_send_w*buffer_width, pe_W, ml_req4 )
+         IF( nbergs_rcvd_from_w > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_w, nbergs_rcvd_from_w)
+            CALL mpprecv( 14, ibuffer_w%data, nbergs_rcvd_from_w*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_w > 0 .AND. l_isend ) CALL mpi_wait( ml_req4, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_w
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_w%data(16,i)),' from west'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_w, i)
+         ENDDO
+      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
+            jne = INT( pt%yj + 0.5 )
+            IF( pe_N >= 0 .AND. jne .GT. njmpp+icbej-1 ) THEN
+               tmpberg => this
+               this => this%next
+               nbergs_to_send_n = nbergs_to_send_n + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',ktberg,' packing berg ',tmpberg%number(:),' for transfer to north'
+                  CALL flush( numicb )
+               ENDIF
+               CALL pack_berg_into_buffer( tmpberg, obuffer_n, nbergs_to_send_n)
+               CALL delete_iceberg_from_list(first_berg, tmpberg)
+            ELSE IF( pe_S >= 0 .AND. jne .LT. njmpp+icbdj-1 ) THEN
+               tmpberg => this
+               this => this%next
+               nbergs_to_send_s = nbergs_to_send_s + 1
+               IF( nn_verbose_level >= 4 ) THEN
+                  WRITE(numicb,*) 'bergstep ',ktberg,' packing berg ',tmpberg%number(:),' for transfer to south'
+                  CALL flush( numicb )
+               ENDIF
+               CALL pack_berg_into_buffer( tmpberg, obuffer_s, nbergs_to_send_s)
+               CALL delete_iceberg_from_list(first_berg, tmpberg)
+            ELSE
+               this => this%next
+            ENDIF
+         ENDDO
+      ENDIF
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',ktberg,' send ns: ', nbergs_to_send_n, nbergs_to_send_s
+         call flush(numicb)
+      endif
+
+      ! send bergs north
+      ! and receive bergs from south (ie ones sent north)
+
+      SELECT CASE ( nbondj )
+      CASE( -1 )
+         snbergs(1) = nbergs_to_send_n
+         CALL mppsend( 16, snbergs(1), 1, pe_N, ml_req1)
+         CALL mpprecv( 15, nsbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err )
+         nbergs_rcvd_from_n = INT( nsbergs(2) )
+      CASE(  0 )
+         nsbergs(1) = nbergs_to_send_s
+         snbergs(1) = nbergs_to_send_n
+         CALL mppsend( 15, nsbergs(1), 1, pe_S, ml_req2)
+         CALL mppsend( 16, snbergs(1), 1, pe_N, ml_req3)
+         CALL mpprecv( 15, nsbergs(2), 1 )
+         CALL mpprecv( 16, snbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err )
+         IF( l_isend ) CALL mpi_wait( ml_req3, ml_stat, ml_err )
+         nbergs_rcvd_from_n = INT( nsbergs(2) )
+         nbergs_rcvd_from_s = INT( snbergs(2) )
+      CASE(  1 )
+         nsbergs(1) = nbergs_to_send_s
+         CALL mppsend( 15, nsbergs(1), 1, pe_S, ml_req4)
+         CALL mpprecv( 16, snbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req4, ml_stat, ml_err )
+         nbergs_rcvd_from_s = INT( snbergs(2) )
+      END SELECT
+      if( nn_verbose_level >= 3) then
+         write(numicb,*) 'bergstep ',ktberg,' recv ns: ', nbergs_rcvd_from_s, nbergs_rcvd_from_n
+         call flush(numicb)
+      endif
+
+      SELECT CASE ( nbondj )
+      CASE( -1 )
+         IF( nbergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, nbergs_to_send_n*buffer_width, pe_N, ml_req1 )
+         IF( nbergs_rcvd_from_n > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_n, nbergs_rcvd_from_n)
+            CALL mpprecv( 17, ibuffer_n%data, nbergs_rcvd_from_n*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_n
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_n, i)
+         ENDDO
+      CASE(  0 )
+         IF( nbergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, nbergs_to_send_s*buffer_width, pe_S, ml_req2 )
+         IF( nbergs_to_send_n > 0 ) CALL mppsend( 18, obuffer_n%data, nbergs_to_send_n*buffer_width, pe_N, ml_req3 )
+         IF( nbergs_rcvd_from_n > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_n, nbergs_rcvd_from_n)
+            CALL mpprecv( 17, ibuffer_n%data, nbergs_rcvd_from_n*buffer_width )
+         ENDIF
+         IF( nbergs_rcvd_from_s > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_s, nbergs_rcvd_from_s)
+            CALL mpprecv( 18, ibuffer_s%data, nbergs_rcvd_from_s*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err )
+         IF( nbergs_to_send_n > 0 .AND. l_isend ) CALL mpi_wait( ml_req3, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_n
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_n%data(16,i)),' from north'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_n, i)
+         ENDDO
+         DO i = 1, nbergs_rcvd_from_s
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_s, i)
+         ENDDO
+      CASE(  1 )
+         IF( nbergs_to_send_s > 0 ) CALL mppsend( 17, obuffer_s%data, nbergs_to_send_s*buffer_width, pe_S, ml_req4 )
+         IF( nbergs_rcvd_from_s > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_s, nbergs_rcvd_from_s)
+            CALL mpprecv( 18, ibuffer_s%data, nbergs_rcvd_from_s*buffer_width )
+         ENDIF
+         IF( nbergs_to_send_s > 0 .AND. l_isend ) CALL mpi_wait( ml_req4, ml_stat, ml_err )
+         DO i = 1, nbergs_rcvd_from_s
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_s%data(16,i)),' from south'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_s, i)
+         ENDDO
+      END SELECT
+
+      IF( nn_verbose_level > 0 ) THEN
+         ! compare the number of icebergs on this processor from the start to the end
+         nbergs_end = count_bergs()
+         i = ( nbergs_rcvd_from_n + nbergs_rcvd_from_s + nbergs_rcvd_from_e + nbergs_rcvd_from_w ) - &
+             ( nbergs_to_send_n + nbergs_to_send_s + nbergs_to_send_e + nbergs_to_send_w )
+         IF( nbergs_end-(nbergs_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: nbergs_end=', &
+                                nbergs_end,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_start=', &
+                                nbergs_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=', &
+                                nbergs_end-(nbergs_start+i),' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_n=', &
+                                nbergs_to_send_n,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_s=', &
+                                nbergs_to_send_s,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_e=', &
+                                nbergs_to_send_e,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_to_send_w=', &
+                                nbergs_to_send_w,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_n=', &
+                                nbergs_rcvd_from_n,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_s=', &
+                                nbergs_rcvd_from_s,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_e=', &
+                                nbergs_rcvd_from_e,' on PE',narea
+            WRITE( numicb,1000) 'send_bergs_to_other_pes: nbergs_rcvd_from_w=', &
+                                nbergs_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 lbc_send_bergs
+      IF( npolj /= 0 .AND. jpni > 1 ) CALL mpp_nfld_bergs( )
+
+      IF( nn_verbose_level > 0 ) THEN
+         i = 0
+         this => first_berg
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            ine = INT( pt%xi + 0.5 )
+            jne = INT( pt%yj + 0.5 )
+!           CALL check_position(grd, this, 'exchange (bot)')
+            IF( ine .LT. nimpp+icbdi-1 .OR. &
+                ine .GT. nimpp+icbei-1 .OR. &
+                jne .LT. njmpp+icbdj-1 .OR. &
+                jne .GT. njmpp+icbej-1) THEN
+               i = i + 1
+               WRITE(numicb,*) 'berg lost in halo: ', this%number(:),ine,jne
+               WRITE(numicb,*) '                   ', nimpp, njmpp
+               WRITE(numicb,*) '                   ', icbdi, icbei, icbdj, icbej
+               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 mpp_send_bergs
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE mpp_nfld_bergs()
+
+      ! Local variables
+      TYPE(iceberg)         , POINTER     :: tmpberg, this
+      TYPE(point)           , POINTER     :: pt
+      INTEGER                             :: nbergs_to_send
+      INTEGER                             :: nbergs_to_rcv
+      INTEGER                             :: jiglo, jjglo, jk, jn
+      INTEGER                             :: jfldproc, kproc, kpts
+      INTEGER                             :: ine, jne
+      REAL(wp), DIMENSION(2)              :: sbergs, nbergs
+      INTEGER                             :: ml_req1, ml_req2, ml_err
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat
+
+      ! set up indices of neighbouring processors
+
+      ! icbfldproc 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( icbfldproc(jn) == -1 ) EXIT
+         jfldproc = icbfldproc(jn)
+         nbergs_to_send = 0
+
+         ! Find number of bergs that need to be exchanged
+         ! Pick out exchanges with processor jfldproc
+         ! if jfldproc is this processor then don't send
+         !
+         IF( ASSOCIATED(first_berg) ) THEN
+            this => first_berg
+            DO WHILE (ASSOCIATED(this))
+               pt => this%current_point
+               ine = INT( pt%xi + 0.5 )
+               jne = INT( pt%yj + 0.5 )
+               kpts  = icbfldpts (ine-nimpp+1)
+               kproc = icbflddest(ine-nimpp+1)
+               IF( jne .GT. njmpp+icbej-1 ) THEN
+                  IF( kproc == jfldproc ) THEN
+                     !
+                     ! moving across the cut line means both position and
+                     ! velocity must change
+                     jjglo = INT( kpts/icbpack )
+                     jiglo = kpts - icbpack*jjglo
+                     pt%xi = jiglo - ( pt%xi - REAL(ine,wp) )
+                     pt%yj = jjglo - ( pt%yj - REAL(jne,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( kproc /= narea ) THEN
+                        tmpberg => this
+                        nbergs_to_send = nbergs_to_send + 1
+                        IF( nn_verbose_level >= 4 ) THEN
+                           WRITE(numicb,*) 'bergstep ',ktberg,' packing berg ',tmpberg%number(:),' for north fold'
+                           CALL flush( numicb )
+                        ENDIF
+                        CALL pack_berg_into_buffer( tmpberg, obuffer_f, nbergs_to_send)
+                        CALL delete_iceberg_from_list(first_berg, tmpberg)
+                     ENDIF
+                     !
+                  ENDIF
+               ENDIF
+               this => this%next
+            ENDDO
+         ENDIF
+         if( nn_verbose_level >= 3) then
+            write(numicb,*) 'bergstep ',ktberg,' send nfld: ', nbergs_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( jfldproc == narea ) CYCLE
+
+         sbergs(1) = nbergs_to_send
+         CALL mppsend( 21, sbergs(1), 1, jfldproc-1, ml_req1)
+         CALL mpprecv( 21, nbergs(2), 1 )
+         IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err )
+         nbergs_to_rcv = INT( nbergs(2) )
+
+         ! send bergs
+
+         IF( nbergs_to_send > 0 )  &
+             CALL mppsend( 12, obuffer_f%data, nbergs_to_send*buffer_width, jfldproc-1, ml_req2 )
+         IF( nbergs_to_rcv  > 0 ) THEN
+            CALL increase_ibuffer(ibuffer_f, nbergs_to_rcv)
+            CALL mpprecv( 12, ibuffer_f%data, nbergs_to_rcv*buffer_width )
+         ENDIF
+         IF( nbergs_to_send > 0 .AND. l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err )
+         DO jk = 1, nbergs_to_rcv
+            IF( nn_verbose_level >= 4 ) THEN
+               WRITE(numicb,*) 'bergstep ',ktberg,' unpacking berg ',INT(ibuffer_f%data(16,jk)),' from north fold'
+               CALL flush( numicb )
+            ENDIF
+            CALL unpack_berg_from_buffer(first_berg, ibuffer_f, jk )
+         ENDDO
+
+      ENDDO
+
+   END SUBROUTINE mpp_nfld_bergs
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE dealloc_buffers()
+
+      CALL dealloc_buffer( obuffer_n )
+      CALL dealloc_buffer( obuffer_s )
+      CALL dealloc_buffer( obuffer_e )
+      CALL dealloc_buffer( obuffer_w )
+      CALL dealloc_buffer( ibuffer_n )
+      CALL dealloc_buffer( ibuffer_s )
+      CALL dealloc_buffer( ibuffer_e )
+      CALL dealloc_buffer( ibuffer_w )
+
+   END SUBROUTINE dealloc_buffers
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE dealloc_buffer(buff)
+      ! Arguments
+      TYPE(buffer), POINTER :: buff
+
+      IF( ASSOCIATED(buff) ) THEN
+         IF( ASSOCIATED(buff%data)) DEALLOCATE(buff%data)
+         DEALLOCATE(buff)
+      ENDIF
+   END SUBROUTINE dealloc_buffer
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE pack_berg_into_buffer(berg, buff, n)
+      ! Arguments
+      TYPE(iceberg),            POINTER :: berg
+      TYPE(buffer) ,            POINTER :: buff
+      INTEGER      , INTENT(in)         :: n
+      ! Local
+      INTEGER                           :: k
+
+      IF( .NOT. ASSOCIATED(buff) ) CALL increase_buffer( buff, delta_buf )
+      IF( n .GT. buff%size ) CALL increase_buffer( buff, delta_buf )
+
+      !! pack points into buffer
+
+      buff%data( 1,n) = berg%current_point%lon
+      buff%data( 2,n) = berg%current_point%lat
+      buff%data( 3,n) = berg%current_point%uvel
+      buff%data( 4,n) = berg%current_point%vvel
+      buff%data( 5,n) = berg%current_point%xi
+      buff%data( 6,n) = berg%current_point%yj
+      buff%data( 7,n) = float(berg%current_point%year)
+      buff%data( 8,n) = berg%current_point%day
+      buff%data( 9,n) = berg%current_point%mass
+      buff%data(10,n) = berg%current_point%thickness
+      buff%data(11,n) = berg%current_point%width
+      buff%data(12,n) = berg%current_point%length
+      buff%data(13,n) = berg%current_point%mass_of_bits
+      buff%data(14,n) = berg%current_point%heat_density
+
+      buff%data(15,n) = berg%mass_scaling
+      DO k=1,nkounts
+         buff%data(15+k,n) = REAL( berg%number(k), wp )
+      END DO
+
+   END SUBROUTINE pack_berg_into_buffer
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE unpack_berg_from_buffer(first, buff, n)
+      ! Arguments
+      TYPE(iceberg),             POINTER :: first
+      TYPE(buffer) ,             POINTER :: buff
+      INTEGER      , INTENT(in)          :: n
+      ! Local variables
+      LOGICAL                            :: lres
+      TYPE(iceberg)                      :: currentberg
+      TYPE(point)                        :: pt
+      INTEGER                           :: k
+
+      pt%lon            =      buff%data( 1,n)
+      pt%lat            =      buff%data( 2,n)
+      pt%uvel           =      buff%data( 3,n)
+      pt%vvel           =      buff%data( 4,n)
+      pt%xi             =      buff%data( 5,n)
+      pt%yj             =      buff%data( 6,n)
+      pt%year           = INT( buff%data( 7,n) )
+      pt%day            =      buff%data( 8,n)
+      pt%mass           =      buff%data( 9,n)
+      pt%thickness      =      buff%data(10,n)
+      pt%width          =      buff%data(11,n)
+      pt%length         =      buff%data(12,n)
+      pt%mass_of_bits   =      buff%data(13,n)
+      pt%heat_density   =      buff%data(14,n)
+
+      currentberg%mass_scaling =      buff%data(15,n)
+      DO k=1,nkounts
+         currentberg%number(k) = INT( buff%data(15+k,n) )
+      END DO
+
+      CALL add_new_berg_to_list(currentberg, pt )
+
+   END SUBROUTINE unpack_berg_from_buffer
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE increase_buffer(old,delta)
+      ! Arguments
+      TYPE(buffer),             POINTER :: old
+      INTEGER     , INTENT(in)          :: delta
+      ! Local variables
+      TYPE(buffer),             POINTER :: new
+      INTEGER                           :: new_size
+
+      IF( .NOT. ASSOCIATED(old) ) THEN
+         new_size = delta
+      ELSE
+         new_size = old%size + delta
+      ENDIF
+      ALLOCATE( new )
+      ALLOCATE( new%data( buffer_width, new_size) )
+      new%size = new_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 increase_buffer
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE increase_ibuffer(old,delta)
+      ! Arguments
+      TYPE(buffer),            POINTER :: old
+      INTEGER     , INTENT(in)         :: delta
+      ! Local variables
+      TYPE(buffer),            POINTER :: new
+      INTEGER                          :: new_size, old_size
+
+      IF( .NOT. ASSOCIATED(old) ) THEN
+         new_size = delta + delta_buf
+         old_size = 0
+      ELSE
+         old_size = old%size
+         IF( delta .LT. old%size ) THEN
+            new_size = old%size + delta
+         ELSE
+            new_size = delta + delta_buf
+         ENDIF
+      ENDIF
+
+      IF( old_size .NE. new_size ) THEN
+         ALLOCATE( new )
+         ALLOCATE( new%data( buffer_width, new_size) )
+         new%size = new_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,*) 'increase_ibuffer',narea,' increased to',new_size
+      ENDIF
+
+   END SUBROUTINE increase_ibuffer
+
+   !!-------------------------------------------------------------------------
+
+#else
+
+   SUBROUTINE mpp_send_bergs()
+      WRITE(numout,*) 'mpp_send_bergs: You should not have seen this message!!'
+   END SUBROUTINE mpp_send_bergs
+
+   SUBROUTINE dealloc_buffers()
+      WRITE(numout,*) 'dealloc_buffers: You should not have seen this message!!'
+   END SUBROUTINE dealloc_buffers
+
+#endif
+
+END MODULE icblbc
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 3339)
@@ -0,0 +1,416 @@
+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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   icebergs_read_restart    : initialise
+   !!   icebergs_write_restart   : generate test icebergs
+   !!----------------------------------------------------------------------
+   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
+
+   REAL(wp), DIMENSION(:,:,:), PRIVATE, POINTER :: griddata => NULL()    ! need 2d array to read in with
+   INTEGER,  DIMENSION(3),     PRIVATE          :: strt3, lngth3
+   INTEGER,                    PRIVATE          :: lonid, latid, xid, yid, uvelid, vvelid
+   INTEGER,                    PRIVATE          :: massid, thicknessid, widthid, lengthid
+   INTEGER,                    PRIVATE          :: yearid, dayid
+   INTEGER,                    PRIVATE          :: scaling_id, mass_of_bits_id, heat_density_id, numberid
+   INTEGER,                    PRIVATE          :: siceid, sheatid, calvid, calvhid, kountid
+
+   PUBLIC   icebergs_read_restart       ! routine called in xxx.F90 module
+   PUBLIC   icebergs_write_restart      ! routine called in xxx.F90 module
+
+CONTAINS
+
+   SUBROUTINE icebergs_read_restart()
+      ! sga - read a restart file
+      !       for this version, we just read back in the restart for this processor
+      !       so we cannot change the processor layout currently with iceberg code
+
+      ! Local variables
+      INTEGER                               :: ndim, nvar, natt
+      INTEGER                               :: k, iret, ncid, unlim_dim, x_dim, y_dim, c_dim, nbergs_in_file
+      INTEGER                               :: mclass
+      INTEGER, DIMENSION(1)                 :: strt, lngth, idata
+      INTEGER, DIMENSION(2)                 :: strt2, lngth2
+      INTEGER, DIMENSION(nkounts)           :: idata2
+      REAL(wp), DIMENSION(1)                :: data                                         ! need 1d array to read in with
+                                                                                            ! start and count arrays
+      LOGICAL                               :: lres, found_restart, multiPErestart=.FALSE.
+      CHARACTER(len=80)                     :: filename
+      CHARACTER(len=NF90_MAX_NAME)          :: dname
+      TYPE(iceberg)                         :: localberg ! NOT a pointer but an actual local variable
+      TYPE(point)                           :: localpt   ! NOT a pointer but an actual local variable
+
+      IF( .NOT. ASSOCIATED(griddata) ) ALLOCATE( griddata(jpi,jpj,1) )
+
+      ! Find a restart file
+      multiPErestart=.FALSE.
+      DO
+         filename = ' '
+         filename = 'restart_icebergs.nc'
+         INQUIRE(file=TRIM(filename),exist=found_restart)
+         IF (found_restart) EXIT
+         filename = ' '
+         WRITE(filename,'("restart_icebergs_",I4.4,".nc")') narea-1
+         INQUIRE(file=TRIM(filename),exist=found_restart)
+         IF (found_restart) THEN
+            multiPErestart=.TRUE.
+            EXIT
+         ENDIF
+         IF (nn_verbose_level >= 0 .AND. lwp) &
+            WRITE(numout,'(a)') 'read_restart_bergs: no restart file found'
+         multiPErestart=.TRUE. ! force sanity checking in a multi-PE mode if no file found on this PE
+         EXIT
+      ENDDO
+
+      IF ( .NOT. 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(filename)
+
+      iret = NF90_OPEN(TRIM(filename), NF90_NOWRITE, ncid)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed')
+
+      iret = nf90_inquire(ncid, ndim, nvar, natt, unlim_dim)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed')
+
+      IF( unlim_dim .NE. -1) THEN
+
+         iret = nf90_inquire_dimension(ncid, unlim_dim, dname, nbergs_in_file)
+         IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed')
+
+         iret = NF90_INQ_VARID(ncid, 'number', numberid)
+         iret = NF90_INQ_VARID(ncid, 'mass_scaling', scaling_id)
+         iret = NF90_INQ_VARID(ncid, 'xi', xid)
+         iret = NF90_INQ_VARID(ncid, 'yj', yid)
+         iret = NF90_INQ_VARID(ncid, 'lon', lonid)
+         iret = NF90_INQ_VARID(ncid, 'lat', latid)
+         iret = NF90_INQ_VARID(ncid, 'uvel', uvelid)
+         iret = NF90_INQ_VARID(ncid, 'vvel', vvelid)
+         iret = NF90_INQ_VARID(ncid, 'mass', massid)
+         iret = NF90_INQ_VARID(ncid, 'thickness', thicknessid)
+         iret = NF90_INQ_VARID(ncid, 'width', widthid)
+         iret = NF90_INQ_VARID(ncid, 'length', lengthid)
+         iret = NF90_INQ_VARID(ncid, 'year', yearid)
+         iret = NF90_INQ_VARID(ncid, 'day', dayid)
+         iret = NF90_INQ_VARID(ncid, 'mass_of_bits', mass_of_bits_id)
+         iret = NF90_INQ_VARID(ncid, 'heat_density', heat_density_id)
+
+         lngth(1) = 1
+         strt2(1) = 1
+         lngth2(1) = nkounts
+         lngth2(2) = 1
+         DO k=1, nbergs_in_file
+
+            strt(1) = k
+            strt2(2) = k
+
+            iret = NF90_GET_VAR(ncid, numberid, idata2, strt2, lngth2 )
+            localberg%number(:) = idata2(:)
+
+            iret = NF90_GET_VAR(ncid, scaling_id, data, strt, lngth )
+            localberg%mass_scaling = data(1)
+
+            iret = NF90_GET_VAR(ncid, lonid, data, strt, lngth)
+            localpt%lon = data(1)
+            iret = NF90_GET_VAR(ncid, latid, data, strt, lngth)
+            localpt%lat = data(1)
+            IF (nn_verbose_level >= 2 .AND. lwp) THEN
+               WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ',k,' is at ', &
+                                              localpt%lon,localpt%lat,' on PE ',narea-1
+            ENDIF
+            iret = NF90_GET_VAR(ncid, xid, data, strt, lngth)
+            localpt%xi = data(1)
+            iret = NF90_GET_VAR(ncid, yid, data, strt, lngth)
+            localpt%yj = data(1)
+            iret = NF90_GET_VAR(ncid, uvelid, data, strt, lngth )
+            localpt%uvel = data(1)
+            iret = NF90_GET_VAR(ncid, vvelid, data, strt, lngth )
+            localpt%vvel = data(1)
+            iret = NF90_GET_VAR(ncid, massid, data, strt, lngth )
+            localpt%mass = data(1)
+            iret = NF90_GET_VAR(ncid, thicknessid, data, strt, lngth )
+            localpt%thickness = data(1)
+            iret = NF90_GET_VAR(ncid, widthid, data, strt, lngth )
+            localpt%width = data(1)
+            iret = NF90_GET_VAR(ncid, lengthid, data, strt, lngth )
+            localpt%length = data(1)
+            iret = NF90_GET_VAR(ncid, yearid, idata, strt, lngth )
+            localpt%year = idata(1)
+            iret = NF90_GET_VAR(ncid, dayid, data, strt, lngth )
+            localpt%day = data(1)
+            iret = NF90_GET_VAR(ncid, mass_of_bits_id, data, strt, lngth )
+            localpt%mass_of_bits = data(1)
+            iret = NF90_GET_VAR(ncid, heat_density_id, data, strt, lngth )
+            localpt%heat_density = data(1)
+            !
+            CALL add_new_berg_to_list( localberg, localpt )
+         ENDDO
+
+      ENDIF
+
+      iret = NF90_INQ_DIMID( ncid, 'c', c_dim )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed')
+
+      iret = NF90_INQUIRE_DIMENSION( ncid, c_dim, dname, mclass )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed')
+
+      iret = NF90_INQ_VARID(ncid, 'kount'       , kountid)
+      iret = NF90_INQ_VARID(ncid, 'calving'     , calvid)
+      iret = NF90_INQ_VARID(ncid, 'calving_hflx', calvhid)
+      iret = NF90_INQ_VARID(ncid, 'stored_ice'  , siceid)
+      iret = NF90_INQ_VARID(ncid, 'stored_heat' , sheatid)
+
+      strt3(1) = 1
+      strt3(2) = 1
+      lngth3(1) = jpi
+      lngth3(2) = jpj
+      lngth3(3) = 1
+
+      DO k=1,mclass
+
+         strt3(3) = k
+
+         iret = NF90_GET_VAR( ncid, siceid , griddata, strt3, lngth3 )
+         berg_grid%stored_ice(:,:,k) = griddata(:,:,1)
+
+      ENDDO
+
+      iret = NF90_GET_VAR( ncid, calvid , p_calving(:,:) )
+      iret = NF90_GET_VAR( ncid, calvhid, p_calving_hflx(:,:) )
+      iret = NF90_GET_VAR( ncid, sheatid, berg_grid%stored_heat(:,:) )
+      iret = NF90_GET_VAR( ncid, kountid, idata2(:) )
+      kount_bergs(:) = idata2(:)
+
+      ! Finish up
+      iret = NF90_CLOSE(ncid)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed')
+
+      ! Sanity check
+      k = count_bergs()
+      IF (nn_verbose_level >= 0)   &
+         WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',k,' on PE',narea-1
+      IF( lk_mpp ) THEN
+         IF (multiPErestart) CALL mpp_sum(nbergs_in_file) ! In case PE 0 didn't open a file
+         CALL mpp_sum(k)
+      ENDIF
+      IF ( lwp ) THEN
+         WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, read_restart_bergs: there were',nbergs_in_file, &
+                                       ' bergs in the restart file and', k,' bergs have been read'
+      ENDIF
+
+      IF ( lwp .and. nn_verbose_level >= 0)  &
+         WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'
+
+   END SUBROUTINE icebergs_read_restart
+
+! ##############################################################################
+
+   SUBROUTINE icebergs_write_restart( kt )
+
+      INTEGER, INTENT( in )                 :: kt
+      ! Local variables
+      INTEGER                               :: iret, ncid, k_dim, n_dim, x_dim, y_dim, c_dim, k
+      CHARACTER(len=80)                     :: filename
+      TYPE(iceberg), POINTER                :: this
+      TYPE(point)  , POINTER                :: pt
+
+      IF( .NOT. ASSOCIATED(griddata) ) ALLOCATE( griddata(jpi,jpj,1) )
+
+      WRITE(filename,'("icebergs_",I8.8,"_restart_",I4.4,".nc")') kt, narea-1
+      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(filename)
+
+      iret = NF90_CREATE(TRIM(filename), NF90_CLOBBER, ncid)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed')
+
+      ! Dimensions
+      iret = NF90_DEF_DIM(ncid, 'x', jpi, x_dim)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim x failed')
+
+      iret = NF90_DEF_DIM(ncid, 'y', jpj, y_dim)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim y failed')
+
+      iret = NF90_DEF_DIM(ncid, 'c', nclasses, c_dim)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim c failed')
+
+      iret = NF90_DEF_DIM(ncid, 'k', nkounts, k_dim)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim k failed')
+
+      IF (associated(first_berg)) then
+         iret = NF90_DEF_DIM(ncid, 'n', NF90_UNLIMITED, n_dim)
+         IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_def_dim n failed')
+      ENDIF
+
+      ! Variables
+      iret = NF90_DEF_VAR(ncid, 'kount'       , NF90_INT   , (/ k_dim /), kountid)
+      iret = NF90_DEF_VAR(ncid, 'calving'     , NF90_DOUBLE, (/ x_dim, y_dim /), calvid)
+      iret = NF90_DEF_VAR(ncid, 'calving_hflx', NF90_DOUBLE, (/ x_dim, y_dim /), calvhid)
+      iret = NF90_DEF_VAR(ncid, 'stored_ice'  , NF90_DOUBLE, (/ x_dim, y_dim, c_dim /), siceid)
+      iret = NF90_DEF_VAR(ncid, 'stored_heat' , NF90_DOUBLE, (/ x_dim, y_dim /), sheatid)
+
+      ! Attributes
+      iret = NF90_PUT_ATT(ncid, calvid , 'long_name', 'iceberg calving')
+      iret = NF90_PUT_ATT(ncid, calvid , 'units', 'some')
+      iret = NF90_PUT_ATT(ncid, calvhid, 'long_name', 'heat flux associated with iceberg calving')
+      iret = NF90_PUT_ATT(ncid, calvhid, 'units', 'some')
+      iret = NF90_PUT_ATT(ncid, siceid , 'long_name', 'stored ice used to calve icebergs')
+      iret = NF90_PUT_ATT(ncid, siceid , 'units', 'kg/s')
+      iret = NF90_PUT_ATT(ncid, sheatid, 'long_name', 'heat in stored ice used to calve icebergs')
+      iret = NF90_PUT_ATT(ncid, sheatid, 'units', 'J/kg/s')
+
+      IF ( ASSOCIATED(first_berg) ) THEN
+
+         ! Only add berg variables for this PE if we have anything to say
+
+         ! Variables
+         iret = NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, n_dim, lonid)
+         iret = NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, n_dim, latid)
+         iret = NF90_DEF_VAR(ncid, 'xi', NF90_DOUBLE, n_dim, xid)
+         iret = NF90_DEF_VAR(ncid, 'yj', NF90_DOUBLE, n_dim, yid)
+         iret = NF90_DEF_VAR(ncid, 'uvel', NF90_DOUBLE, n_dim, uvelid)
+         iret = NF90_DEF_VAR(ncid, 'vvel', NF90_DOUBLE, n_dim, vvelid)
+         iret = NF90_DEF_VAR(ncid, 'mass', NF90_DOUBLE, n_dim, massid)
+         iret = NF90_DEF_VAR(ncid, 'thickness', NF90_DOUBLE, n_dim, thicknessid)
+         iret = NF90_DEF_VAR(ncid, 'width', NF90_DOUBLE, n_dim, widthid)
+         iret = NF90_DEF_VAR(ncid, 'length', NF90_DOUBLE, n_dim, lengthid)
+         iret = NF90_DEF_VAR(ncid, 'number', NF90_INT, (/k_dim,n_dim/), numberid)
+         iret = NF90_DEF_VAR(ncid, 'year', NF90_INT, n_dim, yearid)
+         iret = NF90_DEF_VAR(ncid, 'day', NF90_DOUBLE, n_dim, dayid)
+         iret = NF90_DEF_VAR(ncid, 'mass_scaling', NF90_DOUBLE, n_dim, scaling_id)
+         iret = NF90_DEF_VAR(ncid, 'mass_of_bits', NF90_DOUBLE, n_dim, mass_of_bits_id)
+         iret = NF90_DEF_VAR(ncid, 'heat_density', NF90_DOUBLE, n_dim, heat_density_id)
+
+         ! Attributes
+         iret = NF90_PUT_ATT(ncid, lonid, 'long_name', 'longitude')
+         iret = NF90_PUT_ATT(ncid, lonid, 'units', 'degrees_E')
+         iret = NF90_PUT_ATT(ncid, latid, 'long_name', 'latitude')
+         iret = NF90_PUT_ATT(ncid, latid, 'units', 'degrees_N')
+         iret = NF90_PUT_ATT(ncid, xid, 'long_name', 'x grid box position')
+         iret = NF90_PUT_ATT(ncid, xid, 'units', 'fractional')
+         iret = NF90_PUT_ATT(ncid, yid, 'long_name', 'y grid box position')
+         iret = NF90_PUT_ATT(ncid, yid, 'units', 'fractional')
+         iret = NF90_PUT_ATT(ncid, uvelid, 'long_name', 'zonal velocity')
+         iret = NF90_PUT_ATT(ncid, uvelid, 'units', 'm/s')
+         iret = NF90_PUT_ATT(ncid, vvelid, 'long_name', 'meridional velocity')
+         iret = NF90_PUT_ATT(ncid, vvelid, 'units', 'm/s')
+         iret = NF90_PUT_ATT(ncid, massid, 'long_name', 'mass')
+         iret = NF90_PUT_ATT(ncid, massid, 'units', 'kg')
+         iret = NF90_PUT_ATT(ncid, thicknessid, 'long_name', 'thickness')
+         iret = NF90_PUT_ATT(ncid, thicknessid, 'units', 'm')
+         iret = NF90_PUT_ATT(ncid, widthid, 'long_name', 'width')
+         iret = NF90_PUT_ATT(ncid, widthid, 'units', 'm')
+         iret = NF90_PUT_ATT(ncid, lengthid, 'long_name', 'length')
+         iret = NF90_PUT_ATT(ncid, lengthid, 'units', 'm')
+         iret = NF90_PUT_ATT(ncid, numberid, 'long_name', 'iceberg number on this processor')
+         iret = NF90_PUT_ATT(ncid, numberid, 'units', 'count')
+         iret = NF90_PUT_ATT(ncid, yearid, 'long_name', 'calendar year of calving event')
+         iret = NF90_PUT_ATT(ncid, yearid, 'units', 'years')
+         iret = NF90_PUT_ATT(ncid, dayid, 'long_name', 'year day of calving event')
+         iret = NF90_PUT_ATT(ncid, dayid, 'units', 'days')
+         iret = NF90_PUT_ATT(ncid, scaling_id, 'long_name', 'scaling factor for mass of calving berg')
+         iret = NF90_PUT_ATT(ncid, scaling_id, 'units', 'none')
+         iret = NF90_PUT_ATT(ncid, mass_of_bits_id, 'long_name', 'mass of bergy bits')
+         iret = NF90_PUT_ATT(ncid, mass_of_bits_id, 'units', 'kg')
+         iret = NF90_PUT_ATT(ncid, heat_density_id, 'long_name', 'heat density')
+         iret = NF90_PUT_ATT(ncid, heat_density_id, 'units', 'J/kg')
+
+      ENDIF ! associated(first_berg)
+
+      ! End define mode
+      iret = NF90_ENDDEF(ncid)
+
+      ! --------------------------------
+      ! now write some data
+
+      strt3(1) = 1
+      strt3(2) = 1
+      lngth3(1) = jpi
+      lngth3(2) = jpj
+      lngth3(3) = 1
+
+      DO k=1,nclasses
+         griddata(:,:,1) = berg_grid%stored_ice(:,:,k)
+         strt3(3) = k
+         iret = NF90_PUT_VAR( ncid, siceid, griddata, strt3, lngth3 )
+         IF (iret .ne. NF90_NOERR) THEN
+            IF( lwp ) WRITE(numout,*) TRIM(NF90_STRERROR( iret ))
+            CALL ctl_stop('icebergs, write_restart: nf_put_var stored_ice failed')
+         ENDIF
+      ENDDO
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(filename),' var: stored_ice  written'
+
+      iret = NF90_PUT_VAR( ncid, kountid, kount_bergs(:) )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var kount failed')
+
+      iret = NF90_PUT_VAR( ncid, sheatid, berg_grid%stored_heat(:,:) )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed')
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(filename),' var: stored_heat written'
+
+      iret = NF90_PUT_VAR( ncid, calvid , p_calving(:,:) )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving failed')
+      iret = NF90_PUT_VAR( ncid, calvhid, p_calving_hflx(:,:) )
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed')
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(filename),' var: calving written'
+
+      IF ( ASSOCIATED(first_berg) ) THEN
+
+         ! Write variables
+         ! sga - just write out the current point of the trajectory
+
+         this => first_berg
+         k = 0
+         DO WHILE (ASSOCIATED(this))
+            pt => this%current_point
+            k=k+1
+
+            iret = NF90_PUT_VAR(ncid, numberid, this%number, (/1,k/), (/nkounts,1/) )
+            iret = NF90_PUT_VAR(ncid, scaling_id, this%mass_scaling, (/ k /) )
+
+            iret = NF90_PUT_VAR(ncid, lonid, pt%lon, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, latid, pt%lat, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, xid, pt%xi, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, yid, pt%yj, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, uvelid, pt%uvel, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, vvelid, pt%vvel, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, massid, pt%mass, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, thicknessid, pt%thickness, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, widthid, pt%width, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, lengthid, pt%length, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, yearid, pt%year, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, dayid, pt%day, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, mass_of_bits_id, pt%mass_of_bits, (/ k /) )
+            iret = NF90_PUT_VAR(ncid, heat_density_id, pt%heat_density, (/ k /) )
+
+            this=>this%next
+         ENDDO
+
+      ENDIF ! associated(first_berg)
+
+      ! Finish up
+      iret = NF90_CLOSE(ncid)
+      IF (iret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_close failed')
+
+   END SUBROUTINE icebergs_write_restart
+
+END MODULE icbrst
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrun.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrun.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbrun.F90	(revision 3339)
@@ -0,0 +1,219 @@
+MODULE icbrun
+
+   !!======================================================================
+   !!                       ***  MODULE  icbrun  ***
+   !! Ocean physics:  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 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 xxx.F90 module
+   PUBLIC   icb_end        ! routine called in xxx.F90 module
+
+CONTAINS
+
+   SUBROUTINE icb_stp( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_stp  ***
+      !!
+      !! ** Purpose :   iceberg time stepping.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )           ::  kt
+      !
+      INTEGER                         ::   iyr, imon, iday, ihr, imin, isec
+      LOGICAL                         ::   lerr, sample_traj, l_budget, l_verbose
+      !!----------------------------------------------------------------------
+      !! start of timestep housekeeping
+
+      ktberg = kt
+
+      ! read calving data
+      IF( nn_test_icebergs < 0 ) THEN
+         !
+         CALL fld_read ( kt, 1, sf_icb )
+         p_calving(:,:)      = sf_icb(1)%fnow(:,:,1)    ! calving in km^3/year (water equivalent)
+         p_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_budget_step()
+
+      ! Manage time
+      ! Convert nemo time variables from dom_oce into local versions
+      ! Note that yearday function assumes 365 day year!!
+      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
+      current_year    = iyr
+      current_yearday = yearday(imon, iday, ihr, imin, isec)
+
+      l_verbose = .FALSE.
+      IF( nn_verbose_write .GT. 0 .AND. &
+          MOD(kt-1,nn_verbose_write ) == 0 ) l_verbose = nn_verbose_level >= 0
+      IF( l_verbose ) WRITE(numicb,9100) iyr, imon, iday, ihr, imin, isec, &
+                                        current_year, current_yearday
+ 9100 FORMAT('y,m,d=',3i5,' h,m,s=',3i5,' yr,yrdy=',i5,f8.3)
+
+      ! copy nemo forcing arrays into iceberg versions with extra halo
+      ! only necessary for variables not on T points
+      CALL copy_flds()
+
+      !!----------------------------------------------------------------------
+      !! process icebergs
+
+      ! Accumulate ice from calving
+      CALL accumulate_calving( kt )
+
+      ! Calve excess stored ice into icebergs
+      CALL calve_icebergs()
+
+      !                               !==  For each berg, evolve  ==!
+      !
+      IF( ASSOCIATED(first_berg) )   CALL evolve_icebergs()     ! ice berg dynamics
+
+      IF( lk_mpp ) THEN
+                                      CALL mpp_send_bergs ()     ! Send bergs to other PEs
+      ELSE
+                                      CALL lbc_send_bergs()      ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+
+      IF( ASSOCIATED(first_berg) )   CALL thermodynamics ( kt ) ! Ice berg thermodynamics (melting) + rolling
+
+      !!----------------------------------------------------------------------
+      !! end of timestep housekeeping
+
+      sample_traj = .FALSE.
+      IF( nn_sample_rate .GT. 0 .AND. MOD(kt-1,nn_sample_rate) == 0 )   sample_traj = .TRUE.
+      IF( sample_traj .AND.   &
+          ASSOCIATED(first_berg) )   CALL traj_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'
+
+      ! write out mean budgets - not sure why this happens before they are calculated - sga !
+      CALL icb_budget_put()
+
+      ! Dump icebergs to screen
+      if ( nn_verbose_level >= 2 ) CALL print_bergs( 'icb_stp, status', kt )
+
+      ! Diagnose budgets
+      l_budget = .FALSE.
+      IF( nn_verbose_write .GT. 0 .AND. MOD(kt-1,nn_verbose_write) == 0 )   l_budget = ln_bergdia
+      CALL icb_budget( l_budget )
+
+      IF( MOD(kt,nn_stock) == 0 ) THEN
+         CALL icebergs_write_restart( kt )
+         IF( nn_sample_rate .GT. 0 ) CALL traj_sync()
+      ENDIF
+
+      !
+   END SUBROUTINE icb_stp
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE icb_end( kt )
+
+      ! Arguments
+      INTEGER, INTENT( in )  :: kt
+      ! Local variables
+      TYPE(iceberg), POINTER :: this, next
+
+      ! expanded arrays for bilinear interpolation
+      DEALLOCATE( uo_e )
+      DEALLOCATE( vo_e )
+      DEALLOCATE( ff_e )
+      DEALLOCATE( ua_e )
+      DEALLOCATE( va_e )
+#if defined key_lim2 || defined key_lim3
+      DEALLOCATE( ui_e )
+      DEALLOCATE( vi_e )
+#endif
+      DEALLOCATE( ssh_e )
+
+      DEALLOCATE( icbfldpts )
+
+      IF( lk_mpp ) CALL dealloc_buffers()
+
+      IF (.NOT.ASSOCIATED(berg_grid)) RETURN
+
+      ! only write a restart if not done in icb_stp
+      IF( MOD(kt,nn_stock) .NE. 0 ) CALL icebergs_write_restart( kt )
+
+      ! finish with trajectories if they were written
+      IF( nn_sample_rate .GT. 0 ) CALL traj_end()
+
+      ! Delete bergs and structures
+      this=>first_berg
+      DO WHILE (ASSOCIATED(this))
+        next=>this%next
+        CALL destroy_iceberg(this)
+        this=>next
+      ENDDO
+
+      CALL icb_budget_end()
+
+      DEALLOCATE(berg_grid%calving)
+      DEALLOCATE(berg_grid%calving_hflx)
+      DEALLOCATE(berg_grid%stored_heat)
+      DEALLOCATE(berg_grid%floating_melt)
+      DEALLOCATE(berg_grid%maxclass)
+      DEALLOCATE(berg_grid%tmp)
+      DEALLOCATE(berg_grid%stored_ice)
+      DEALLOCATE(berg_grid)
+
+      DEALLOCATE(initial_width)
+      DEALLOCATE(initial_length)
+
+      IF (lwp) WRITE(numout,'(a,i6)') 'icebergs: icb_end complete',narea
+      CALL flush( numicb )
+      CLOSE( numicb )
+
+   END SUBROUTINE icb_end
+
+   !!-------------------------------------------------------------------------
+
+END MODULE icbrun
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 3339)
@@ -0,0 +1,209 @@
+MODULE icbthm
+
+   !!======================================================================
+   !!                       ***  MODULE  icbthm  ***
+   !! Ocean physics:  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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   thermodynamics : 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   thermodynamics ! routine called in xxx.F90 module
+
+CONTAINS
+
+   SUBROUTINE thermodynamics( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE thermodynamics  ***
+      !!
+      !! ** Purpose :   compute the iceberg thermodynamics.
+      !!
+      !! ** Method  : - blah blah
+      !!----------------------------------------------------------------------
+      INTEGER                         ::   kt          ! timestep number, just passed to print_berg
+      !
+      REAL(wp)                        ::   M, T, W, L, SST, Vol, Ln, Wn, Tn, nVol, IC, Dn
+      REAL(wp)                        ::   Mv, Me, Mb, melt, dvo, dva, dM, Ss, dMe, dMb, dMv
+      REAL(wp)                        ::   Mnew, Mnew1, Mnew2, heat
+      REAL(wp)                        ::   Mbits, nMbits, dMbitsE, dMbitsM, Lbits, Abits, Mbb
+      REAL(wp)                        ::   xi, yj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
+      INTEGER                         ::   ii, ij
+      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
+         knberg = this%number(1)
+         CALL interp_flds( 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 )
+         !
+         SST = pt%sst
+         IC  = MIN( 1._wp, pt%cn + rn_sicn_shift )     ! Shift sea-ice concentration       !!gm ???
+         M   = pt%mass
+         T   = pt%thickness                               ! total thickness
+       ! D   = (rn_rho_bergs/rho_seawater)*T ! draught (keel depth)
+       ! F   = T - D ! freeboard
+         W   = pt%width
+         L   = pt%length
+         xi  = pt%xi                                      ! position in (i,j) referential
+         yj  = pt%yj
+         ii  = INT( xi + 0.5 ) - nimpp + 1                    ! t-cell of the berg
+         ij  = INT( yj + 0.5 ) - njmpp + 1
+         Vol = T * W * L
+         zdt = berg_dt   ;   z1_dt = 1._wp / zdt
+
+         ! Environment
+         dvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
+         dva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
+         Ss  = 1.5 * SQRT( dva ) + 0.1 * dva                ! Sea state      (eqn M.A9)
+
+         ! Melt rates in m/s (i.e. division by rday)
+         Mv = MAX( 7.62e-3*SST+1.29e-3*(SST**2)            , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10)
+         Mb = MAX( 0.58*(dvo**0.8)*(SST+4.0)/(L**0.2)      , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 )
+         Me = MAX( 1./12.*(SST+2.)*Ss*(1+cos(rpi*(IC**3))) , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
+
+         IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
+            Tn    = MAX( T - Mb*zdt , 0._wp )         ! new total thickness (m)
+            nVol  = Tn * W * L                        ! new volume (m^3)
+            Mnew1 = (nVol/Vol) * M                    ! new mass (kg)
+            dMb   = M - Mnew1                         ! mass lost to basal melting (>0) (kg)
+            !
+            Ln    = MAX( L - Mv*zdt , 0._wp )         ! new length (m)
+            Wn    = MAX( W - Mv*zdt , 0._wp )         ! new width (m)
+            nVol  = Tn * Wn * Ln                      ! new volume (m^3)
+            Mnew2 = (nVol/Vol) * M                    ! new mass (kg)
+            dMv   = Mnew1 - Mnew2                     ! mass lost to buoyant convection (>0) (kg)
+            !
+            Ln    = MAX( Ln - Me*zdt , 0._wp )        ! new length (m)
+            Wn    = MAX( Wn - Me*zdt , 0._wp )        ! new width (m)
+            nVol  = Tn * Wn * Ln                      ! new volume (m^3)
+            Mnew  = ( nVol / Vol ) * M                ! new mass (kg)
+            dMe   = Mnew2 - Mnew                      ! mass lost to erosion (>0) (kg)
+            dM    = M - Mnew                          ! mass lost to all erosion and melting (>0) (kg)
+            !
+         ELSE                                         ! Update dimensions of berg
+            Ln = MAX( L -(Mv+Me)*zdt ,0._wp )         ! (m)
+            Wn = MAX( W -(Mv+Me)*zdt ,0._wp )         ! (m)
+            Tn = MAX( T - Mb    *zdt ,0._wp )         ! (m)
+            ! Update volume and mass of berg
+            nVol = Tn*Wn*Ln                           ! (m^3)
+            Mnew = (nVol/Vol)*M                       ! (kg)
+            dM   = M - Mnew                           ! (kg)
+            dMb = (M/Vol) * (W*   L ) *Mb*zdt         ! approx. mass loss to basal melting (kg)
+            dMe = (M/Vol) * (T*(W+L)) *Me*zdt         ! approx. mass lost to erosion (kg)
+            dMv = (M/Vol) * (T*(W+L)) *Mv*zdt         ! approx. mass loss to buoyant convection (kg)
+         ENDIF
+
+         IF( rn_bits_erosion_fraction > 0._wp ) THEN      ! Bergy bits
+            !
+            Mbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
+            dMbitsE = rn_bits_erosion_fraction * dMe                        ! change in mass of bits (kg)
+            nMbits  = Mbits + dMbitsE                                               ! add new bergy bits to mass (kg)
+            Lbits   = MIN( L, W, T, 40._wp )                                        ! assume bergy bits are smallest dimension or 40 meters
+            Abits   = ( Mbits / rn_rho_bergs ) / Lbits                           ! Effective bottom area (assuming T=Lbits)
+            Mbb     = MAX( 0.58*(dvo**0.8)*(SST+2.0)/(Lbits**0.2), 0.) * z1_rday    ! Basal turbulent melting (for bits)
+            Mbb     = rn_rho_bergs * Abits * Mbb                                 ! in kg/s
+            dMbitsM = MIN( Mbb*zdt , nMbits )                                       ! bergy bits mass lost to melting (kg)
+            nMbits  = nMbits-dMbitsM                                                ! remove mass lost to bergy bits melt
+            IF( Mnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
+               dMbitsM = dMbitsM + nMbits                                           ! instantly melt all the bergy bits
+               nMbits  = 0._wp
+            ENDIF
+         ELSE                                                     ! No bergy bits
+            Abits   = 0._wp
+            dMbitsE = 0._wp
+            dMbitsM = 0._wp
+            nMbits  = 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
+            melt    = ( dM - ( dMbitsE - dMbitsM ) ) * z1_dt   ! kg/s
+            berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + melt    * z1_e1e2    ! kg/m2/s
+            heat = melt * pt%heat_density              ! kg/s x J/kg = J/s
+            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + heat    * z1_e1e2    ! W/m2
+            CALL melt_budget(ii, ij, Mnew, heat, this%mass_scaling, dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2 )
+         ELSE
+            WRITE(numout,*) 'thermodynamics: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
+            CALL print_berg( this, kt )
+            WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
+            CALL ctl_stop('thermodynamics', 'berg appears to have grounded!')
+         ENDIF
+
+         ! Rolling
+         Dn = ( rn_rho_bergs / rho_seawater ) * Tn       ! draught (keel depth)
+         IF( Dn > 0._wp .AND. MAX(Wn,Ln) < SQRT( 0.92*(Dn**2) + 58.32*Dn ) ) THEN
+            T  = Tn
+            Tn = Wn
+            Wn = T
+         endif
+
+         ! Store the new state of iceberg (with L>W)
+         pt%mass         = Mnew
+         pt%mass_of_bits = nMbits
+         pt%thickness    = Tn
+         pt%width        = min(Wn,Ln)
+         pt%length       = max(Wn,Ln)
+
+         next=>this%next
+
+!!gm  add a test to avoid over melting ?
+
+         IF( Mnew <= 0._wp ) THEN       ! Delete the berg if completely melted
+            CALL delete_iceberg_from_list( first_berg, this )
+            !
+         ELSE                            ! Diagnose mass distribution on grid
+            z1_e1e2                 = 1._wp / e1e2t(ii,ij) * this%mass_scaling
+            CALL size_budget(ii, ij, Wn, Ln, Abits, this%mass_scaling, Mnew, nMbits, 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(:,:)  !! heat flux not yet properly coded
+      ENDIF
+      !
+   END  SUBROUTINE thermodynamics
+
+END MODULE icbthm
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 3339)
@@ -0,0 +1,252 @@
+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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   write_trajectory      :
+   !!----------------------------------------------------------------------
+   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
+
+   INTEGER, PRIVATE   ::   kount_traj
+   INTEGER, PRIVATE   ::   n_dim, k_dim
+   INTEGER, PRIVATE   ::   trajid
+   INTEGER, PRIVATE   ::   numberid, stepid, scaling_id
+   INTEGER, PRIVATE   ::   lonid, latid, xid, yid, uvelid, vvelid, massid
+   INTEGER, PRIVATE   ::   uoid, void, uaid, vaid, uiid, viid
+   INTEGER, PRIVATE   ::   sshxid, sshyid, sstid, cntid, thkid
+   INTEGER, PRIVATE   ::   thicknessid, widthid, lengthid
+   INTEGER, PRIVATE   ::   yearid, dayid
+   INTEGER, PRIVATE   ::   mass_of_bits_id, heat_density_id
+
+   PUBLIC   traj_init           ! routine called in xxx.F90 module
+   PUBLIC   traj_write          ! routine called in xxx.F90 module
+   PUBLIC   traj_sync           ! routine called in xxx.F90 module
+   PUBLIC   traj_end            ! routine called in xxx.F90 module
+
+CONTAINS
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE traj_init( nitend )
+
+      ! local variables
+      INTEGER, INTENT( in )                 :: nitend
+      ! Local variables
+      INTEGER                               :: iret, k
+      CHARACTER(len=80)                     :: filename
+      TYPE(iceberg), POINTER                :: this
+      TYPE(point)  , POINTER                :: pt
+
+      IF( lk_mpp ) THEN
+         WRITE(filename,'("trajectory_icebergs_",I6.6,"_",I4.4,".nc")') nitend, narea-1
+      ELSE
+         WRITE(filename,'("trajectory_icebergs_",I6.6,".nc")') nitend
+      ENDIF
+      IF ( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, traj_init: creating ',TRIM(filename)
+
+      iret = NF90_CREATE(TRIM(filename), NF90_CLOBBER, trajid)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_create failed')
+
+      ! Dimensions
+      iret = NF90_DEF_DIM(trajid, 'n', NF90_UNLIMITED, n_dim)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_def_dim n failed')
+      iret = NF90_DEF_DIM(trajid, 'k', nkounts, k_dim)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_init: nf_def_dim k failed')
+
+      ! Variables
+      iret = NF90_DEF_VAR(trajid, 'iceberg_number', NF90_INT, (/k_dim,n_dim/), numberid)
+      iret = NF90_DEF_VAR(trajid, 'timestep', NF90_INT, n_dim, stepid)
+      iret = NF90_DEF_VAR(trajid, 'mass_scaling', NF90_DOUBLE, n_dim, scaling_id)
+      iret = NF90_DEF_VAR(trajid, 'lon', NF90_DOUBLE, n_dim, lonid)
+      iret = NF90_DEF_VAR(trajid, 'lat', NF90_DOUBLE, n_dim, latid)
+      iret = NF90_DEF_VAR(trajid, 'xi', NF90_DOUBLE, n_dim, xid)
+      iret = NF90_DEF_VAR(trajid, 'yj', NF90_DOUBLE, n_dim, yid)
+      iret = NF90_DEF_VAR(trajid, 'uvel', NF90_DOUBLE, n_dim, uvelid)
+      iret = NF90_DEF_VAR(trajid, 'vvel', NF90_DOUBLE, n_dim, vvelid)
+      iret = NF90_DEF_VAR(trajid, 'uto', NF90_DOUBLE, n_dim, uoid)
+      iret = NF90_DEF_VAR(trajid, 'vto', NF90_DOUBLE, n_dim, void)
+      iret = NF90_DEF_VAR(trajid, 'uta', NF90_DOUBLE, n_dim, uaid)
+      iret = NF90_DEF_VAR(trajid, 'vta', NF90_DOUBLE, n_dim, vaid)
+      iret = NF90_DEF_VAR(trajid, 'uti', NF90_DOUBLE, n_dim, uiid)
+      iret = NF90_DEF_VAR(trajid, 'vti', NF90_DOUBLE, n_dim, viid)
+      iret = NF90_DEF_VAR(trajid, 'ssh_x', NF90_DOUBLE, n_dim, sshyid)
+      iret = NF90_DEF_VAR(trajid, 'ssh_y', NF90_DOUBLE, n_dim, sshxid)
+      iret = NF90_DEF_VAR(trajid, 'sst', NF90_DOUBLE, n_dim, sstid)
+      iret = NF90_DEF_VAR(trajid, 'icnt', NF90_DOUBLE, n_dim, cntid)
+      iret = NF90_DEF_VAR(trajid, 'ithk', NF90_DOUBLE, n_dim, thkid)
+      iret = NF90_DEF_VAR(trajid, 'mass', NF90_DOUBLE, n_dim, massid)
+      iret = NF90_DEF_VAR(trajid, 'thickness', NF90_DOUBLE, n_dim, thicknessid)
+      iret = NF90_DEF_VAR(trajid, 'width', NF90_DOUBLE, n_dim, widthid)
+      iret = NF90_DEF_VAR(trajid, 'length', NF90_DOUBLE, n_dim, lengthid)
+      iret = NF90_DEF_VAR(trajid, 'year', NF90_INT, n_dim, yearid)
+      iret = NF90_DEF_VAR(trajid, 'day', NF90_DOUBLE, n_dim, dayid)
+      iret = NF90_DEF_VAR(trajid, 'mass_of_bits', NF90_DOUBLE, n_dim, mass_of_bits_id)
+      iret = NF90_DEF_VAR(trajid, 'heat_density', NF90_DOUBLE, n_dim, heat_density_id)
+
+      ! Attributes
+      iret = NF90_PUT_ATT(trajid, numberid, 'long_name', 'iceberg number on this processor')
+      iret = NF90_PUT_ATT(trajid, numberid, 'units', 'count')
+      iret = NF90_PUT_ATT(trajid, stepid, 'long_name', 'timestep number kt')
+      iret = NF90_PUT_ATT(trajid, stepid, 'units', 'count')
+      iret = NF90_PUT_ATT(trajid, lonid, 'long_name', 'longitude')
+      iret = NF90_PUT_ATT(trajid, lonid, 'units', 'degrees_E')
+      iret = NF90_PUT_ATT(trajid, latid, 'long_name', 'latitude')
+      iret = NF90_PUT_ATT(trajid, latid, 'units', 'degrees_N')
+      iret = NF90_PUT_ATT(trajid, xid, 'long_name', 'x grid box position')
+      iret = NF90_PUT_ATT(trajid, xid, 'units', 'fractional')
+      iret = NF90_PUT_ATT(trajid, yid, 'long_name', 'y grid box position')
+      iret = NF90_PUT_ATT(trajid, yid, 'units', 'fractional')
+      iret = NF90_PUT_ATT(trajid, uvelid, 'long_name', 'zonal velocity')
+      iret = NF90_PUT_ATT(trajid, uvelid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, vvelid, 'long_name', 'meridional velocity')
+      iret = NF90_PUT_ATT(trajid, vvelid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, uoid, 'long_name', 'ocean u component')
+      iret = NF90_PUT_ATT(trajid, uoid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, void, 'long_name', 'ocean v component')
+      iret = NF90_PUT_ATT(trajid, void, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, uaid, 'long_name', 'atmosphere u component')
+      iret = NF90_PUT_ATT(trajid, uaid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, vaid, 'long_name', 'atmosphere v component')
+      iret = NF90_PUT_ATT(trajid, vaid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, uiid, 'long_name', 'sea ice u component')
+      iret = NF90_PUT_ATT(trajid, uiid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, viid, 'long_name', 'sea ice v component')
+      iret = NF90_PUT_ATT(trajid, viid, 'units', 'm/s')
+      iret = NF90_PUT_ATT(trajid, sshxid, 'long_name', 'sea surface height gradient from x points')
+      iret = NF90_PUT_ATT(trajid, sshxid, 'units', 'm/m')
+      iret = NF90_PUT_ATT(trajid, sshyid, 'long_name', 'sea surface height gradient from y points')
+      iret = NF90_PUT_ATT(trajid, sshyid, 'units', 'm/m')
+      iret = NF90_PUT_ATT(trajid, sstid, 'long_name', 'sea surface temperature')
+      iret = NF90_PUT_ATT(trajid, sstid, 'units', 'degC')
+      iret = NF90_PUT_ATT(trajid, cntid, 'long_name', 'sea ice concentration')
+      iret = NF90_PUT_ATT(trajid, cntid, 'units', 'degC')
+      iret = NF90_PUT_ATT(trajid, thkid, 'long_name', 'sea ice thickness')
+      iret = NF90_PUT_ATT(trajid, thkid, 'units', 'm')
+      iret = NF90_PUT_ATT(trajid, massid, 'long_name', 'mass')
+      iret = NF90_PUT_ATT(trajid, massid, 'units', 'kg')
+      iret = NF90_PUT_ATT(trajid, thicknessid, 'long_name', 'thickness')
+      iret = NF90_PUT_ATT(trajid, thicknessid, 'units', 'm')
+      iret = NF90_PUT_ATT(trajid, widthid, 'long_name', 'width')
+      iret = NF90_PUT_ATT(trajid, widthid, 'units', 'm')
+      iret = NF90_PUT_ATT(trajid, lengthid, 'long_name', 'length')
+      iret = NF90_PUT_ATT(trajid, lengthid, 'units', 'm')
+      iret = NF90_PUT_ATT(trajid, yearid, 'long_name', 'calendar year')
+      iret = NF90_PUT_ATT(trajid, yearid, 'units', 'years')
+      iret = NF90_PUT_ATT(trajid, dayid, 'long_name', 'day of year')
+      iret = NF90_PUT_ATT(trajid, dayid, 'units', 'days')
+      iret = NF90_PUT_ATT(trajid, scaling_id, 'long_name', 'scaling factor for mass of berg')
+      iret = NF90_PUT_ATT(trajid, scaling_id, 'units', 'none')
+      iret = NF90_PUT_ATT(trajid, mass_of_bits_id, 'long_name', 'mass of bergy bits')
+      iret = NF90_PUT_ATT(trajid, mass_of_bits_id, 'units', 'kg')
+      iret = NF90_PUT_ATT(trajid, heat_density_id, 'long_name', 'heat density')
+      iret = NF90_PUT_ATT(trajid, heat_density_id, 'units', 'J/kg')
+
+      ! End define mode
+      iret = NF90_ENDDEF(trajid)
+
+   END SUBROUTINE traj_init
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE traj_write( kt )
+      !! 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
+      ! Local variables
+      INTEGER                               :: iret, k
+      CHARACTER(len=80)                     :: filename
+      TYPE(iceberg), POINTER                :: this
+      TYPE(point)  , POINTER                :: pt
+
+      ! Write variables
+      ! sga - just write out the current point of the trajectory
+
+      this => first_berg
+      k = kount_traj
+      DO WHILE (ASSOCIATED(this))
+         pt => this%current_point
+         k=k+1
+
+         iret = NF90_PUT_VAR(trajid, numberid, this%number, (/1,k/), (/nkounts,1/) )
+         iret = NF90_PUT_VAR(trajid, stepid, kt, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, scaling_id, this%mass_scaling, (/ k /) )
+
+         iret = NF90_PUT_VAR(trajid, lonid, pt%lon, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, latid, pt%lat, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, xid, pt%xi, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, yid, pt%yj, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, uvelid, pt%uvel, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, vvelid, pt%vvel, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, uoid, pt%uo, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, void, pt%vo, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, uaid, pt%ua, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, vaid, pt%va, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, uiid, pt%ui, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, viid, pt%vi, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, sshxid, pt%ssh_x, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, sshyid, pt%ssh_y, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, sstid, pt%sst, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, cntid, pt%cn, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, thkid, pt%hi, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, massid, pt%mass, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, thicknessid, pt%thickness, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, widthid, pt%width, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, lengthid, pt%length, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, yearid, pt%year, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, dayid, pt%day, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, mass_of_bits_id, pt%mass_of_bits, (/ k /) )
+         iret = NF90_PUT_VAR(trajid, heat_density_id, pt%heat_density, (/ k /) )
+
+         this=>this%next
+      ENDDO
+      IF( lwp .and. nn_verbose_level > 0 ) WRITE(numout,*) 'trajectory write to frame ', k
+      kount_traj = k
+
+   END SUBROUTINE traj_write
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE traj_sync()
+      ! Local variables
+      INTEGER                               :: iret
+
+      ! flush to file
+      iret = NF90_SYNC(trajid)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_sync: nf_sync failed')
+
+   END SUBROUTINE traj_sync
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE traj_end()
+      ! Local variables
+      INTEGER                               :: iret
+
+      ! Finish up
+      iret = NF90_CLOSE(trajid)
+      IF (iret .NE. NF90_NOERR) CALL ctl_stop('icebergs, traj_end: nf_close failed')
+
+   END SUBROUTINE traj_end
+
+   !!-------------------------------------------------------------------------
+
+END MODULE icbtrj
Index: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
===================================================================
--- branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 3339)
+++ branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 3339)
@@ -0,0 +1,667 @@
+MODULE icbutl
+
+   !!======================================================================
+   !!                       ***  MODULE  icbutl  ***
+   !! Ocean physics:  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
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   interp_flds   :
+   !!   bilin         :
+   !!   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: hi => 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   copy_flds                ! routine called in xxx.F90 module
+   PUBLIC   interp_flds              ! routine called in xxx.F90 module
+   PUBLIC   bilin
+   PUBLIC   bilin_x
+   PUBLIC   bilin_e
+   PUBLIC   add_new_berg_to_list
+   PUBLIC   insert_berg_into_list
+   PUBLIC   delete_iceberg_from_list
+   PUBLIC   destroy_iceberg
+   PUBLIC   track_berg
+   PUBLIC   print_berg
+   PUBLIC   print_bergs
+   PUBLIC   count_bergs
+   PUBLIC   increment_kounter
+   PUBLIC   yearday
+   PUBLIC   sum_mass
+   PUBLIC   sum_heat
+
+   PRIVATE  create_iceberg
+
+   !!-------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE copy_flds()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE copy_flds  ***
+      !!
+      !! ** 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 copy_flds
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE interp_flds( pi, pe1, puo, pui, pua, pssh_i,                     &
+      &                    pj, pe2, pvo, pvi, pva, pssh_j, psst, pcn, phi, pff )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interp_flds  ***
+      !!
+      !! ** Purpose :   iceberg initialization.
+      !!
+      !! ** Method  : - blah blah
+      !!
+      !!       !!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 = bilin_e( e1t, e1u, e1v, e1f, pi, pj )                 ! scale factors
+      pe2 = bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+      !
+      puo = bilin( uo_e, pi, pj, 'U', 1, 1 )                      ! ocean velocities
+      pvo = bilin( vo_e, pi, pj, 'V', 1, 1 )
+      psst = bilin( sst_m, pi, pj, 'T', 0, 0 )                    ! SST
+      pcn  = bilin( fr_i , pi, pj, 'T', 0, 0 )                    ! ice concentration
+      pff  = bilin( ff_e , pi, pj, 'F', 1, 1 )                    ! Coriolis parameter
+      !
+      pua = bilin( ua_e , pi, pj, 'U', 1, 1 )                     ! 10m wind
+      pva = bilin( va_e , pi, pj, 'V', 1, 1 )                     ! 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 = bilin( ui_e, pi, pj, 'U', 1, 1 )                      ! sea-ice velocities
+      pvi = bilin( vi_e, pi, pj, 'V', 1, 1 )
+      phi = bilin( hi  , pi, pj, 'T', 0, 0 )                      ! 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 = ( bilin( ssh_e, pi+0.1_wp, pj, 'T', 1, 1 ) - bilin( ssh_e, pi-0.1_wp, pj, 'T', 1, 1 )  ) / ( 0.2_wp * pe1 )
+      pssh_j = ( bilin( ssh_e, pi, pj+0.1_wp, 'T', 1, 1 ) - bilin( ssh_e, pi, pj-0.1_wp, 'T', 1, 1 )  ) / ( 0.2_wp * pe2 )
+      !
+   END SUBROUTINE interp_flds
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION bilin( pfld, pi, pj, cd_type, jdi, jdj )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION 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
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                         , INTENT(in) ::   jdi, jdj  ! extra halo on grid
+      REAL(wp), DIMENSION(1-jdi:jpi+jdi,1-jdj:jpj+jdj), 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 = ii - nimpp + 1
+      ij = ij - njmpp + 1
+      !
+      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 bilin
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION bilin_x( pfld, pi, pj )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION 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), 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
+      ii = ii - nimpp + 1
+      ij = ij - njmpp + 1
+      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
+      !
+      bilin_x = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
+      IF( bilin_x > 180._wp ) bilin_x = bilin_x - 360._wp
+      !
+   END FUNCTION bilin_x
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION 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 point 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)
+      zj = pj - REAL(ij,wp)
+
+      ! find position in this processor
+      ii = ii - nimpp + 1
+      ij = ij - njmpp + 1
+
+      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
+      !
+      bilin_e = ( ze01 * (1.-zi) + ze11 * zi ) *     zj    &
+         &    + ( ze00 * (1.-zi) + ze10 * zi ) * (1.-zj)
+      !
+   END FUNCTION bilin_e
+
+   !!-------------------------------------------------------------------------
+
+
+   SUBROUTINE add_new_berg_to_list( bergvals, ptvals )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE add_new_berg_to_list  ***
+      !!
+      !! ** Purpose :   add a new berg to the iceberg list
+      !!
+      !! ** method  : - ?
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), INTENT(in)           ::   bergvals
+      TYPE(point)  , INTENT(in)           ::   ptvals
+      !
+      TYPE(iceberg), POINTER ::   new => NULL()
+      !!----------------------------------------------------------------------
+      !
+      new => NULL()
+      CALL create_iceberg( new, bergvals, ptvals )
+      CALL insert_berg_into_list( new )
+      new => NULL()     ! Clear new
+      !
+   END SUBROUTINE add_new_berg_to_list
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE create_iceberg(berg, bergvals, ptvals )
+      ! Arguments
+      TYPE(iceberg), INTENT(in) :: bergvals
+      TYPE(point)  , INTENT(in) :: ptvals
+      TYPE(iceberg), POINTER    :: berg
+      ! local variables
+      TYPE(point)  , POINTER    :: pt
+      INTEGER                   :: istat
+
+      IF ( ASSOCIATED(berg) ) then
+         CALL ctl_stop( 'icebergs, create_iceberg: berg already associated' )
+      ENDIF
+      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 create_iceberg
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE insert_berg_into_list( newberg )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE insert_berg_into_list  ***
+      !!
+      !! ** Purpose :   add a new berg to the iceberg list
+      !!
+      !! ** method  : - ?
+      !!----------------------------------------------------------------------
+      TYPE(iceberg), POINTER  ::   newberg
+      !
+      TYPE(iceberg), POINTER  ::   this, prev, last
+      !!----------------------------------------------------------------------
+
+      IF( ASSOCIATED( first_berg ) ) THEN
+!        last = last_berg()
+         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 insert_berg_into_list
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION yearday(imon, iday, ihr, imin, isec)
+      ! sga - improved but still only applies to 365 day year, need to do this properly
+      ! Arguments
+      INTEGER, intent(in)     :: imon, iday, ihr, imin, isec
+      ! Local variables
+      INTEGER                 :: i
+      INTEGER, DIMENSION(12)  :: months = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
+
+      yearday = FLOAT( SUM( months(1:imon) ) )
+      yearday = yearday + FLOAT(iday-1) + (FLOAT(ihr) + (FLOAT(imin) + FLOAT(isec)/60.)/60.)/24.
+
+   END FUNCTION yearday
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE delete_iceberg_from_list( first, berg )
+      ! Arguments
+      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
+
+      ! Bye-bye berg
+      CALL destroy_iceberg(berg)
+
+   END SUBROUTINE delete_iceberg_from_list
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE destroy_iceberg(berg)
+      ! Arguments
+      TYPE(iceberg), POINTER :: berg
+
+      ! Remove any points
+      IF( ASSOCIATED( berg%current_point ) ) DEALLOCATE( berg%current_point )
+
+      ! Bye-bye berg
+      DEALLOCATE(berg)
+
+   END SUBROUTINE destroy_iceberg
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE track_berg( num, label, kt )
+      ! Arguments
+      INTEGER, DIMENSION(nkounts)    :: num            ! iceberg number
+      CHARACTER(len=*)               :: label
+      INTEGER                        :: kt             ! timestep number
+      ! Local variables
+      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) /= num(k) ) match = .FALSE.
+         END DO
+         IF( match ) CALL print_berg(this, kt)
+         this => this%next
+      ENDDO
+
+   END SUBROUTINE track_berg
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE print_berg( berg, kt )
+      ! Arguments
+      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 print_berg
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE print_bergs( label, kt )
+      ! Arguments
+      CHARACTER(len=*)       :: label
+      INTEGER                :: kt             ! timestep number
+      ! Local variables
+      INTEGER                :: nbergs, nnbergs
+      TYPE(iceberg), POINTER :: this
+
+      this => first_berg
+      IF( ASSOCIATED(this) ) THEN
+         WRITE(numicb,'(a," pe=(",i3,")")' ) 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 print_berg(this, kt)
+        this => this%next
+      ENDDO
+      nbergs = count_bergs()
+      nnbergs = nbergs
+      IF( lk_mpp ) CALL mpp_sum(nnbergs)
+      IF ( nbergs .GT. 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &
+                                           label, nbergs, nnbergs, narea
+
+   END SUBROUTINE print_bergs
+
+   !!-------------------------------------------------------------------------
+
+   SUBROUTINE increment_kounter()
+      ! Small routine for coping with very large integer values labelling icebergs
+      ! kount_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
+
+      ! Local variables
+      INTEGER                :: i, ibig
+
+      ibig = HUGE(kount_bergs(1))
+      IF( ibig-jpnij < kount_bergs(1) ) THEN
+         kount_bergs(1) = narea
+         DO i = 2,nkounts
+            IF( kount_bergs(i) == ibig ) THEN
+               kount_bergs(i) = 0
+               IF( i == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
+            ELSE
+               kount_bergs(i) = kount_bergs(i) + 1
+               EXIT
+            ENDIF
+         END DO
+      ELSE
+         kount_bergs(1) = kount_bergs(1) + jpnij
+      ENDIF
+
+   END SUBROUTINE increment_kounter
+
+   !!-------------------------------------------------------------------------
+
+   INTEGER function count_bergs()
+      ! Local variables
+      TYPE(iceberg), POINTER :: this
+
+      count_bergs = 0
+      this => first_berg
+      DO WHILE( ASSOCIATED(this) )
+         count_bergs = count_bergs+1
+         this => this%next
+      ENDDO
+
+   END FUNCTION count_bergs
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION sum_mass( first, justbits, justbergs )
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION sum_mass  ***
+      !!
+      !! ** Purpose :   compute the mass all iceberg, all bergies or all bergs.
+      !!----------------------------------------------------------------------
+      TYPE(iceberg)      , POINTER  ::   first
+      TYPE(point)        , POINTER  ::   pt
+      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
+      !
+      TYPE(iceberg), POINTER ::   this
+      !!----------------------------------------------------------------------
+      sum_mass = 0._wp
+      this => first
+      !
+      IF( PRESENT( justbergs  ) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            sum_mass = sum_mass + pt%mass         * this%mass_scaling
+            this => this%next
+         END DO
+      ELSEIF( PRESENT(justbits) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            sum_mass = sum_mass + pt%mass_of_bits * this%mass_scaling
+            this => this%next
+         END DO
+      ELSE
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            sum_mass = sum_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
+            this => this%next
+         END DO
+      ENDIF
+      !
+   END FUNCTION sum_mass
+
+   !!-------------------------------------------------------------------------
+
+   REAL(wp) FUNCTION sum_heat( first, justbits, justbergs )
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION sum_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
+      !!----------------------------------------------------------------------
+      sum_heat = 0._wp
+      this => first
+      !
+      IF( PRESENT( justbergs  ) ) THEN
+         DO WHILE( ASSOCIATED( this ) )
+            pt => this%current_point
+            sum_heat = sum_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
+            sum_heat = sum_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
+            sum_heat = sum_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
+            this => this%next
+         END DO
+      ENDIF
+      !
+   END FUNCTION sum_heat
+
+   !!-------------------------------------------------------------------------
+
+END MODULE icbutl
