!=======================================================================
! Driver for ice mechanical redistribution (ridging)
!
! See these references:
!
! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength
! in modeling the thickness distribution of Arctic sea ice,
! J. Geophys. Res., 100, 18,611-18,626.
!
! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice
! cover, Mon. Wea. Rev., 108, 1943-1973, 1980.
!
! Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki, 2007:
! Improving ridging schemes for high-resolution sea ice models.
! J. Geophys. Res. 112, C03S91, doi:10.1029/2005JC003355.
!
! Rothrock, D. A., 1975: The energetics of the plastic deformation of
! pack ice by ridging, J. Geophys. Res., 80, 4514-4519.
!
! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,
! 1975: The thickness distribution of sea ice, J. Geophys. Res.,
! 80, 4501-4513.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
! 2004: Block structure added by William Lipscomb
! 2006: New options for participation and redistribution (WHL)
! 2006: Streamlined for efficiency by Elizabeth Hunke
! Converted to free source form (F90)
module icepack_mechred
use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c10, c25, Cf, Cp, Pstar, Cstar
use icepack_parameters, only: p05, p15, p25, p333, p5
use icepack_parameters, only: puny, Lfresh, rhoi, rhos
use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
use icepack_parameters, only: heat_capacity, conserv_check
use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice
use icepack_tracers, only: nt_apnd, nt_hpnd
use icepack_tracers, only: n_iso
use icepack_tracers, only: icepack_compute_tracers
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_itd, only: column_sum
use icepack_itd, only: column_conservation_check
use icepack_itd, only: cleanup_itd
implicit none
private
public :: ridge_ice, &
asum_ridging, &
ridge_itd, &
icepack_ice_strength, &
icepack_step_ridge
real (kind=dbl_kind), parameter :: &
exp_argmax = 100.0_dbl_kind, & ! maximum argument of exponential for underflow
Cs = p25 , & ! fraction of shear energy contrbtng to ridging
fsnowrdg = p5 , & ! snow fraction that survives in ridging
Gstar = p15 , & ! max value of G(h) that participates
! (krdg_partic = 0)
astar = p05 , & ! e-folding scale for G(h) participation
!echmod astar = p1 , & ! e-folding scale for G(h) participation
! (krdg_partic = 1)
maxraft= c1 , & ! max value of hrmin - hi = max thickness
! of ice that rafts (m)
Hstar = c25 ! determines mean thickness of ridged ice (m)
! (krdg_redist = 0)
! Flato & Hibler (1995) have Hstar = 100
!=======================================================================
contains
!=======================================================================
! Compute changes in the ice thickness distribution due to divergence
! and shear.
!
! author: William H. Lipscomb, LANL
subroutine ridge_ice (dt, ndtd, &
ncat, n_aero, &
nilyr, nslyr, &
ntrcr, hin_max, &
rdg_conv, rdg_shear, &
aicen, trcrn, &
vicen, vsnon, &
aice0, &
trcr_depend, trcr_base, &
n_trcr_strata, &
nt_strata, &
krdg_partic, krdg_redist,&
mu_rdg, tr_brine, &
dardg1dt, dardg2dt, &
dvirdgdt, opening, &
fpond, &
fresh, fhocn, &
faero_ocn, fiso_ocn, &
aparticn, krdgn, &
aredistn, vredistn, &
dardg1ndt, dardg2ndt, &
dvirdgndt, &
araftn, vraftn, &
closing_flag,closing )
integer (kind=int_kind), intent(in) :: &
ndtd , & ! number of dynamics subcycles
ncat , & ! number of thickness categories
nilyr , & ! number of ice layers
nslyr , & ! number of snow layers
n_aero, & ! number of aerosol tracers
ntrcr ! number of tracers in use
real (kind=dbl_kind), intent(in) :: &
mu_rdg , & ! gives e-folding scale of ridged ice (m^.5)
dt ! time step
real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max ! category limits (m)
real (kind=dbl_kind), intent(in) :: &
rdg_conv , & ! normalized energy dissipation due to convergence (1/s)
rdg_shear ! normalized energy dissipation due to shear (1/s)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn ! ice tracers
real (kind=dbl_kind), intent(inout) :: &
aice0 ! concentration of open water
integer (kind=int_kind), dimension (:), intent(in) :: &
trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
n_trcr_strata ! number of underlying tracer layers
real (kind=dbl_kind), dimension (:,:), intent(in) :: &
trcr_base ! = 0 or 1 depending on tracer dependency
! argument 2: (1) aice, (2) vice, (3) vsno
integer (kind=int_kind), dimension (:,:), intent(in) :: &
nt_strata ! indices of underlying tracer layers
integer (kind=int_kind), intent(in) :: &
krdg_partic, & ! selects participation function
krdg_redist ! selects redistribution function
logical (kind=log_kind), intent(in) :: &
closing_flag, &! flag if closing is valid
tr_brine ! if .true., brine height differs from ice thickness
! optional history fields
real (kind=dbl_kind), intent(inout), optional :: &
dardg1dt , & ! rate of fractional area loss by ridging ice (1/s)
dardg2dt , & ! rate of fractional area gain by new ridges (1/s)
dvirdgdt , & ! rate of ice volume ridged (m/s)
opening , & ! rate of opening due to divergence/shear (1/s)
closing , & ! rate of closing due to divergence/shear (1/s)
fpond , & ! fresh water flux to ponds (kg/m^2/s)
fresh , & ! fresh water flux to ocean (kg/m^2/s)
fhocn ! net heat flux to ocean (W/m^2)
real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s)
dardg2ndt , & ! rate of fractional area gain by new ridges (1/s)
dvirdgndt , & ! rate of ice volume ridged (m/s)
aparticn , & ! participation function
krdgn , & ! mean ridge thickness/thickness of ridging ice
araftn , & ! rafting ice area
vraftn , & ! rafting ice volume
aredistn , & ! redistribution function: fraction of new ridge area
vredistn ! redistribution function: fraction of new ridge volume
real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
faero_ocn ! aerosol flux to ocean (kg/m^2/s)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)
! local variables
real (kind=dbl_kind), dimension (ncat) :: &
eicen ! energy of melting for each ice layer (J/m^2)
real (kind=dbl_kind), dimension (ncat) :: &
esnon, & ! energy of melting for each snow layer (J/m^2)
vbrin, & ! ice volume with defined by brine height (m)
sicen ! Bulk salt in h ice (ppt*m)
real (kind=dbl_kind) :: &
asum , & ! sum of ice and open water area
aksum , & ! ratio of area removed to area ridged
msnow_mlt , & ! mass of snow added to ocean (kg m-2)
esnow_mlt , & ! energy needed to melt snow in ocean (J m-2)
mpond , & ! mass of pond added to ocean (kg m-2)
closing_net, & ! net rate at which area is removed (1/s)
! (ridging ice area - area of new ridges) / dt
divu_adv , & ! divu as implied by transport scheme (1/s)
opning , & ! rate of opening due to divergence/shear
! opning is a local variable;
! opening is the history diagnostic variable
ardg1 , & ! fractional area loss by ridging ice
ardg2 , & ! fractional area gain by new ridges
virdg , & ! ice volume ridged
aopen ! area opening due to divergence/shear
real (kind=dbl_kind), dimension (n_aero) :: &
maero ! aerosol mass added to ocean (kg m-2)
real (kind=dbl_kind), dimension (n_iso) :: &
miso ! isotope mass added to ocean (kg m-2)
real (kind=dbl_kind), dimension (0:ncat) :: &
apartic ! participation function; fraction of ridging
! and closing associated w/ category n
real (kind=dbl_kind), dimension (ncat) :: &
hrmin , & ! minimum ridge thickness
hrmax , & ! maximum ridge thickness (krdg_redist = 0)
hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
krdg , & ! mean ridge thickness/thickness of ridging ice
ardg1n , & ! area of ice ridged
ardg2n , & ! area of new ridges
virdgn , & ! ridging ice volume
mraftn ! rafting ice mask
real (kind=dbl_kind) :: &
vice_init, vice_final, & ! ice volume summed over categories
vsno_init, vsno_final, & ! snow volume summed over categories
eice_init, eice_final, & ! ice energy summed over layers
vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories
sice_init ,sice_final, & ! ice bulk salinity summed over categories
esno_init, esno_final ! snow energy summed over layers
integer (kind=int_kind), parameter :: &
nitermax = 20 ! max number of ridging iterations
integer (kind=int_kind) :: &
n , & ! thickness category index
niter , & ! iteration counter
k , & ! vertical index
it ! tracer index
real (kind=dbl_kind) :: &
dti ! 1 / dt
logical (kind=log_kind) :: &
iterate_ridging ! if true, repeat the ridging
character (len=char_len) :: &
fieldid ! field identifier
character(len=*),parameter :: subname='(ridge_ice)'
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
msnow_mlt = c0
esnow_mlt = c0
maero (:) = c0
miso (:) = c0
mpond = c0
ardg1 = c0
ardg2 = c0
virdg = c0
ardg1n(:) = c0
ardg2n(:) = c0
virdgn(:) = c0
mraftn(:) = c0
aopen = c0
!-----------------------------------------------------------------
! Compute area of ice plus open water before ridging.
!-----------------------------------------------------------------
call asum_ridging (ncat, aicen, aice0, asum)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Compute the area opening and closing.
!-----------------------------------------------------------------
if (closing_flag) then
opning = opening
closing_net = closing
divu_adv = opening - closing
else
call ridge_prep (dt, &
ncat, hin_max, &
rdg_conv, rdg_shear, &
asum, closing_net, &
divu_adv, opning)
endif
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Compute initial values of conserved quantities.
!-----------------------------------------------------------------
if (conserv_check) then
do n = 1, ncat
eicen(n) = c0
esnon(n) = c0
sicen(n) = c0
vbrin(n) = c0
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
do k = 1, nslyr
esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
* vsnon(n)/real(nslyr,kind=dbl_kind)
enddo
vbrin(n) = vicen(n)
if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vicen(n)
enddo ! n
call column_sum (ncat, &
vicen, vice_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
vsnon, vsno_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
eicen, eice_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
esnon, esno_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
sicen, sice_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
vbrin, vbri_init)
if (icepack_warnings_aborted(subname)) return
endif
rdg_iteration: do niter = 1, nitermax
!-----------------------------------------------------------------
! Compute the thickness distribution of ridging ice
! and various quantities associated with the new ridged ice.
!-----------------------------------------------------------------
call ridge_itd (ncat, aice0, &
aicen, vicen, &
krdg_partic, krdg_redist, &
mu_rdg, &
aksum, apartic, &
hrmin, hrmax, &
hrexp, krdg, &
aparticn, krdgn, &
mraftn)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Redistribute area, volume, and energy.
!-----------------------------------------------------------------
call ridge_shift (ntrcr, dt, &
ncat, hin_max, &
aicen, trcrn, &
vicen, vsnon, &
aice0, trcr_depend, &
trcr_base, n_trcr_strata,&
nt_strata, krdg_redist, &
aksum, apartic, &
hrmin, hrmax, &
hrexp, krdg, &
closing_net, opning, &
ardg1, ardg2, &
virdg, aopen, &
ardg1n, ardg2n, &
virdgn, &
nslyr, n_aero, &
msnow_mlt, esnow_mlt, &
maero, miso, &
mpond, &
aredistn, vredistn)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Make sure the new area = 1. If not (because the closing
! and opening rates were reduced above), prepare to ridge again
! with new rates.
!-----------------------------------------------------------------
call asum_ridging (ncat, aicen, aice0, asum)
if (icepack_warnings_aborted(subname)) return
if (abs(asum - c1) < puny) then
iterate_ridging = .false.
closing_net = c0
opning = c0
else
iterate_ridging = .true.
divu_adv = (c1 - asum) / dt
closing_net = max(c0, -divu_adv)
opning = max(c0, divu_adv)
endif
!-----------------------------------------------------------------
! If done, exit. If not, prepare to ridge again.
!-----------------------------------------------------------------
if (iterate_ridging) then
write(warnstr,*) subname, 'Repeat ridging, niter =', niter
call icepack_warnings_add(warnstr)
else
exit rdg_iteration
endif
if (niter == nitermax) then
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Exceeded max number of ridging iterations'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'max =',nitermax
call icepack_warnings_add(warnstr)
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" ridge_ice: Exceeded max number of ridging iterations" )
return
endif
enddo rdg_iteration ! niter
!-----------------------------------------------------------------
! Compute final values of conserved quantities.
! Check for conservation (allowing for snow thrown into ocean).
!-----------------------------------------------------------------
if (conserv_check) then
do n = 1, ncat
eicen(n) = c0
esnon(n) = c0
sicen(n) = c0
vbrin(n) = c0
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
do k = 1, nslyr
esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
* vsnon(n)/real(nslyr,kind=dbl_kind)
enddo
vbrin(n) = vicen(n)
if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vbrin(n)
enddo ! n
call column_sum (ncat, &
vicen, vice_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
vsnon, vsno_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
eicen, eice_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
esnon, esno_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
sicen, sice_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, &
vbrin, vbri_final)
if (icepack_warnings_aborted(subname)) return
vsno_final = vsno_final + msnow_mlt/rhos
esno_final = esno_final + esnow_mlt
fieldid = subname//':vice'
call column_conservation_check (fieldid, &
vice_init, vice_final, &
puny)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':vsno'
call column_conservation_check (fieldid, &
vsno_init, vsno_final, &
puny)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':eice'
call column_conservation_check (fieldid, &
eice_init, eice_final, &
puny*Lfresh*rhoi)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':esno'
call column_conservation_check (fieldid, &
esno_init, esno_final, &
puny*Lfresh*rhos)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':sice'
call column_conservation_check (fieldid, &
sice_init, sice_final, &
puny)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':vbrin'
call column_conservation_check (fieldid, &
vbri_init, vbri_final, &
puny*c10)
if (icepack_warnings_aborted(subname)) return
endif ! conserv_check
!-----------------------------------------------------------------
! Compute ridging diagnostics.
!-----------------------------------------------------------------
dti = c1/dt
if (present(dardg1dt)) then
dardg1dt = ardg1*dti
endif
if (present(dardg2dt)) then
dardg2dt = ardg2*dti
endif
if (present(dvirdgdt)) then
dvirdgdt = virdg*dti
endif
if (present(opening)) then
opening = aopen*dti
endif
if (present(dardg1ndt)) then
do n = 1, ncat
dardg1ndt(n) = ardg1n(n)*dti
enddo
endif
if (present(dardg2ndt)) then
do n = 1, ncat
dardg2ndt(n) = ardg2n(n)*dti
enddo
endif
if (present(dvirdgndt)) then
do n = 1, ncat
dvirdgndt(n) = virdgn(n)*dti
enddo
endif
if (present(araftn)) then
do n = 1, ncat
araftn(n) = mraftn(n)*ardg2n(n)
! araftn(n) = mraftn(n)*ardg1n(n)*p5
enddo
endif
if (present(vraftn)) then
do n = 1, ncat
vraftn(n) = mraftn(n)*virdgn(n)
enddo
endif
!-----------------------------------------------------------------
! Update fresh water and heat fluxes due to snow melt.
!-----------------------------------------------------------------
! use thermodynamic time step (ndtd*dt here) to average properly
dti = c1/(ndtd*dt)
if (present(fresh)) then
fresh = fresh + msnow_mlt*dti
endif
if (present(fhocn)) then
fhocn = fhocn + esnow_mlt*dti
endif
if (present(faero_ocn)) then
do it = 1, n_aero
faero_ocn(it) = faero_ocn(it) + maero(it)*dti
enddo
endif
if (tr_iso) then
! check size fiso_ocn vs n_iso ???
do it = 1, n_iso
fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti
enddo
endif
if (present(fpond)) then
fpond = fpond - mpond ! units change later
endif
!-----------------------------------------------------------------
! Check for fractional ice area > 1.
!-----------------------------------------------------------------
if (abs(asum - c1) > puny) then
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" total area > 1" )
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Ridging error: total area > 1'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'area:', asum
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'n, aicen:'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 0, aice0
call icepack_warnings_add(warnstr)
do n = 1, ncat
write(warnstr,*) subname, n, aicen(n)
call icepack_warnings_add(warnstr)
enddo
return
endif
end subroutine ridge_ice
!=======================================================================
! Find the total area of ice plus open water in each grid cell.
!
! This is similar to the aggregate_area subroutine except that the
! total area can be greater than 1, so the open water area is
! included in the sum instead of being computed as a residual.
!
! author: William H. Lipscomb, LANL
subroutine asum_ridging (ncat, aicen, aice0, asum)
integer (kind=int_kind), intent(in) :: &
ncat ! number of thickness categories
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen ! concentration of ice in each category
real (kind=dbl_kind), intent(in) :: &
aice0 ! concentration of open water
real (kind=dbl_kind), intent(out):: &
asum ! sum of ice and open water area
! local variables
integer (kind=int_kind) :: n
character(len=*),parameter :: subname='(asum_ridging)'
asum = aice0
do n = 1, ncat
asum = asum + aicen(n)
enddo
end subroutine asum_ridging
!=======================================================================
! Initialize arrays, compute area of closing and opening
!
! author: William H. Lipscomb, LANL
subroutine ridge_prep (dt, &
ncat, hin_max, &
rdg_conv, rdg_shear, &
asum, closing_net, &
divu_adv, opning)
integer (kind=int_kind), intent(in) :: &
ncat ! number of thickness categories
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max ! category limits (m)
real (kind=dbl_kind), intent(in) :: &
rdg_conv , & ! normalized energy dissipation due to convergence (1/s)
rdg_shear ! normalized energy dissipation due to shear (1/s)
real (kind=dbl_kind), intent(inout):: &
asum ! sum of ice and open water area
real (kind=dbl_kind), &
intent(out):: &
closing_net, & ! net rate at which area is removed (1/s)
divu_adv , & ! divu as implied by transport scheme (1/s)
opning ! rate of opening due to divergence/shear
! local variables
real (kind=dbl_kind), parameter :: &
big = 1.0e+8_dbl_kind
character(len=*),parameter :: subname='(ridge_prep)'
! Set hin_max(ncat) to a big value to ensure that all ridged ice
! is thinner than hin_max(ncat).
hin_max(ncat) = big
!-----------------------------------------------------------------
! Compute the net rate of closing due to convergence
! and shear, based on Flato and Hibler (1995).
!
! For the elliptical yield curve:
! rdg_conv = -min (divu, 0)
! rdg_shear = (1/2) * (Delta - abs(divu))
! Note that the shear term also accounts for divergence.
!
! The energy dissipation rate is equal to the net closing rate
! times the ice strength.
!
! NOTE: The NET closing rate is equal to the rate that open water
! area is removed, plus the rate at which ice area is removed by
! ridging, minus the rate at which area is added in new ridges.
! The GROSS closing rate is equal to the first two terms (open
! water closing and thin ice ridging) without the third term
! (thick, newly ridged ice).
!
! rdg_conv is calculated differently in EAP (update_ice_rdg) and
! represents closing_net directly. In that case, rdg_shear=0.
!-----------------------------------------------------------------
closing_net = Cs*rdg_shear + rdg_conv
!-----------------------------------------------------------------
! Compute divu_adv, the divergence rate given by the transport/
! advection scheme, which may not be equal to divu as computed
! from the velocity field.
!
! If divu_adv < 0, make sure the closing rate is large enough
! to give asum = 1.0 after ridging.
!-----------------------------------------------------------------
divu_adv = (c1-asum) / dt
if (divu_adv < c0) closing_net = max(closing_net, -divu_adv)
!-----------------------------------------------------------------
! Compute the (non-negative) opening rate that will give
! asum = 1.0 after ridging.
!-----------------------------------------------------------------
opning = closing_net + divu_adv
end subroutine ridge_prep
!=======================================================================
! Compute the thickness distribution of the ice and open water
! participating in ridging and of the resulting ridges.
!
! This version includes new options for ridging participation and
! redistribution.
! The new participation scheme (krdg_partic = 1) improves stability
! by increasing the time scale for large changes in ice strength.
! The new exponential redistribution function (krdg_redist = 1) improves
! agreement between ITDs of modeled and observed ridges.
!
! author: William H. Lipscomb, LANL
!
! 2006: Changed subroutine name to ridge_itd
! Added new options for ridging participation and redistribution.
subroutine ridge_itd (ncat, aice0, &
aicen, vicen, &
krdg_partic, krdg_redist, &
mu_rdg, &
aksum, apartic, &
hrmin, hrmax, &
hrexp, krdg, &
aparticn, krdgn, &
mraft)
integer (kind=int_kind), intent(in) :: &
ncat ! number of thickness categories
real (kind=dbl_kind), intent(in) :: &
mu_rdg , & ! gives e-folding scale of ridged ice (m^.5)
aice0 ! concentration of open water
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
integer (kind=int_kind), intent(in) :: &
krdg_partic , & ! selects participation function
krdg_redist ! selects redistribution function
real (kind=dbl_kind), intent(out):: &
aksum ! ratio of area removed to area ridged
real (kind=dbl_kind), dimension (0:ncat), intent(out) :: &
apartic ! participation function; fraction of ridging
! and closing associated w/ category n
real (kind=dbl_kind), dimension (:), intent(out) :: &
hrmin , & ! minimum ridge thickness
hrmax , & ! maximum ridge thickness (krdg_redist = 0)
hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
krdg ! mean ridge thickness/thickness of ridging ice
! diagnostic, category values
real (kind=dbl_kind), dimension(:), intent(out), optional :: &
aparticn, & ! participation function
krdgn ! mean ridge thickness/thickness of ridging ice
real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
mraft ! rafting ice mask
! local variables
integer (kind=int_kind) :: &
n ! thickness category index
real (kind=dbl_kind), parameter :: &
Gstari = c1/Gstar, &
astari = c1/astar
real (kind=dbl_kind), dimension(-1:ncat) :: &
Gsum ! Gsum(n) = sum of areas in categories 0 to n
real (kind=dbl_kind) :: &
work ! temporary work variable
real (kind=dbl_kind) :: &
hi , & ! ice thickness for each cat (m)
hrmean , & ! mean ridge thickness (m)
xtmp ! temporary variable
character(len=*),parameter :: subname='(ridge_itd)'
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
Gsum (-1:ncat) = c0 ! initialize
Gsum (-1) = c0 ! by definition
if (aice0 > puny) then
Gsum(0) = aice0
else
Gsum(0) = Gsum(-1)
endif
apartic(0) = c0
do n = 1, ncat
apartic(n) = c0
hrmin (n) = c0
hrmax (n) = c0
hrexp (n) = c0
krdg (n) = c1
!-----------------------------------------------------------------
! Compute the thickness distribution of ice participating in ridging.
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! First compute the cumulative thickness distribution function Gsum,
! where Gsum(n) is the fractional area in categories 0 to n.
! Ignore categories with very small areas.
!-----------------------------------------------------------------
if (aicen(n) > puny) then
Gsum(n) = Gsum(n-1) + aicen(n)
else
Gsum(n) = Gsum(n-1)
endif
enddo
! normalize
if (Gsum(ncat) > c0) then
work = c1 / Gsum(ncat)
do n = 0, ncat
Gsum(n) = Gsum(n) * work
enddo
endif
!-----------------------------------------------------------------
! Compute the participation function apartic; this is analogous to
! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
!
! area lost from category n due to ridging/closing
! apartic(n) = --------------------------------------------------
! total area lost due to ridging/closing
!
!-----------------------------------------------------------------
if (krdg_partic == 0) then ! Thornike et al. 1975 formulation
!-----------------------------------------------------------------
! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
! The expressions for apartic are found by integrating b(h)g(h) between
! the category boundaries.
!-----------------------------------------------------------------
do n = 0, ncat
if (Gsum(n) < Gstar) then
apartic(n) = Gstari*(Gsum(n ) - Gsum(n-1)) * &
(c2 - Gstari*(Gsum(n-1) + Gsum(n )))
elseif (Gsum(n-1) < Gstar) then
apartic(n) = Gstari*(Gstar - Gsum(n-1)) * &
(c2 - Gstari*(Gstar + Gsum(n-1)))
endif
enddo ! n
elseif (krdg_partic==1) then ! exponential dependence on G(h)
!-----------------------------------------------------------------
! b(h) = exp(-G(h)/astar)
! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)].
! The expression for apartic is found by integrating b(h)g(h)
! between the category boundaries.
!-----------------------------------------------------------------
! precompute exponential terms using Gsum as work array
xtmp = c1 / (c1 - exp(-astari))
Gsum(-1) = exp(-Gsum(-1)*astari) * xtmp
do n = 0, ncat
Gsum(n) = exp(-Gsum(n)*astari) * xtmp
apartic(n) = Gsum(n-1) - Gsum(n)
enddo ! n
endif ! krdg_partic
!-----------------------------------------------------------------
! Compute variables related to ITD of ridged ice:
!
! krdg = mean ridge thickness / thickness of ridging ice
! hrmin = min ridge thickness
! hrmax = max ridge thickness (krdg_redist = 0)
! hrexp = ridge e-folding scale (krdg_redist = 1)
!----------------------------------------------------------------
if (krdg_redist == 0) then ! Hibler 1980 formulation
!-----------------------------------------------------------------
! Assume ridged ice is uniformly distributed between hrmin and hrmax.
!
! This parameterization is a modified version of Hibler (1980).
! In the original paper the min ridging thickness is hrmin = 2*hi,
! and the max thickness is hrmax = 2*sqrt(hi*Hstar).
!
! Here the min thickness is hrmin = min(2*hi, hi+maxraft),
! so thick ridging ice is not required to raft.
!
!-----------------------------------------------------------------
do n = 1, ncat
if (aicen(n) > puny) then
hi = vicen(n) / aicen(n)
hrmin(n) = min(c2*hi, hi + maxraft)
hrmax(n) = c2*sqrt(Hstar*hi)
hrmax(n) = max(hrmax(n), hrmin(n)+puny)
hrmean = p5 * (hrmin(n) + hrmax(n))
krdg(n) = hrmean / hi
! diagnostic rafting mask not implemented
endif
enddo ! n
else ! krdg_redist = 1; exponential redistribution
!-----------------------------------------------------------------
! The ridge ITD is a negative exponential:
!
! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin
!
! where hrmin is the minimum thickness of ridging ice and
! hrexp is the e-folding thickness.
!
! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
! That is, the minimum ridge thickness results from rafting,
! unless the ice is thicker than maxraft.
!
! Also, assume that hrexp = mu_rdg*sqrt(hi).
! The parameter mu_rdg is tuned to give e-folding scales mostly
! in the range 2-4 m as observed by upward-looking sonar.
!
! Values of mu_rdg in the right column give ice strengths
! roughly equal to values of Hstar in the left column
! (within ~10 kN/m for typical ITDs):
!
! Hstar mu_rdg
!
! 25 3.0
! 50 4.0
! 75 5.0
! 100 6.0
!-----------------------------------------------------------------
do n = 1, ncat
if (aicen(n) > puny) then
hi = vicen(n) / aicen(n)
hi = max(hi,puny)
hrmin(n) = min(c2*hi, hi + maxraft)
hrexp(n) = mu_rdg * sqrt(hi)
krdg(n) = (hrmin(n) + hrexp(n)) / hi
!echmod: check computational efficiency
! diagnostic rafting mask
if (present(mraft)) then
mraft(n) = max(c0, sign(c1, hi+maxraft-hrmin(n)))
xtmp = mraft(n)*((c2*hi+hrexp(n))/hi - krdg(n))
mraft(n) = max(c0, sign(c1, puny-abs(xtmp)))
endif
endif
enddo
endif ! krdg_redist
!----------------------------------------------------------------
! Compute aksum = net ice area removed / total area participating.
! For instance, if a unit area of ice with h = 1 participates in
! ridging to form a ridge with a = 1/3 and h = 3, then
! aksum = 1 - 1/3 = 2/3.
!----------------------------------------------------------------
aksum = apartic(0) ! area participating = area removed
do n = 1, ncat
! area participating > area removed
aksum = aksum + apartic(n) * (c1 - c1/krdg(n))
enddo
! diagnostics
if (present(aparticn)) then
do n = 1, ncat
aparticn(n) = apartic(n)
enddo
endif
if (present(krdgn)) then
do n = 1, ncat
krdgn(n) = krdg(n)
enddo
endif
end subroutine ridge_itd
!=======================================================================
! Remove area, volume, and energy from each ridging category
! and add to thicker ice categories.
!
! Tracers: Ridging conserves ice volume and therefore conserves volume
! tracers. It does not conserve ice area, and therefore a portion of area
! tracers are lost (corresponding to the net closing). Area tracers on
! ice that participates in ridging are carried onto the resulting ridged
! ice (except the portion that are lost due to closing). Therefore,
! tracers must be decremented if they are lost to the ocean during ridging
! (e.g. snow, ponds) or if they are being carried only on the level ice
! area.
!
! author: William H. Lipscomb, LANL
subroutine ridge_shift (ntrcr, dt, &
ncat, hin_max, &
aicen, trcrn, &
vicen, vsnon, &
aice0, trcr_depend, &
trcr_base, n_trcr_strata, &
nt_strata, krdg_redist, &
aksum, apartic, &
hrmin, hrmax, &
hrexp, krdg, &
closing_net, opning, &
ardg1, ardg2, &
virdg, aopen, &
ardg1nn, ardg2nn, &
virdgnn, &
nslyr, n_aero, &
msnow_mlt, esnow_mlt, &
maero, miso, &
mpond, &
aredistn, vredistn)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nslyr , & ! number of snow layers
ntrcr , & ! number of tracers in use
n_aero, & ! number of aerosol tracers
krdg_redist ! selects redistribution function
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
integer (kind=int_kind), dimension (:), intent(in) :: &
trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
n_trcr_strata ! number of underlying tracer layers
real (kind=dbl_kind), dimension (:,:), intent(in) :: &
trcr_base ! = 0 or 1 depending on tracer dependency
! argument 2: (1) aice, (2) vice, (3) vsno
integer (kind=int_kind), dimension (:,:), intent(in) :: &
nt_strata ! indices of underlying tracer layers
real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
hin_max ! category limits (m)
real (kind=dbl_kind), intent(inout) :: &
aice0 ! concentration of open water
real (kind=dbl_kind), dimension (:), intent(inout) :: &
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn ! ice tracers
real (kind=dbl_kind), intent(in) :: &
aksum ! ratio of area removed to area ridged
real (kind=dbl_kind), dimension (0:ncat), intent(in) :: &
apartic ! participation function; fraction of ridging
! and closing associated w/ category n
real (kind=dbl_kind), dimension (:), intent(in) :: &
hrmin , & ! minimum ridge thickness
hrmax , & ! maximum ridge thickness (krdg_redist = 0)
hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
krdg ! mean ridge thickness/thickness of ridging ice
real (kind=dbl_kind), intent(inout) :: &
closing_net, & ! net rate at which area is removed (1/s)
opning , & ! rate of opening due to divergence/shear (1/s)
ardg1 , & ! fractional area loss by ridging ice
ardg2 , & ! fractional area gain by new ridges
virdg , & ! ice volume ridged (m)
aopen ! area opened due to divergence/shear
real (kind=dbl_kind), dimension(:), intent(inout) :: &
ardg1nn , & ! area of ice ridged
ardg2nn , & ! area of new ridges
virdgnn ! ridging ice volume
real (kind=dbl_kind), intent(inout) :: &
msnow_mlt , & ! mass of snow added to ocean (kg m-2)
esnow_mlt , & ! energy needed to melt snow in ocean (J m-2)
mpond ! mass of pond added to ocean (kg m-2)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
maero ! aerosol mass added to ocean (kg m-2)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
miso ! isotope mass added to ocean (kg m-2)
real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
aredistn , & ! redistribution function: fraction of new ridge area
vredistn ! redistribution function: fraction of new ridge volume
! local variables
integer (kind=int_kind) :: &
n, nr , & ! thickness category indices
k , & ! ice layer index
it , & ! tracer index
ntr , & ! tracer index
itl ! loop index
real (kind=dbl_kind), dimension (ncat) :: &
aicen_init , & ! ice area before ridging
vicen_init , & ! ice volume before ridging
vsnon_init ! snow volume before ridging
real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
atrcrn ! aicen*trcrn
real (kind=dbl_kind), dimension(3) :: &
trfactor ! base quantity on which tracers are carried
real (kind=dbl_kind) :: &
work , & ! temporary variable
expL_arg , & ! temporary exp arg values
expR_arg , & ! temporary exp arg values
closing_gross ! rate at which area removed, not counting
! area of new ridges
! ECH note: the following arrays only need be defined on iridge cells
real (kind=dbl_kind) :: &
afrac , & ! fraction of category area ridged
ardg1n , & ! area of ice ridged
ardg2n , & ! area of new ridges
virdgn , & ! ridging ice volume
vsrdgn , & ! ridging snow volume
dhr , & ! hrmax - hrmin
dhr2 , & ! hrmax^2 - hrmin^2
farea , & ! fraction of new ridge area going to nr
fvol ! fraction of new ridge volume going to nr
real (kind=dbl_kind) :: &
esrdgn ! ridging snow energy
real (kind=dbl_kind) :: &
hi1 , & ! thickness of ridging ice
hexp , & ! ridge e-folding thickness
hL, hR , & ! left and right limits of integration
expL, expR , & ! exponentials involving hL, hR
tmpfac , & ! factor by which opening/closing rates are cut
wk1 ! work variable
character(len=*),parameter :: subname='(ridge_shift)'
do n = 1, ncat
!-----------------------------------------------------------------
! Save initial state variables
!-----------------------------------------------------------------
aicen_init(n) = aicen(n)
vicen_init(n) = vicen(n)
vsnon_init(n) = vsnon(n)
!-----------------------------------------------------------------
! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
!-----------------------------------------------------------------
do it = 1, ntrcr
atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
+ trcr_base(it,2) * vicen(n) &
+ trcr_base(it,3) * vsnon(n))
if (n_trcr_strata(it) > 0) then ! additional tracer layers
do itl = 1, n_trcr_strata(it)
ntr = nt_strata(it,itl)
atrcrn(it,n) = atrcrn(it,n) * trcrn(ntr,n)
enddo
endif
enddo
enddo ! ncat
!-----------------------------------------------------------------
! Based on the ITD of ridging and ridged ice, convert the net
! closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
!-----------------------------------------------------------------
closing_gross = closing_net / aksum
!-----------------------------------------------------------------
! Reduce the closing rate if more than 100% of the open water
! would be removed. Reduce the opening rate proportionately.
!-----------------------------------------------------------------
if (apartic(0) > c0) then
wk1 = apartic(0) * closing_gross * dt
if (wk1 > aice0) then
tmpfac = aice0 / wk1
closing_gross = closing_gross * tmpfac
opning = opning * tmpfac
endif
endif
!-----------------------------------------------------------------
! Reduce the closing rate if more than 100% of any ice category
! would be removed. Reduce the opening rate proportionately.
!-----------------------------------------------------------------
do n = 1, ncat
if (aicen(n) > puny .and. apartic(n) > c0) then
wk1 = apartic(n) * closing_gross * dt
if (wk1 > aicen(n)) then
tmpfac = aicen(n) / wk1
closing_gross = closing_gross * tmpfac
opning = opning * tmpfac
endif
endif
enddo ! n
!-----------------------------------------------------------------
! Compute change in open water area due to closing and opening.
!-----------------------------------------------------------------
aice0 = aice0 - apartic(0)*closing_gross*dt + opning*dt
if (aice0 < -puny) then
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//' Ridging error: aice0 < 0')
write(warnstr,*) subname, 'aice0:', aice0
call icepack_warnings_add(warnstr)
return
elseif (aice0 < c0) then ! roundoff error
aice0 = c0
endif
aopen = opning*dt ! optional diagnostic
!-----------------------------------------------------------------
! Compute the area, volume, and energy of ice ridging in each
! category, along with the area of the resulting ridge.
!-----------------------------------------------------------------
do n = 1, ncat
!-----------------------------------------------------------------
! Identify grid cells with nonzero ridging
!-----------------------------------------------------------------
if (aicen_init(n) > puny .and. apartic(n) > c0 &
.and. closing_gross > c0) then
!-----------------------------------------------------------------
! Compute area of ridging ice (ardg1n) and of new ridge (ardg2n).
! Make sure ridging fraction <=1. (Roundoff errors can give
! ardg1 slightly greater than aicen.)
!-----------------------------------------------------------------
ardg1n = apartic(n)*closing_gross*dt
if (ardg1n > aicen_init(n) + puny) then
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//' Ridging error: ardg > aicen')
write(warnstr,*) subname, 'n, ardg, aicen:', &
n, ardg1n, aicen_init(n)
call icepack_warnings_add(warnstr)
return
else
ardg1n = min(aicen_init(n), ardg1n)
endif
ardg2n = ardg1n / krdg(n)
afrac = ardg1n / aicen_init(n)
!-----------------------------------------------------------------
! Subtract area, volume, and energy from ridging category n.
! Note: Tracer values are unchanged.
!-----------------------------------------------------------------
virdgn = vicen_init(n) * afrac
vsrdgn = vsnon_init(n) * afrac
aicen(n) = aicen(n) - ardg1n
vicen(n) = vicen(n) - virdgn
vsnon(n) = vsnon(n) - vsrdgn
!-----------------------------------------------------------------
! Increment ridging diagnostics
!-----------------------------------------------------------------
ardg1 = ardg1 + ardg1n
ardg2 = ardg2 + ardg2n
virdg = virdg + virdgn
ardg1nn(n) = ardg1n
ardg2nn(n) = ardg2n
virdgnn(n) = virdgn
!-----------------------------------------------------------------
! Place part of the snow and tracer lost by ridging into the ocean.
!-----------------------------------------------------------------
msnow_mlt = msnow_mlt + rhos*vsrdgn*(c1-fsnowrdg)
if (tr_aero) then
do it = 1, n_aero
maero(it) = maero(it) &
+ vsrdgn*(c1-fsnowrdg) &
*(trcrn(nt_aero +4*(it-1),n) &
+ trcrn(nt_aero+1+4*(it-1),n))
enddo
endif
if (tr_iso) then
do it = 1, n_iso
miso(it) = miso(it) + vsrdgn*(c1-fsnowrdg) &
* (trcrn(nt_isosno+it-1,n) &
+ trcrn(nt_isoice+it-1,n))
enddo
endif
if (tr_pond_topo) then
mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
* trcrn(nt_hpnd,n)
endif
!-----------------------------------------------------------------
! Compute quantities used to apportion ice among categories
! in the nr loop below
!-----------------------------------------------------------------
dhr = hrmax(n) - hrmin(n)
dhr2 = hrmax(n) * hrmax(n) - hrmin(n) * hrmin(n)
!-----------------------------------------------------------------
! Increment energy needed to melt snow in ocean.
! Note that esnow_mlt < 0; the ocean must cool to melt snow.
!-----------------------------------------------------------------
do k = 1, nslyr
esrdgn = vsrdgn * trcrn(nt_qsno+k-1,n) &
/ real(nslyr,kind=dbl_kind)
esnow_mlt = esnow_mlt + esrdgn*(c1-fsnowrdg)
enddo
!-----------------------------------------------------------------
! Subtract area- and volume-weighted tracers from category n.
!-----------------------------------------------------------------
do it = 1, ntrcr
trfactor(1) = trcr_base(it,1)*ardg1n
trfactor(2) = trcr_base(it,2)*virdgn
trfactor(3) = trcr_base(it,3)*vsrdgn
work = c0
do k = 1, 3
work = work + trfactor(k)*trcrn(it,n)
enddo
if (n_trcr_strata(it) > 0) then ! additional tracer layers
do itl = 1, n_trcr_strata(it)
ntr = nt_strata(it,itl)
work = work * trcrn(ntr,n)
enddo
endif
atrcrn(it,n) = atrcrn(it,n) - work
enddo ! ntrcr
!-----------------------------------------------------------------
! Add area, volume, and energy of new ridge to each category nr.
!-----------------------------------------------------------------
do nr = 1, ncat
if (krdg_redist == 0) then ! Hibler 1980 formulation
!-----------------------------------------------------------------
! Compute the fraction of ridged ice area and volume going to
! thickness category nr.
!-----------------------------------------------------------------
if (hrmin(n) >= hin_max(nr) .or. &
hrmax(n) <= hin_max(nr-1)) then
hL = c0
hR = c0
else
hL = max (hrmin(n), hin_max(nr-1))
hR = min (hrmax(n), hin_max(nr))
endif
farea = (hR-hL) / dhr
fvol = (hR*hR - hL*hL) / dhr2
else ! krdg_redist = 1; 2005 exponential formulation
!-----------------------------------------------------------------
! Compute the fraction of ridged ice area and volume going to
! thickness category nr.
!-----------------------------------------------------------------
if (nr < ncat) then
hi1 = hrmin(n)
hexp = hrexp(n)
if (hi1 >= hin_max(nr)) then
farea = c0
fvol = c0
else
hL = max (hi1, hin_max(nr-1))
hR = hin_max(nr)
expL_arg = min(((hL-hi1)/hexp),exp_argmax)
expR_arg = min(((hR-hi1)/hexp),exp_argmax)
expL = exp(-(expL_arg))
expR = exp(-(expR_arg))
farea = expL - expR
fvol = ((hL + hexp)*expL &
- (hR + hexp)*expR) / (hi1 + hexp)
endif
else ! nr = ncat
hi1 = hrmin(n)
hexp = hrexp(n)
hL = max (hi1, hin_max(nr-1))
expL_arg = min(((hL-hi1)/hexp),exp_argmax)
expL = exp(-(expL_arg))
farea = expL
fvol = (hL + hexp)*expL / (hi1 + hexp)
endif ! nr < ncat
! diagnostics
if (n ==1) then ! only for thinnest ridging ice
if (present(aredistn)) then
aredistn(nr) = farea*ardg2n
endif
if (present(vredistn)) then
vredistn(nr) = fvol*virdgn
endif
endif
endif ! krdg_redist
!-----------------------------------------------------------------
! Transfer ice area, ice volume, and snow volume to category nr.
!-----------------------------------------------------------------
aicen(nr) = aicen(nr) + farea*ardg2n
vicen(nr) = vicen(nr) + fvol *virdgn
vsnon(nr) = vsnon(nr) + fvol *vsrdgn*fsnowrdg
!-----------------------------------------------------------------
! Transfer area-weighted and volume-weighted tracers to category nr.
! Note: The global sum aicen*trcrn of ice area tracers
! (trcr_depend = 0) is not conserved by ridging.
! However, ridging conserves the global sum of volume
! tracers (trcr_depend = 1 or 2).
! Tracers associated with level ice, or that are otherwise lost
! from ridging ice, are not transferred.
! We assume that all pond water is lost from ridging ice.
!-----------------------------------------------------------------
do it = 1, ntrcr
if (it /= nt_alvl .and. it /= nt_vlvl) then
trfactor(1) = trcr_base(it,1)*ardg2n*farea
trfactor(2) = trcr_base(it,2)*virdgn*fvol
trfactor(3) = trcr_base(it,3)*vsrdgn*fvol*fsnowrdg
else
trfactor(1) = c0
trfactor(2) = c0
trfactor(3) = c0
endif
work = c0
do k = 1, 3
work = work + trfactor(k)*trcrn(it,n)
enddo
if (n_trcr_strata(it) > 0) then ! additional tracer layers
do itl = 1, n_trcr_strata(it)
ntr = nt_strata(it,itl)
if (ntr == nt_fbri) then ! brine fraction only
work = work * trcrn(ntr,n)
else
work = c0
endif
enddo
endif
atrcrn(it,nr) = atrcrn(it,nr) + work
enddo ! ntrcr
enddo ! nr (new ridges)
endif ! nonzero ridging
enddo ! n (ridging categories)
!-----------------------------------------------------------------
! Compute new tracers
!-----------------------------------------------------------------
do n = 1, ncat
call icepack_compute_tracers (ntrcr, trcr_depend, &
atrcrn(:,n), aicen(n), &
vicen(n), vsnon(n), &
trcr_base, n_trcr_strata, &
nt_strata, trcrn(:,n))
if (icepack_warnings_aborted(subname)) return
enddo
end subroutine ridge_shift
!=======================================================================
!autodocument_start icepack_ice_strength
! Compute the strength of the ice pack, defined as the energy (J m-2)
! dissipated per unit area removed from the ice pack under compression,
! and assumed proportional to the change in potential energy caused
! by ridging.
!
! See Rothrock (1975) and Hibler (1980).
!
! For simpler strength parameterization, see this reference:
! Hibler, W. D. III, 1979: A dynamic-thermodynamic sea ice model,
! J. Phys. Oceanog., 9, 817-846.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
subroutine icepack_ice_strength (ncat, &
aice, vice, &
aice0, aicen, &
vicen, &
strength)
integer (kind=int_kind), intent(in) :: &
ncat ! number of thickness categories
real (kind=dbl_kind), intent(in) :: &
aice , & ! concentration of ice
vice , & ! volume per unit area of ice (m)
aice0 ! concentration of open water
real (kind=dbl_kind), dimension(:), intent(in) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), intent(inout) :: &
strength ! ice strength (N/m)
!autodocument_end
! local variables
real (kind=dbl_kind) :: &
asum , & ! sum of ice and open water area
aksum ! ratio of area removed to area ridged
real (kind=dbl_kind), dimension (0:ncat) :: &
apartic ! participation function; fraction of ridging
! and closing associated w/ category n
real (kind=dbl_kind), dimension (ncat) :: &
hrmin , & ! minimum ridge thickness
hrmax , & ! maximum ridge thickness (krdg_redist = 0)
hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
krdg ! mean ridge thickness/thickness of ridging ice
integer (kind=int_kind) :: &
n ! thickness category index
real (kind=dbl_kind) :: &
hi , & ! ice thickness (m)
h2rdg , & ! mean value of h^2 for new ridge
dh2rdg ! change in mean value of h^2 per unit area
! consumed by ridging
character(len=*),parameter :: subname='(icepack_ice_strength)'
if (kstrength == 1) then ! Rothrock 1975 formulation
!-----------------------------------------------------------------
! Compute thickness distribution of ridging and ridged ice.
!-----------------------------------------------------------------
call asum_ridging (ncat, aicen, aice0, asum)
if (icepack_warnings_aborted(subname)) return
call ridge_itd (ncat, aice0, &
aicen, vicen, &
krdg_partic, krdg_redist, &
mu_rdg, &
aksum, apartic, &
hrmin, hrmax, &
hrexp, krdg)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Compute ice strength based on change in potential energy,
! as in Rothrock (1975)
!-----------------------------------------------------------------
if (krdg_redist==0) then ! Hibler 1980 formulation
do n = 1, ncat
if (aicen(n) > puny .and. apartic(n) > c0)then
hi = vicen(n) / aicen(n)
h2rdg = p333 * (hrmax(n)**3 - hrmin(n)**3) &
/ (hrmax(n) - hrmin(n))
dh2rdg = -hi*hi + h2rdg/krdg(n)
strength = strength + apartic(n) * dh2rdg
endif ! aicen > puny
enddo ! n
elseif (krdg_redist==1) then ! exponential formulation
do n = 1, ncat
if (aicen(n) > puny .and. apartic(n) > c0) then
hi = vicen(n) / aicen(n)
h2rdg = hrmin(n)*hrmin(n) &
+ c2*hrmin(n)*hrexp(n) &
+ c2*hrexp(n)*hrexp(n)
dh2rdg = -hi*hi + h2rdg/krdg(n)
strength = strength + apartic(n) * dh2rdg
endif
enddo ! n
endif ! krdg_redist
strength = Cf * Cp * strength / aksum
! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow)
! Cf accounts for frictional dissipation
else ! kstrength /= 1: Hibler (1979) form
!-----------------------------------------------------------------
! Compute ice strength as in Hibler (1979)
!-----------------------------------------------------------------
strength = Pstar*vice*exp(-Cstar*(c1-aice))
endif ! kstrength
end subroutine icepack_ice_strength
!=======================================================================
!autodocument_start icepack_step_ridge
! Computes sea ice mechanical deformation
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
subroutine icepack_step_ridge (dt, ndtd, &
nilyr, nslyr, &
nblyr, &
ncat, hin_max, &
rdg_conv, rdg_shear, &
aicen, &
trcrn, &
vicen, vsnon, &
aice0, trcr_depend, &
trcr_base, n_trcr_strata, &
nt_strata, &
dardg1dt, dardg2dt, &
dvirdgdt, opening, &
fpond, &
fresh, fhocn, &
n_aero, &
faero_ocn, fiso_ocn, &
aparticn, krdgn, &
aredistn, vredistn, &
dardg1ndt, dardg2ndt, &
dvirdgndt, &
araftn, vraftn, &
aice, fsalt, &
first_ice, fzsal, &
flux_bio, closing )
real (kind=dbl_kind), intent(in) :: &
dt ! time step
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
ndtd , & ! number of dynamics supercycles
nblyr , & ! number of bio layers
nilyr , & ! number of ice layers
nslyr , & ! number of snow layers
n_aero ! number of aerosol tracers
real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max ! category limits (m)
integer (kind=int_kind), dimension (:), intent(in) :: &
trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
n_trcr_strata ! number of underlying tracer layers
real (kind=dbl_kind), dimension (:,:), intent(in) :: &
trcr_base ! = 0 or 1 depending on tracer dependency
! argument 2: (1) aice, (2) vice, (3) vsno
integer (kind=int_kind), dimension (:,:), intent(in) :: &
nt_strata ! indices of underlying tracer layers
real (kind=dbl_kind), intent(inout) :: &
aice , & ! sea ice concentration
aice0 , & ! concentration of open water
rdg_conv , & ! convergence term for ridging (1/s)
rdg_shear, & ! shear term for ridging (1/s)
dardg1dt , & ! rate of area loss by ridging ice (1/s)
dardg2dt , & ! rate of area gain by new ridges (1/s)
dvirdgdt , & ! rate of ice volume ridged (m/s)
opening , & ! rate of opening due to divergence/shear (1/s)
fpond , & ! fresh water flux to ponds (kg/m^2/s)
fresh , & ! fresh water flux to ocean (kg/m^2/s)
fsalt , & ! salt flux to ocean (kg/m^2/s)
fhocn , & ! net heat flux to ocean (W/m^2)
fzsal ! zsalinity flux to ocean(kg/m^2/s)
real (kind=dbl_kind), intent(inout), optional :: &
closing ! rate of closing due to divergence/shear (1/s)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon , & ! volume per unit area of snow (m)
dardg1ndt, & ! rate of area loss by ridging ice (1/s)
dardg2ndt, & ! rate of area gain by new ridges (1/s)
dvirdgndt, & ! rate of ice volume ridged (m/s)
aparticn , & ! participation function
krdgn , & ! mean ridge thickness/thickness of ridging ice
araftn , & ! rafting ice area
vraftn , & ! rafting ice volume
aredistn , & ! redistribution function: fraction of new ridge area
vredistn , & ! redistribution function: fraction of new ridge volume
faero_ocn, & ! aerosol flux to ocean (kg/m^2/s)
flux_bio ! all bio fluxes to ocean
real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
trcrn ! tracers
!logical (kind=log_kind), intent(in) :: &
!tr_pond_topo,& ! if .true., use explicit topography-based ponds
!tr_aero ,& ! if .true., use aerosol tracers
!tr_brine !,& ! if .true., brine height differs from ice thickness
!heat_capacity ! if true, ice has nonzero heat capacity
logical (kind=log_kind), dimension(:), intent(inout) :: &
first_ice ! true until ice forms
!autodocument_end
! local variables
real (kind=dbl_kind) :: &
dtt ! thermo time step
real (kind=dbl_kind), dimension(:), allocatable :: &
l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
real (kind=dbl_kind) :: &
l_closing ! local rate of closing due to divergence/shear (1/s)
logical (kind=log_kind) :: &
l_closing_flag ! flag if closing is passed
character(len=*),parameter :: subname='(icepack_step_ridge)'
!-----------------------------------------------------------------
! Identify ice-ocean cells.
! Note: We can not limit the loop here using aice>puny because
! aice has not yet been updated since the transport (and
! it may be out of whack, which the ridging helps fix).-ECH
!-----------------------------------------------------------------
if (present(fiso_ocn)) then
allocate(l_fiso_ocn(size(fiso_ocn)))
l_fiso_ocn = fiso_ocn
else
! check tr_iso = true ???
allocate(l_fiso_ocn(1))
l_fiso_ocn = c0
endif
if (present(closing)) then
l_closing_flag = .true.
l_closing = closing
else
l_closing_flag = .false.
l_closing = c0
endif
call ridge_ice (dt, ndtd, &
ncat, n_aero, &
nilyr, nslyr, &
ntrcr, hin_max, &
rdg_conv, rdg_shear, &
aicen, &
trcrn, &
vicen, vsnon, &
aice0, &
trcr_depend, &
trcr_base, &
n_trcr_strata, &
nt_strata, &
krdg_partic, krdg_redist, &
mu_rdg, tr_brine, &
dardg1dt, dardg2dt, &
dvirdgdt, opening, &
fpond, &
fresh, fhocn, &
faero_ocn, l_fiso_ocn, &
aparticn, krdgn, &
aredistn, vredistn, &
dardg1ndt, dardg2ndt, &
dvirdgndt, &
araftn, vraftn, &
l_closing_flag, &
l_closing )
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! ITD cleanup: Rebin thickness categories if necessary, and remove
! categories with very small areas.
!-----------------------------------------------------------------
dtt = dt * ndtd ! for proper averaging over thermo timestep
call cleanup_itd (dtt, ntrcr, &
nilyr, nslyr, &
ncat, hin_max, &
aicen, trcrn, &
vicen, vsnon, &
aice0, aice, &
n_aero, &
nbtrcr, nblyr, &
tr_aero, &
tr_pond_topo, heat_capacity, &
first_ice, &
trcr_depend, trcr_base, &
n_trcr_strata, nt_strata, &
fpond, fresh, &
fsalt, fhocn, &
faero_ocn, l_fiso_ocn, &
fzsal, &
flux_bio)
if (icepack_warnings_aborted(subname)) return
if (present(fiso_ocn)) fiso_ocn = l_fiso_ocn
deallocate(l_fiso_ocn)
end subroutine icepack_step_ridge
!=======================================================================
end module icepack_mechred
!=======================================================================