Forked from
Consortium Members / UKMO / GOSI / GOSI
223 commits behind the upstream repository.
-
Guillaume Samson authored89746a6d
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
icbdyn.F90 20.32 KiB
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
!!----------------------------------------------------------------------
!
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 )
!
! 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 ) RETURN
!
ELSE
IF( tmask(ii,ij,1) /= 0._wp ) 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