ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/ICB/icbdyn.F90

MODULE icbdyn

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

   IMPLICIT NONE
   PRIVATE

   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
   !! $Id: icbdyn.F90 2355 2015-05-20 07:11:50Z ufla $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE icb_dyn( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_dyn  ***
      !!
      !! ** Purpose :   iceberg evolution.
      !!
      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
      !!----------------------------------------------------------------------
      REAL(wp)                        ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
      REAL(wp)                        ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
      REAL(wp)                        ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
      REAL(wp)                        ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
      REAL(wp)                        ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
      REAL(wp)                        ::   zdt, zdt_2, zdt_6, ze1, ze2
      LOGICAL                         ::   ll_bounced
      TYPE(iceberg), POINTER          ::   berg
      TYPE(point)  , POINTER          ::   pt
      INTEGER                         ::   kt
      !!----------------------------------------------------------------------

      ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
      !                    with I.C.'s:   X=X1 and V=V1
      !
      !                                    ; A1=A(X1,V1)
      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
      !
      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6

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

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

         ll_bounced = .false.


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


         !                                         !**   A1 = A(X1,V1)
         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
         !
         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
         zv1 = zvvel1 / ze2

         ! STEP 2 !
         ! ====== !
         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
         ! position using di/dt & djdt   !   V2  in m/s
         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
         !
         CALL icb_ground( zxi2, zxi1, zu1,   &
         &                zyj2, zyj1, zv1, ll_bounced )

         !                                         !**   A2 = A(X2,V2)
         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
         !
         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
         zv2 = zvvel2 / ze2
         !
         ! STEP 3 !
         ! ====== !
         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
         !
         CALL icb_ground( zxi3, zxi1, zu3,   &
         &                zyj3, zyj1, zv3, ll_bounced )

         !                                         !**   A3 = A(X3,V3)
         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
         !
         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
         zv3 = zvvel3 / ze2

         ! STEP 4 !
         ! ====== !
         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3

         CALL icb_ground( zxi4, zxi1, zu4,   &
         &                zyj4, zyj1, zv4, ll_bounced )

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

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

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

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

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

         ! update actual position
         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )

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


   SUBROUTINE icb_ground( pi, pi0, pu,   &
      &                         pj, pj0, pv, ld_bounced )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_ground  ***
      !!
      !! ** Purpose :   iceberg grounding.
      !!
      !! ** Method  : - adjust velocity and then put iceberg back to start position
      !!                NB two possibilities available one of which is hard-coded here
      !!----------------------------------------------------------------------
      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
      !
      INTEGER  ::   ii, ii0
      INTEGER  ::   ij, ij0
      INTEGER  ::   ibounce_method
      !!----------------------------------------------------------------------
      !
      ld_bounced = .FALSE.
      !
      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
      !
      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
      !
      ! map into current processor
      ii0 = mi1( ii0 )
      ij0 = mj1( ij0 )
      ii  = mi1( ii  )
      ij  = mj1( ij  )
      !
      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
      !
      ! From here, berg have reach land: treat grounding/bouncing
      ! -------------------------------
      ld_bounced = .TRUE.

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

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


   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE icb_accel  ***
      !!
      !! ** Purpose :   compute the iceberg acceleration.
      !!
      !! ** Method  : - sum the terms in the momentum budget
      !!----------------------------------------------------------------------
      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
      !
      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
      !
      INTEGER  ::   itloop
      REAL(wp) ::   zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi
      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
      REAL(wp) ::   z_ocn, z_atm, z_ice
      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
      !!----------------------------------------------------------------------

      ! Interpolate gridded fields to berg
      nknberg = berg%number(1)
      CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     &
         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff )

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

      zhi   = MIN( zhi   , zD    )
      zD_hi = MAX( 0._wp, zD-zhi )

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

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

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

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

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

      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
         IF( zspeed > 0._wp ) THEN
            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
            zspeed_new = zloc_dx / pdt * rn_speed_limit     ! Speed limit as a factor of dx / dt
            IF( zspeed_new < zspeed ) THEN
               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
               CALL icb_dia_speed()
            ENDIF
         ENDIF
      ENDIF
      !                                      ! check the speed and acceleration limits
      IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
      IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
      !
   END SUBROUTINE icb_accel

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