Index: /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdia.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdia.F90	(revision 14166)
+++ /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdia.F90	(revision 14167)
@@ -83,8 +83,10 @@
    REAL(wp)                      ::  heat_to_bergs_net, heat_to_ocean_net, melt_net
    REAL(wp)                      ::  calving_to_bergs_net
+   REAL(wp)                      ::  vel_factor_min
 
    INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved
    INTEGER                       ::  nbergs_melted
-   INTEGER                       ::  nspeeding_tickets
+   INTEGER , DIMENSION(4)        ::  nspeeding_tickets
+   INTEGER                       ::  icount_max
    INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class
 
@@ -124,5 +126,7 @@
       nbergs_calved             = 0
       nbergs_calved_by_class(:) = 0
-      nspeeding_tickets         = 0
+      nspeeding_tickets(:)      = 0
+      icount_max                = 0
+      vel_factor_min            = 1._wp
       stored_heat_end           = 0._wp
       floating_heat_end         = 0._wp
@@ -157,5 +161,5 @@
       IF( lk_mpp ) THEN
          ALLOCATE( rsumbuf(23) )          ; rsumbuf(:) = 0._wp
-         ALLOCATE( nsumbuf(4+nclasses) )  ; nsumbuf(:) = 0
+         ALLOCATE( nsumbuf(7+nclasses) )  ; nsumbuf(:) = 0
          rsumbuf(1) = floating_mass_start
          rsumbuf(2) = bergs_mass_start
@@ -265,18 +269,20 @@
             nsumbuf(2) = nbergs_calved
             nsumbuf(3) = nbergs_melted
-            nsumbuf(4) = nspeeding_tickets
+            nsumbuf(4:7) = nspeeding_tickets(:)
             DO ik = 1, nclasses
-               nsumbuf(4+ik) = nbergs_calved_by_class(ik)
+               nsumbuf(7+ik) = nbergs_calved_by_class(ik)
             END DO
-            CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+4), nclasses+4 )
+            CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+7), nclasses+7 )
             !
             nbergs_end        = nsumbuf(1)
             nbergs_calved     = nsumbuf(2)
             nbergs_melted     = nsumbuf(3)
-            nspeeding_tickets = nsumbuf(4)
+            nspeeding_tickets(:) = nsumbuf(4:7)
             DO ik = 1,nclasses
-               nbergs_calved_by_class(ik)= nsumbuf(4+ik)
+               nbergs_calved_by_class(ik)= nsumbuf(7+ik)
             END DO
             !
+            CALL mpp_min( 'icbdia', vel_factor_min, 1 )
+            CALL mpp_max( 'icbdia', icount_max, 1 )
          ENDIF
          !
@@ -329,5 +335,9 @@
          IF (nn_verbose_level > 0) THEN
             WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses)
-            IF( nspeeding_tickets > 0 )   WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets
+            WRITE( numicb, '("n speeding tickets by RK4 stage = ",i6,3(",",i6))') (nspeeding_tickets(ik),ik=1,4)
+            IF( SUM(nspeeding_tickets) > 0 ) THEN
+               WRITE( numicb, '("max number of iterations = ",i6)') icount_max
+               WRITE( numicb, '("min velocity reduction factor = ",f12.8)') vel_factor_min
+            ENDIF
          ENDIF
          !
@@ -337,5 +347,7 @@
          nbergs_calved             = 0
          nbergs_calved_by_class(:) = 0
-         nspeeding_tickets         = 0
+         nspeeding_tickets(:)      = 0
+         icount_max                = 0
+         vel_factor_min            = 1._wp
          stored_heat_start         = stored_heat_end
          floating_heat_start       = floating_heat_end
@@ -474,10 +486,15 @@
 
 
-   SUBROUTINE icb_dia_speed()
-      !!----------------------------------------------------------------------
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT.ln_bergdia )   RETURN
-      nspeeding_tickets = nspeeding_tickets + 1
+   SUBROUTINE icb_dia_speed(pvel_factor, pn_count, pn_stage)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pvel_factor   ! factor by which velocity reduced
+      INTEGER , INTENT(in) ::   pn_count      ! which velocity-limiting stage are we on
+      INTEGER , INTENT(in) ::   pn_stage  ! which stage of the RK4 calculation are we on
+      !
+      IF( (.NOT.ln_bergdia) .OR. pn_stage .lt. 1 .OR. pn_stage .gt. 4 )   RETURN
+      nspeeding_tickets(pn_stage) = nspeeding_tickets(pn_stage) + 1
+      vel_factor_min = MIN(vel_factor_min,pvel_factor) ! keep track of minimum reduction factor
+      icount_max = MAX(icount_max,pn_count)              ! keep track of maximum number of iterations
       !
    END SUBROUTINE icb_dia_speed
Index: /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdyn.F90	(revision 14166)
+++ /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdyn.F90	(revision 14167)
@@ -41,5 +41,5 @@
       INTEGER, INTENT(in) ::   kt   !
       !
-      LOGICAL  ::   ll_bounced
+      LOGICAL  ::   ll_bounced, l_loop
       REAL(wp) ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
       REAL(wp) ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
@@ -47,6 +47,9 @@
       REAL(wp) ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
       REAL(wp) ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
+      REAL(wp) ::   zvel_factor, zuvel_rk4, zvvel_rk4
+      REAL(wp) ::   zvel_factor_keep(100)
       REAL(wp) ::   zspeed, zspeed_new, zloc_dx
       REAL(wp) ::   zdt, zdt_2, zdt_6, ze1, ze2
+      INTEGER  ::   icount
       TYPE(iceberg), POINTER ::   berg
       TYPE(point)  , POINTER ::   pt
@@ -77,91 +80,111 @@
          ll_bounced = .FALSE.
 
-
-         ! STEP 1 !
-         ! ====== !
-         zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
-         zyj1 = pt%yj   ;   zvvel1 = pt%vvel
-
-
-         !                                         !**   A1 = A(X1,V1)
-         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
-            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
-         !
-         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
-         zv1 = zvvel1 / ze2
-
-         ! STEP 2 !
-         ! ====== !
-         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
-         ! position using di/dt & djdt   !   V2  in m/s
-         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
-         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
-         !
-         CALL icb_ground( zxi2, zxi1, zu1,   &
-            &             zyj2, zyj1, zv1, ll_bounced )
-
-         !                                         !**   A2 = A(X2,V2)
-         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
-            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
-         !
-         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
-         zv2 = zvvel2 / ze2
-         !
-         ! STEP 3 !
-         ! ====== !
-         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
-         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
-         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
-         !
-         CALL icb_ground( zxi3, zxi1, zu3,   &
-            &                zyj3, zyj1, zv3, ll_bounced )
-
-         !                                         !**   A3 = A(X3,V3)
-         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
-            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
-         !
-         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
-         zv3 = zvvel3 / ze2
-
-         ! STEP 4 !
-         ! ====== !
-         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
-         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
-         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
-
-         CALL icb_ground( zxi4, zxi1, zu4,   &
-            &             zyj4, zyj1, zv4, ll_bounced )
-
-         !                                         !**   A4 = A(X4,V4)
-         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
-            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
-
-         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
-         zv4 = zvvel4 / ze2
-
-         ! FINAL STEP !
-         ! ========== !
-         !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
-         !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
-         zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
-         zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
-         zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 )
-         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
+         zvel_factor = 1.0
+         zvel_factor_keep(:) = 0.0
+         l_loop = .true.
+         icount = 0
+
+         DO WHILE(l_loop)
+
+            l_loop = .false.
+            icount=icount+1
+            zvel_factor_keep(icount) = zvel_factor
+
+            ! STEP 1 !
+            ! ====== !
+            zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
+            zyj1 = pt%yj   ;   zvvel1 = pt%vvel
+            !                                         !**   A1 = A(X1,V1)
+            CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
+                 &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+            !
+            zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
+            zv1 = zvvel1 / ze2
+
+            ! STEP 2 !
+            ! ====== !
+            !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
+            ! position using di/dt & djdt   !   V2  in m/s
+            zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zvel_factor * ( zuvel1 + zdt_2 * zax1 )
+            zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvel_factor * ( zvvel1 + zdt_2 * zay1 )
+            !
+            IF( rn_speed_limit > 0._wp ) THEN  ! Limit speed of bergs based on a CFL criteria (if asked)
+               CALL icb_speed_limit(ze1, ze2, zdt_2, zuvel2, zvvel2, rn_speed_limit, l_loop, zvel_factor, icount, 1)
+               IF(l_loop) CYCLE
+            ENDIF
+            !
+            CALL icb_ground( zxi2, zxi1, zu1,   &
+                 &             zyj2, zyj1, zv1, ll_bounced )
+
+            !                                         !**   A2 = A(X2,V2)
+            CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
+                 &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+            !
+            zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
+            zv2 = zvvel2 / ze2
+            !
+            ! STEP 3 !
+            ! ====== !
+            !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
+            zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zvel_factor * ( zuvel1 + zdt_2 * zax2 )
+            zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvel_factor * ( zvvel1 + zdt_2 * zay2 )
+            !
+            IF( rn_speed_limit > 0._wp ) THEN  ! Limit speed of bergs based on a CFL criteria (if asked)
+               CALL icb_speed_limit(ze1, ze2, zdt, zuvel3, zvvel3, rn_speed_limit, l_loop, zvel_factor, icount, 2)
+               IF(l_loop) CYCLE
+            ENDIF
+            !
+            CALL icb_ground( zxi3, zxi1, zu3,   &
+                 &                zyj3, zyj1, zv3, ll_bounced )
+
+            !                                         !**   A3 = A(X3,V3)
+            CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
+                 &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+            !
+            zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
+            zv3 = zvvel3 / ze2
+
+            ! STEP 4 !
+            ! ====== !
+            !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
+            zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zvel_factor * ( zuvel1 + zdt * zax3 )
+            zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvel_factor * ( zvvel1 + zdt * zay3 )
+            !
+            IF( rn_speed_limit > 0._wp ) THEN  ! Limit speed of bergs based on a CFL criteria (if asked)
+               zuvel_rk4 = zuvel1 + 2.*(zuvel2 + zuvel3) + zuvel4
+               zvvel_rk4 = zvvel1 + 2.*(zvvel2 + zvvel3) + zuvel4
+               CALL icb_speed_limit(ze1, ze2, zdt_6, zuvel_rk4, zvvel_rk4, rn_speed_limit, l_loop, zvel_factor, icount, 3)
+               IF(l_loop) CYCLE
+            ENDIF
+
+            CALL icb_ground( zxi4, zxi1, zu4,   &
+                 &           zyj4, zyj1, zv4, ll_bounced )
+
+            !                                         !**   A4 = A(X4,V4)
+            CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
+                 &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+
+            zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
+            zv4 = zvvel4 / ze2
+
+            ! FINAL STEP !
+            ! ========== !
+            !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
+            !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
+            zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
+            zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
+            zuvel_n = zvel_factor * ( pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 ) )
+            zvvel_n = zvel_factor * ( pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 ) )
+
+            IF( rn_speed_limit > 0._wp ) THEN  ! Limit speed of bergs based on a CFL criteria (if asked)
+               CALL icb_speed_limit(ze1, ze2, zdt_2, zuvel_n, zvvel_n, rn_speed_limit, l_loop, zvel_factor, icount, 4)
+            ENDIF
+
+         END DO
+
+         IF( icount > 5 ) WRITE(numicb, '("Large icount : vel_factors = ",i4,":",30(",",f12.8))') icount,zvel_factor_keep(1:icount)
 
          CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
             &             zyj_n, zyj1, zvvel_n, ll_bounced )
-
-         IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
-            zspeed = SQRT( zuvel_n*zuvel_n + zvvel_n*zvvel_n )    ! Speed of berg
-            IF( zspeed > 0._wp ) THEN
-               zloc_dx = MIN( ze1, ze2 )                          ! minimum grid spacing
-               zspeed_new = zloc_dx / zdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
-               IF( zspeed_new < zspeed ) THEN
-                  zuvel_n = zuvel_n * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
-                  zvvel_n = zvvel_n * ( zspeed_new / zspeed )        ! without changing the direction
-                  CALL icb_dia_speed()
-               ENDIF
-            ENDIF
-         ENDIF
 
          pt%uvel = zuvel_n                        !** save in berg structure
@@ -381,4 +404,38 @@
    END SUBROUTINE icb_accel
 
+   SUBROUTINE icb_speed_limit( pe1, pe2, pdt, pu, pv, pn_speed_limit, pl_loop, pvel_factor, pn_count, pn_stage )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE icb_speed_limit  ***
+      !!
+      !! ** Purpose :   calculate CFL-based speed limit for icebergs.
+      !!
+      !! ** Method  : - if iceberg velocity exceeds CFL-based limit, calculate
+      !!                a reduction factor for the velocity.
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pe1, pe2       ! scale factors for this iceberg grid cell
+      REAL(wp), INTENT(in   ) ::   pdt            ! timestep (for this RK4 stage)
+      REAL(wp), INTENT(inout) ::   pu  , pv       ! current iceberg velocities
+      REAL(wp), INTENT(in   ) ::   pn_speed_limit ! CFL-based speed limit
+      LOGICAL , INTENT(inout) ::   pl_loop        ! flag to trigger recalculation of RK4 step
+      REAL(wp), INTENT(inout) ::   pvel_factor    ! factor to multiply velocities
+      INTEGER,  INTENT(in   ) ::   pn_count       ! speed limiting iteration that we are on
+      INTEGER,  INTENT(in   ) ::   pn_stage   ! RK4 stage that we are on
+       !
+      REAL  ::   zspeed, zspeed_new, zloc_dx
+      !!----------------------------------------------------------------------
+
+      zspeed = SQRT( pu*pu + pv*pv )    ! Speed of berg
+      IF( zspeed > 0._wp ) THEN
+         zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
+         zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
+         IF( zspeed_new < zspeed ) THEN
+            pl_loop = .true.
+            ! arbitrary 0.99 factor to speed up convergence
+            pvel_factor = pvel_factor * 0.99_wp * zspeed_new/zspeed
+            CALL icb_dia_speed(pvel_factor, pn_count, pn_stage)
+         ENDIF
+      ENDIF
+   END SUBROUTINE icb_speed_limit
+
    !!======================================================================
 END MODULE icbdyn
