NEMOCFG/eORCA025/eORCA025.L121-Lemaitre4/cfgs/eORCA025_ICE/MY_SRC/icbdyn.F90

MODULE icbdyn
   !!======================================================================
   !!                       ***  MODULE  icbdyn  ***
   !! Iceberg:  time stepping routine for iceberg tracking
   !!======================================================================
   !! History :  3.3  !  2010-01  (Martin&Adcroft)  Original code
   !!             -   !  2011-03  (Madec)  Part conversion to NEMO form
   !!             -   !                    Removal of mapping from another grid
   !!             -   !  2011-04  (Alderson)  Split into separate modules
   !!             -   !  2011-05  (Alderson)  Replace broken grounding routine with one of
   !!             -   !                       Gurvan's suggestions (just like the broken one)
   !!----------------------------------------------------------------------
   USE par_oce        ! NEMO parameters
   USE dom_oce        ! NEMO ocean domain
   USE phycst         ! NEMO physical constants
   USE in_out_manager                      ! IO parameters
   !
   USE icb_oce        ! define iceberg arrays
   USE icbutl         ! iceberg utility routines
   USE icbdia         ! iceberg budget routines

   IMPLICIT NONE
   PRIVATE

   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: icbdyn.F90 15088 2021-07-06 13:03:34Z acc $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE icb_dyn( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_dyn  ***
      !!
      !! ** Purpose :   iceberg evolution.
      !!
      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   !
      !
      LOGICAL  ::   ll_bounced
      REAL(wp) ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
      REAL(wp) ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
      REAL(wp) ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
      REAL(wp) ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
      REAL(wp) ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
      REAL(wp) ::   zdt, zdt_2, zdt_6, ze1, ze2
      TYPE(iceberg), POINTER ::   berg
      TYPE(point)  , POINTER ::   pt
      !!----------------------------------------------------------------------
      !
      ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
      !                    with I.C.'s:   X=X1 and V=V1
      !
      !                                    ; A1=A(X1,V1)
      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
      !
      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6

      ! time steps
      zdt   = berg_dt
      zdt_2 = zdt * 0.5_wp
      zdt_6 = zdt / 6._wp

      berg => first_berg                    ! start from the first berg
      !
      DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
         !
         pt => berg%current_point

         ll_bounced = .FALSE.


         ! STEP 1 !
         ! ====== !
         zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
         zyj1 = pt%yj   ;   zvvel1 = pt%vvel


         !                                         !**   A1 = A(X1,V1)
         CALL icb_accel( kt, berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 0.5_wp )
         !
         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( berg, zxi2, zxi1, zu1,   &
            &                   zyj2, zyj1, zv1, ll_bounced )

         !                                         !**   A2 = A(X2,V2)
         CALL icb_accel( kt, berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 0.5_wp )
         !
         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( berg, zxi3, zxi1, zu2,   &
            &                   zyj3, zyj1, zv2, ll_bounced )

         !                                         !**   A3 = A(X3,V3)
         CALL icb_accel( kt, berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1._wp )
         !
         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( berg, zxi4, zxi1, zu3,   &
            &                   zyj4, zyj1, zv3, ll_bounced )

         !                                         !**   A4 = A(X4,V4)
         CALL icb_accel( kt, berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1._wp )

         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
         zv4 = zvvel4 / ze2

         ! FINAL STEP !
         ! ========== !
         !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
         !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
         zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
         zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
         zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 )
         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )

         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   &
            &                   zyj_n, zyj1, zvvel_n, ll_bounced )

         pt%uvel = zuvel_n                        !** save in berg structure
         pt%vvel = zvvel_n
         pt%xi   = zxi_n
         pt%yj   = zyj_n

         berg => berg%next                         ! switch to the next berg
         !
      END DO                                  !==  end loop over all bergs  ==!
      !
   END SUBROUTINE icb_dyn


   SUBROUTINE icb_ground( berg, pi, pi0, pu,   &
      &                         pj, pj0, pv, ld_bounced )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_ground  ***
      !!
      !! ** Purpose :   iceberg grounding.
      !!
      !! ** Method  : - adjust velocity and then put iceberg back to start position
      !!                NB two possibilities available one of which is hard-coded here
      !!----------------------------------------------------------------------
      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
      !
      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
      !
      INTEGER  ::   ii, ii0
      INTEGER  ::   ij, ij0
      INTEGER  ::   ikb
      INTEGER  ::   ibounce_method
      !
      REAL(wp) :: zD 
      REAL(wp), DIMENSION(jpk) :: ze3t
      !
      LOGICAL  :: licb_free
      !!----------------------------------------------------------------------
      !
      ld_bounced = .FALSE.
      !
      ii0 = INT( pi0+0.5 ) + (nn_hls-1)   ;   ij0 = INT( pj0+0.5 ) + (nn_hls-1)      ! initial gridpoint position (T-cell)
      ii  = INT( pi +0.5 ) + (nn_hls-1)   ;   ij  = INT( pj +0.5 ) + (nn_hls-1)      ! current     -         -
      !
      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
      !
      ! map into current processor
      ii0 = mi1( ii0 )
      ij0 = mj1( ij0 )
      ii  = mi1( ii  )
      ij  = mj1( ij  )
      !
      ! first check if enough room to go in the new cell (use before virtual area to avoid reproducibility issue)
      ! case icb stay in the same cell => icb free to move (prevent issue with icb with icb area > cell area).
      licb_free=.TRUE.
      IF ( ln_icb_area_mask ) THEN                                   ! 
         IF ( ii0 /= ii .OR. ij0 /= ij ) THEN                           ! icb enter in a new cell
            IF ( virtual_area_e(ii,ij) > e1e2t(ii,ij) ) licb_free=.FALSE. ! the new cell is full, icb cannot go in
         END IF
      END IF
         
      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
      IF ( ln_M2016 .AND. ln_icb_grd ) THEN
         !
         ! draught (keel depth)
         zD = rho_berg_1_oce * berg%current_point%thickness
         !
         ! interpol needed data
         CALL icb_utl_interp( pi, pj, pe3t=ze3t )
         ! 
         !compute bottom level
         CALL icb_utl_getkb( ikb, ze3t, zD )
         !
         ! berg reach a new t-cell, but an ocean one
         ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0)
         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp .AND. licb_free ) RETURN
         !
      ELSE
         IF(  tmask(ii,ij,1)  /=   0._wp .AND. licb_free )   RETURN           ! berg reach a new t-cell, but an ocean one
      END IF
      !
      ! From here, berg have reach land: treat grounding/bouncing
      ! -------------------------------
      ld_bounced = .TRUE.

      !! not obvious what should happen now
      !! if berg tries to enter a land box, the only location we can return it to is the start 
      !! position (pi0,pj0), since it has to be in a wet box to do any melting;
      !! first option is simply to set whole velocity to zero and move back to start point
      !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
      !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast

      ibounce_method = 2
      SELECT CASE ( ibounce_method )
      CASE ( 1 )
         pi = pi0
         pj = pj0
         pu = 0._wp
         pv = 0._wp
      CASE ( 2 )
         IF( ii0 /= ii ) THEN
            pi = pi0                   ! return back to the initial position
            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
         ENDIF
         IF( ij0 /= ij ) THEN
            pj = pj0                   ! return back to the initial position
            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
         ENDIF
      END SELECT
      !
   END SUBROUTINE icb_ground


   SUBROUTINE icb_accel( kt, berg , pxi, pe1, puvel, puvel0, pax,                 &
      &                             pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_accel  ***
      !!
      !! ** Purpose :   compute the iceberg acceleration.
      !!
      !! ** Method  : - sum the terms in the momentum budget
      !!----------------------------------------------------------------------
      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
      INTEGER                , INTENT(in   ) ::   kt               ! time step
      REAL(wp)               , INTENT(in   ) ::   pcfl_scale
      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
      !
      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
      !
      INTEGER  ::   itloop, ikb, jk
      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi
      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y
      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep
      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), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw
      !!----------------------------------------------------------------------

      ! Interpolate gridded fields to berg
      nknberg = berg%number(1)
      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor
         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient
         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis

      zM = berg%current_point%mass
      zT = berg%current_point%thickness               ! total thickness
      zD = rho_berg_1_oce * zT                        ! draught (keel depth)
      zF = zT - zD                                    ! freeboard
      zW = berg%current_point%width
      zL = berg%current_point%length

      zhi   = MIN( zhi   , zD    )
      zD_hi = MAX( 0._wp, zD-zhi )
 
     ! Wave radiation
      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model
      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
      zampl        = 0.5_wp * 0.02025_wp * zwmod      ! This is "a", the wave amplitude
      zLwavelength =       0.32_wp    * zwmod         ! Surface wave length fitted to data in table at
      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
      zLcutoff     = 0.125_wp * zLwavelength
      zLtop        = 0.25_wp  * zLwavelength
      zCr          = pp_Cr0 * MIN(  MAX( 0._wp, (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1._wp)  ! Wave radiation coefficient
      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
      zwave_rad    = 0.5_wp * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2._wp*zW*zL) / (zW+zL)
      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
      IF( zwmod /= 0._wp ) THEN
         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
         zvwave = zva/zwmod
      ELSE
         zuwave = 0._wp   ;    zvwave=0._wp   ;    zwave_rad=0._wp ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
      ENDIF

      ! Weighted drag coefficients
      z_ocn = pp_rho_seawater / zM * (0.5_wp*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
      z_atm = pp_rho_air      / zM * (0.5_wp*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
      z_ice = pp_rho_ice      / zM * (0.5_wp*pp_Cd_iv*zW*zhi              )
      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp

      ! lateral velocities
      ! default ssu and ssv
      ! ln_M2016: mean velocity along the profile
      IF ( ln_M2016 ) THEN
         ! interpol needed data
         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities
        
         !compute bottom level
         CALL icb_utl_getkb( ikb, ze3t, zD )
         
         ! compute mean velocity 
         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb)
         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb)
      ELSE
         zuo = zssu
         zvo = zssv
      END IF

      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
      !
      DO itloop = 1, 2  ! Iterate on drag coefficients
         !
         zus = 0.5_wp * ( zuveln + puvel )
         zvs = 0.5_wp * ( zvveln + pvvel )
         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
         !
         ! Explicit accelerations
         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
         zaxe = -grav * zssh_x + zwave_rad * zuwave
         zaye = -grav * zssh_y + zwave_rad * zvwave
         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
            zaxe = zaxe + zff*pvvel0
            zaye = zaye - zff*puvel0
         ELSE
            zaxe = zaxe + zff*pvvel
            zaye = zaye - zff*puvel
         ENDIF
         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
         ELSE
            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
         ENDIF

         ! Solve for implicit accelerations
         IF( pp_alpha + pp_beta > 0._wp ) THEN
            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
            zA11    = 1._wp + pp_beta *pdt*zlambda
            zA12    =         pp_alpha*pdt*zff
            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
         ELSE
            pax = zaxe   ;   pay = zaye
         ENDIF

         zuveln = puvel0 + pdt*pax
         zvveln = pvvel0 + pdt*pay
         !
      END DO      ! itloop

      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
         IF( zspeed > 0._wp ) THEN
            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing
            ! cfl scale is function of the RK4 step
            zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_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
               pax = (zuveln - puvel0)/pdt
               pay = (zvveln - pvvel0)/pdt
               !
               ! print speeding ticket
               IF (nn_verbose_level > 0) THEN
                  WRITE(numicb, 9200) 'icb speeding : ',kt, nknberg, zspeed, &
                       &                pxi, pyj, zuo, zvo, zua, zva, zui, zvi
                  9200 FORMAT(a,i9,i6,f9.2,1x,4(1x,2f9.2))
               END IF
               !
               CALL icb_dia_speed()
            ENDIF
         ENDIF
      ENDIF
      !                                      ! check the speed and acceleration limits
      IF (nn_verbose_level > 0) THEN
         IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
         IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
      ENDIF
      !
   END SUBROUTINE icb_accel

   !!======================================================================
END MODULE icbdyn