Index: /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdia.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdia.F90	(revision 14273)
+++ /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdia.F90	(revision 14274)
@@ -83,8 +83,9 @@
    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 , DIMENSION(nclasses) ::  nbergs_calved_by_class
 
@@ -124,5 +125,6 @@
       nbergs_calved             = 0
       nbergs_calved_by_class(:) = 0
-      nspeeding_tickets         = 0
+      nspeeding_tickets(:)      = 0
+      vel_factor_min            = 1._wp
       stored_heat_end           = 0._wp
       floating_heat_end         = 0._wp
@@ -157,5 +159,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 +267,19 @@
             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 )
          ENDIF
          !
@@ -329,5 +332,8 @@
          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, '("min velocity reduction factor = ",f12.8)') vel_factor_min
+            ENDIF
          ENDIF
          !
@@ -337,5 +343,6 @@
          nbergs_calved             = 0
          nbergs_calved_by_class(:) = 0
-         nspeeding_tickets         = 0
+         nspeeding_tickets(:)      = 0
+         vel_factor_min            = 1._wp
          stored_heat_start         = stored_heat_end
          floating_heat_start       = floating_heat_end
@@ -474,10 +481,13 @@
 
 
-   SUBROUTINE icb_dia_speed()
-      !!----------------------------------------------------------------------
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT.ln_bergdia )   RETURN
-      nspeeding_tickets = nspeeding_tickets + 1
+   SUBROUTINE icb_dia_speed(pvel_factor, pn_stage)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pvel_factor   ! factor by which velocity reduced
+      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
       !
    END SUBROUTINE icb_dia_speed
Index: /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdyn.F90	(revision 14273)
+++ /NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdyn.F90	(revision 14274)
@@ -85,5 +85,5 @@
          !                                         !**   A1 = A(X1,V1)
          CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
-            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 1 )
          !
          zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
@@ -102,5 +102,5 @@
          !                                         !**   A2 = A(X2,V2)
          CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
-            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 2 )
          !
          zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
@@ -118,5 +118,5 @@
          !                                         !**   A3 = A(X3,V3)
          CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
-            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 3 )
          !
          zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
@@ -134,5 +134,5 @@
          !                                         !**   A4 = A(X4,V4)
          CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
-            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 4 )
 
          zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
@@ -235,5 +235,5 @@
 
    SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
-      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
+      &                         pyj, pe2, pvvel, pvvel0, pay, pdt, pstage )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_accel  ***
@@ -250,4 +250,5 @@
       REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
       REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
+      INTEGER                , INTENT(in   ) ::   pstage           ! RK4 stage (for speed limiter)
       !
       REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
@@ -265,5 +266,5 @@
       REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
       REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
-      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
+      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale
       !!----------------------------------------------------------------------
 
@@ -361,10 +362,21 @@
          zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
          IF( zspeed > 0._wp ) THEN
-            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
-            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
+            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing
+            ! cfl scale is function of the RK4 step
+            IF( pstage .le. 2 ) THEN
+               zcfl_scale=0.5
+            ELSE
+               zcfl_scale=1.0
+            ENDIF
+            zspeed_new = zloc_dx / pdt * rn_speed_limit * zcfl_scale ! Speed limit as a factor of dx / dt
             IF( zspeed_new < zspeed ) THEN
-               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
-               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
-               CALL icb_dia_speed()
+               zvel_factor = zspeed_new / zspeed
+               zuveln = zuveln * ( zvel_factor )             ! Scale velocity to reduce speed
+               zvveln = zvveln * ( zvel_factor )             ! without changing the direction
+               pax = (zuveln - puvel0)/pdt
+               pay = (zvveln - pvvel0)/pdt
+               WRITE(numicb,*) 'speeding ticket (zspeed_new/zspeed): ',zvel_factor, pdt, zcfl_scale
+               CALL FLUSH(numicb)
+               CALL icb_dia_speed(zvel_factor,pstage)
             ENDIF
          ENDIF
