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