!=========================================================================
!
! This module contains the subroutines required to define
! a floe size distribution tracer for sea ice
!
! Theory based on:
!
! Horvat, C., & Tziperman, E. (2015). A prognostic model of the sea-ice
! floe size and thickness distribution. The Cryosphere, 9(6), 2119-2134.
! doi:10.5194/tc-9-2119-2015
!
! and implementation described in:
!
! Roach, L. A., Horvat, C., Dean, S. M., & Bitz, C. M. (2018). An emergent
! sea ice floe size distribution in a global coupled ocean--sea ice model.
! Journal of Geophysical Research: Oceans, 123(6), 4322-4337.
! doi:10.1029/2017JC013692
!
! with some modifications.
!
! For floe welding parameter and tensile mode parameter values, see
!
! Roach, L. A., Smith, M. M., & Dean, S. M. (2018). Quantifying
! growth of pancake sea ice floes using images from drifting buoys.
! Journal of Geophysical Research: Oceans, 123(4), 2851-2866.
! doi: 10.1002/2017JC013693
!
! Variable naming convention:
! for k = 1, nfsd and n = 1, ncat
! afsdn(k,n) = trcrn(:,:,nt_nfsd+k-1,n,:)
! afsd (k) is per thickness category or averaged over n
! afs is the associated scalar value for (k,n)
!
! authors: Lettie Roach, VUW/NIWA
! C. M. Bitz, UW
!
! 2016: CMB started
! 2016-8: LR worked on most of it
! 2019: ECH ported to Icepack
!-----------------------------------------------------------------
module icepack_fsd
use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c4, p01, p1, p5, puny
use icepack_parameters, only: pi, floeshape, wave_spec, bignum, gravit, rhoi
use icepack_tracers, only: nt_fsd, tr_fsd
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
implicit none
private
public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, &
fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd
real(kind=dbl_kind), dimension(:), allocatable :: &
floe_rad_h, & ! fsd size higher bound in m (radius)
floe_area_l, & ! fsd area at lower bound (m^2)
floe_area_h, & ! fsd area at higher bound (m^2)
floe_area_c, & ! fsd area at bin centre (m^2)
floe_area_binwidth ! floe area bin width (m^2)
integer(kind=int_kind), dimension(:,:), allocatable, public :: &
iweld ! floe size categories that can combine
! during welding (dimensionless)
!=======================================================================
contains
!=======================================================================
! Note that radius widths cannot be larger than twice previous
! category width or floe welding will not have an effect
!
! Note also that the bound of the lowest floe size category is used
! to define the lead region width and the domain spacing for wave fracture
!
!autodocument_start icepack_init_fsd_bounds
! Initialize ice fsd bounds (call whether or not restarting)
! Define the bounds, midpoints and widths of floe size
! categories in area and radius
!
! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
subroutine icepack_init_fsd_bounds(nfsd, &
floe_rad_l, & ! fsd size lower bound in m (radius)
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth, & ! fsd size bin width in m (radius)
c_fsd_range, & ! string for history output
write_diags ) ! flag for writing diagnostics
integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories
real(kind=dbl_kind), dimension(:), intent(inout) :: &
floe_rad_l, & ! fsd size lower bound in m (radius)
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
character (len=35), intent(out) :: &
c_fsd_range(nfsd) ! string for history output
logical (kind=log_kind), intent(in), optional :: &
write_diags ! write diags flag
!autodocument_end
! local variables
integer (kind=int_kind) :: n, m, k
integer (kind=int_kind) :: ierr
real (kind=dbl_kind) :: test
real (kind=dbl_kind), dimension (nfsd+1) :: &
area_lims, area_lims_scaled
real (kind=dbl_kind), dimension (0:nfsd) :: &
floe_rad
real (kind=dbl_kind), dimension(:), allocatable :: &
lims
logical (kind=log_kind) :: &
l_write_diags ! local write diags
character(len=8) :: c_fsd1,c_fsd2
character(len=2) :: c_nf
character(len=*), parameter :: subname='(icepack_init_fsd_bounds)'
l_write_diags = .true.
if (present(write_diags)) then
l_write_diags = write_diags
endif
if (nfsd.eq.24) then
allocate(lims(24+1))
lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, &
3.35434988e+03, 4.55051413e+03, 6.17323164e+03, 8.37461170e+03, &
1.13610059e+04, 1.54123510e+04, 2.09084095e+04, 2.83643675e+04, &
3.84791270e+04 /)
elseif (nfsd.eq.16) then
allocate(lims(16+1))
lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, &
3.35434988e+03 /)
else if (nfsd.eq.12) then
allocate(lims(12+1))
lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, &
5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, &
3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, &
9.45812834e+02 /)
else if (nfsd.eq.1) then ! default case
allocate(lims(1+1))
lims = (/ 6.65000000e-02, 3.0e+02 /)
else
call icepack_warnings_add(subname//&
' floe size categories not defined for nfsd')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
end if
allocate( &
floe_rad_h (nfsd), & ! fsd size higher bound in m (radius)
floe_area_l (nfsd), & ! fsd area at lower bound (m^2)
floe_area_h (nfsd), & ! fsd area at higher bound (m^2)
floe_area_c (nfsd), & ! fsd area at bin centre (m^2)
floe_area_binwidth (nfsd), & ! floe area bin width (m^2)
iweld (nfsd, nfsd), & ! fsd categories that can weld
stat=ierr)
if (ierr/=0) then
call icepack_warnings_add(subname//' Out of Memory fsd')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
floe_rad_l = lims(1:nfsd )
floe_rad_h = lims(2:nfsd+1)
floe_rad_c = (floe_rad_h+floe_rad_l)/c2
floe_area_l = c4*floeshape*floe_rad_l**2
floe_area_c = c4*floeshape*floe_rad_c**2
floe_area_h = c4*floeshape*floe_rad_h**2
floe_binwidth = floe_rad_h - floe_rad_l
floe_area_binwidth = floe_area_h - floe_area_l
! floe size categories that can combine during welding
iweld(:,:) = -999
do n = 1, nfsd
do m = 1, nfsd
test = floe_area_c(n) + floe_area_c(m)
do k = 1, nfsd-1
if ((test >= floe_area_l(k)) .and. (test < floe_area_h(k))) &
iweld(n,m) = k
end do
if (test >= floe_area_l(nfsd)) iweld(n,m) = nfsd
end do
end do
if (allocated(lims)) deallocate(lims)
! write fsd bounds
floe_rad(0) = floe_rad_l(1)
do n = 1, nfsd
floe_rad(n) = floe_rad_h(n)
! Save character string to write to history file
write (c_nf, '(i2)') n
write (c_fsd1, '(f6.3)') floe_rad(n-1)
write (c_fsd2, '(f6.3)') floe_rad(n)
c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m'
enddo
if (l_write_diags) then
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname
call icepack_warnings_add(warnstr)
write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
call icepack_warnings_add(warnstr)
do n = 1, nfsd
write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
call icepack_warnings_add(warnstr)
enddo
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
endif
end subroutine icepack_init_fsd_bounds
!=======================================================================
! When growing from no-ice conditions, initialize to zero.
! This allows the FSD to emerge, as described in Roach, Horvat et al. (2018)
!
! Otherwise initalize with a power law, following Perovich
! & Jones (2014). The basin-wide applicability of such a
! prescribed power law has not yet been tested.
!
! Perovich, D. K., & Jones, K. F. (2014). The seasonal evolution of
! sea ice floe size distribution. Journal of Geophysical Research: Oceans,
! 119(12), 8767–8777. doi:10.1002/2014JC010136
!
!autodocument_start icepack_init_fsd
!
! Initialize the FSD
!
! authors: Lettie Roach, NIWA/VUW
subroutine icepack_init_fsd(nfsd, ice_ic, &
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth, & ! fsd size bin width in m (radius)
afsd) ! floe size distribution tracer
integer(kind=int_kind), intent(in) :: &
nfsd
character(len=char_len_long), intent(in) :: &
ice_ic ! method of ice cover initialization
real(kind=dbl_kind), dimension(:), intent(inout) :: &
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
afsd ! floe size tracer: fraction distribution of floes
!autodocument_end
! local variables
real (kind=dbl_kind) :: alpha, totfrac
integer (kind=int_kind) :: k
real (kind=dbl_kind), dimension (nfsd) :: &
num_fsd ! number distribution of floes
if (trim(ice_ic) == 'none') then
afsd(:) = c0
else ! Perovich (2014)
! fraction of ice in each floe size and thickness category
! same for ALL cells (even where no ice) initially
alpha = 2.1_dbl_kind
totfrac = c0 ! total fraction of floes
do k = 1, nfsd
num_fsd(k) = (2*floe_rad_c(k))**(-alpha-c1) ! number distribution of floes
afsd (k) = num_fsd(k)*floe_area_c(k)*floe_binwidth(k) ! fraction distribution of floes
totfrac = totfrac + afsd(k)
enddo
afsd = afsd/totfrac ! normalize
endif ! ice_ic
end subroutine icepack_init_fsd
!=======================================================================
!autodocument_start icepack_cleanup_fsd
!
! Clean up small values and renormalize
!
! authors: Elizabeth Hunke, LANL
!
subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd ! number of floe size categories
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
afsdn ! floe size distribution tracer
!autodocument_end
! local variables
integer (kind=int_kind) :: &
n ! thickness category index
character(len=*), parameter :: subname='(icepack_cleanup_fsd)'
if (tr_fsd) then
do n = 1, ncat
call icepack_cleanup_fsdn(nfsd, afsdn(:,n))
if (icepack_warnings_aborted(subname)) return
enddo
endif ! tr_fsd
end subroutine icepack_cleanup_fsd
!=======================================================================
!
! Clean up small values and renormalize -- per category
!
! authors: Elizabeth Hunke, LANL
!
subroutine icepack_cleanup_fsdn (nfsd, afsd)
integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories
real (kind=dbl_kind), dimension(:), intent(inout) :: &
afsd ! floe size distribution tracer
! local variables
integer (kind=int_kind) :: &
k ! floe size category index
real (kind=dbl_kind) :: &
tot
do k = 1, nfsd
if (afsd(k) < puny) afsd(k) = c0
enddo
tot = sum(afsd(:))
if (tot > puny) then
do k = 1, nfsd
afsd(k) = afsd(k) / tot ! normalize
enddo
else ! represents ice-free cell, so set to zero
afsd(:) = c0
endif
end subroutine icepack_cleanup_fsdn
!=======================================================================
!
! Given the joint ice thickness and floe size distribution, calculate
! the lead region and the total lateral surface area following Horvat
! and Tziperman (2015).
!
! authors: Lettie Roach, NIWA/VUW
subroutine partition_area (ncat, nfsd, &
floe_rad_c, aice, &
aicen, vicen, &
afsdn, lead_area, &
latsurf_area)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd ! number of floe size categories
real (kind=dbl_kind), dimension(:), intent(in) :: &
floe_rad_c ! fsd size bin centre in m (radius)
real (kind=dbl_kind), intent(in) :: &
aice ! ice concentration
real (kind=dbl_kind), dimension(:), intent(in) :: &
aicen , & ! fractional area of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
afsdn ! floe size distribution tracer
real (kind=dbl_kind), intent(out) :: &
lead_area , & ! the fractional area of the lead region
latsurf_area ! the area covered by lateral surface of floes
! local variables
integer (kind=int_kind) :: &
n , & ! thickness category index
k ! floe size index
real (kind=dbl_kind) :: &
width_leadreg, & ! width of lead region (m)
thickness ! actual thickness of ice in thickness cat (m)
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
lead_area = c0
latsurf_area = c0
! Set the width of the lead region to be the smallest
! floe size category, as per Horvat & Tziperman (2015)
width_leadreg = floe_rad_c(1)
! Only calculate these areas where there is some ice
if (aice > puny) then
! lead area = sum of areas of annuli around all floes
do n = 1, ncat
do k = 1, nfsd
lead_area = lead_area + aicen(n) * afsdn(k,n) &
* (c2*width_leadreg /floe_rad_c(k) &
+ width_leadreg**2/floe_rad_c(k)**2)
enddo ! k
enddo ! n
! cannot be greater than the open water fraction
lead_area=MIN(lead_area, c1-aice)
if (lead_area < c0) then
if (lead_area < -puny) then
! stop 'lead_area lt0 in partition_area'
else
lead_area=MAX(lead_area,c0)
end if
end if
! Total fractional lateral surface area in each grid (per unit ocean area)
do n = 1, ncat
thickness = c0
if (aicen(n) > c0) thickness = vicen(n)/aicen(n)
do k = 1, nfsd
latsurf_area = latsurf_area &
+ afsdn(k,n) * aicen(n) * c2 * thickness/floe_rad_c(k)
end do
end do
! check
! if (latsurf_area < c0) stop 'negative latsurf_area'
! if (latsurf_area /= latsurf_area) stop 'latsurf_area NaN'
end if ! aice
end subroutine partition_area
!=======================================================================
!
! Lateral growth at the edges of floes
!
! Compute the portion of new ice growth that occurs at the edges of
! floes. The remainder will grow as new ice frazil ice in open water
! (away from floes).
!
! See Horvat & Tziperman (2015) and Roach, Horvat et al. (2018).
!
! authors: Lettie Roach, NIWA/VUW
!
subroutine fsd_lateral_growth (ncat, nfsd, &
dt, aice, &
aicen, vicen, &
vi0new, &
frazil, floe_rad_c, &
afsdn, &
lead_area, latsurf_area, &
G_radial, d_an_latg, &
tot_latg)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd ! number of floe size categories
real (kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
aice ! total concentration of ice
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
afsdn ! floe size distribution tracer
real (kind=dbl_kind), intent(inout) :: &
vi0new , & ! volume of new ice added to cat 1 (m)
frazil ! frazil ice growth (m/step-->cm/day)
! floe size distribution
real (kind=dbl_kind), dimension (:), intent(in) :: &
floe_rad_c ! fsd size bin centre in m (radius)
real (kind=dbl_kind), dimension(ncat), intent(out) :: &
d_an_latg ! change in aicen occuring due to lateral growth
real (kind=dbl_kind), intent(out) :: &
G_radial , & ! lateral melt rate (m/s)
tot_latg ! total change in aice due to
! lateral growth at the edges of floes
! local variables
integer (kind=int_kind) :: &
n, k ! ice category indices
real (kind=dbl_kind) :: &
vi0new_lat ! volume of new ice added laterally to FSD (m)
real (kind=dbl_kind), intent(out) :: &
lead_area , & ! the fractional area of the lead region
latsurf_area ! the area covered by lateral surface of floes
character(len=*),parameter :: subname='(fsd_lateral_growth)'
lead_area = c0
latsurf_area = c0
G_radial = c0
tot_latg = c0
d_an_latg = c0
! partition volume into lateral growth and frazil
call partition_area (ncat, nfsd, &
floe_rad_c, aice, &
aicen, vicen, &
afsdn, lead_area, &
latsurf_area)
if (icepack_warnings_aborted(subname)) return
vi0new_lat = c0
if (latsurf_area > puny) then
vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area)
end if
! for history/diagnostics
frazil = vi0new - vi0new_lat
! lateral growth increment
if (vi0new_lat > puny) then
G_radial = vi0new_lat/dt
do n = 1, ncat
do k = 1, nfsd
d_an_latg(n) = d_an_latg(n) &
+ c2*aicen(n)*afsdn(k,n)*G_radial*dt/floe_rad_c(k)
end do
end do ! n
! cannot expand ice laterally beyond lead region
if (SUM(d_an_latg(:)).ge.lead_area) then
d_an_latg(:) = d_an_latg(:)/SUM(d_an_latg(:))
d_an_latg(:) = d_an_latg(:)*lead_area
end if
endif ! vi0new_lat > 0
! Use remaining ice volume as in standard model,
! but ice cannot grow into the area that has grown laterally
vi0new = vi0new - vi0new_lat
tot_latg = SUM(d_an_latg(:))
end subroutine fsd_lateral_growth
!=======================================================================
!
! Evolve the FSD subject to lateral growth and the growth of new ice
! in thickness category 1.
!
! If ocean surface wave forcing is provided, the floe size of new ice
! (grown away from floe edges) can computed from the wave field
! based on the tensile (stretching) stress limitation following
! Shen et al. (2001). Otherwise, new floes all grow in the smallest
! floe size category, representing pancake ice formation.
!
! Shen, H., Ackley, S., & Hopkins, M. (2001). A conceptual model
! for pancake-ice formation in a wave field.
! Annals of Glaciology, 33, 361-367. doi:10.3189/172756401781818239
!
! authors: Lettie Roach, NIWA/VUW
!
subroutine fsd_add_new_ice (ncat, n, nfsd, &
dt, ai0new, &
d_an_latg, d_an_newi, &
floe_rad_c, floe_binwidth, &
G_radial, area2, &
wave_sig_ht, &
wave_spectrum, &
wavefreq, &
dwavefreq, &
d_afsd_latg, &
d_afsd_newi, &
afsdn, aicen_init, &
aicen, trcrn)
integer (kind=int_kind), intent(in) :: &
n , & ! thickness category number
ncat , & ! number of thickness categories
nfsd ! number of floe size categories
real (kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
ai0new , & ! area of new ice added to cat 1
G_radial , & ! lateral melt rate (m/s)
wave_sig_ht ! wave significant height (everywhere) (m)
real (kind=dbl_kind), dimension(:), intent(in) :: &
wave_spectrum ! ocean surface wave spectrum as a function of frequency
! power spectral density of surface elevation, E(f) (units m^2 s)
real(kind=dbl_kind), dimension(:), intent(in) :: &
wavefreq , & ! wave frequencies (s^-1)
dwavefreq ! wave frequency bin widths (s^-1)
real (kind=dbl_kind), dimension(:), intent(in) :: &
d_an_latg , & ! change in aicen due to lateral growth
d_an_newi ! change in aicen due to frazil ice formation
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen_init , & ! fractional area of ice
aicen , & ! after update
floe_rad_c , & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
real (kind=dbl_kind), dimension (:,:), intent(in) :: &
afsdn ! floe size distribution tracer
real (kind=dbl_kind), dimension (:), intent(in) :: &
area2 ! area after lateral growth, before new ice formation
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn ! ice tracers
real (kind=dbl_kind), dimension(:), intent(inout) :: &
! change in floe size distribution (area)
d_afsd_latg, & ! due to fsd lateral growth
d_afsd_newi ! new ice formation
integer (kind=int_kind) :: &
k, & ! floe size category index
new_size, & ! index for floe size of new ice
nsubt ! number of subtimesteps
real (kind=dbl_kind) :: &
elapsed_t, subdt ! elapsed time, subtimestep (s)
real (kind=dbl_kind), dimension (nfsd,ncat) :: &
afsdn_latg ! fsd after lateral growth
real (kind=dbl_kind), dimension (nfsd) :: &
dafsd_tmp, & ! tmp FSD
df_flx , & ! finite differences for fsd
afsd_ni ! fsd after new ice added
real (kind=dbl_kind), dimension(nfsd+1) :: &
f_flx ! finite differences in floe size
character(len=*),parameter :: subname='(fsd_add_new_ice)'
afsdn_latg(:,n) = afsdn(:,n) ! default
if (d_an_latg(n) > puny) then ! lateral growth
! adaptive timestep
elapsed_t = c0
nsubt = 0
DO WHILE (elapsed_t.lt.dt)
nsubt = nsubt + 1
if (nsubt.gt.100) print *, 'latg not converging'
! finite differences
df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radial * afsdn_latg(k-1,n) / floe_binwidth(k-1)
end do
do k = 1, nfsd
df_flx(k) = f_flx(k+1) - f_flx(k)
end do
! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0'
dafsd_tmp(:) = c0
do k = 1, nfsd
dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) &
* (c1/floe_rad_c(k) - SUM(afsdn_latg(:,n)/floe_rad_c(:))) )
end do
! timestep required for this
subdt = get_subdt_fsd(nfsd, afsdn_latg(:,n), dafsd_tmp(:))
subdt = MIN(subdt, dt)
! update fsd and elapsed time
afsdn_latg(:,n) = afsdn_latg(:,n) + subdt*dafsd_tmp(:)
elapsed_t = elapsed_t + subdt
END DO
call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n))
if (icepack_warnings_aborted(subname)) return
trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n)
end if ! lat growth
new_size = nfsd
if (n == 1) then
! add new frazil ice to smallest thickness
if (d_an_newi(n) > puny) then
afsd_ni(:) = c0
if (SUM(afsdn_latg(:,n)) > puny) then ! fsd exists
if (wave_spec) then
if (wave_sig_ht > puny) then
call wave_dep_growth (nfsd, wave_spectrum, &
wavefreq, dwavefreq, &
new_size)
if (icepack_warnings_aborted(subname)) return
end if
! grow in new_size category
afsd_ni(new_size) = (afsdn_latg(new_size,n)*area2(n) + ai0new) &
/ (area2(n) + ai0new)
do k = 1, new_size-1 ! diminish other floe cats accordingly
afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new)
end do
do k = new_size+1, nfsd ! diminish other floe cats accordingly
afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new)
end do
else ! grow in smallest floe size category
afsd_ni(1) = (afsdn_latg(1,n)*area2(n) + ai0new) &
/ (area2(n) + ai0new)
do k = 2, nfsd ! diminish other floe cats accordingly
afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n)+ai0new)
enddo
end if ! wave spec
else ! no fsd, so entirely new ice
if (wave_spec) then
if (wave_sig_ht > puny) then
call wave_dep_growth (nfsd, wave_spectrum, &
wavefreq, dwavefreq, &
new_size)
if (icepack_warnings_aborted(subname)) return
end if
afsd_ni(new_size) = c1
else
afsd_ni(1) = c1
endif ! wave forcing
endif ! entirely new ice or not
trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_ni(:)
call icepack_cleanup_fsdn (nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,n))
if (icepack_warnings_aborted(subname)) return
endif ! d_an_newi > puny
endif ! n = 1
! history/diagnostics
do k = 1, nfsd
! sum over n
d_afsd_latg(k) = d_afsd_latg(k) &
+ area2(n)*afsdn_latg(k,n) & ! after latg
- aicen_init(n)*afsdn(k,n) ! at start
d_afsd_newi(k) = d_afsd_newi(k) &
+ aicen(n)*trcrn(nt_fsd+k-1,n) & ! after latg and newi
- area2(n)*afsdn_latg(k,n) ! after latg
enddo ! k
end subroutine fsd_add_new_ice
!=======================================================================
!
! Given a wave spectrum, calculate size of new floes based on
! tensile failure, following Shen et al. (2001)
!
! The tensile mode parameter is based on in-situ measurements
! by Roach, Smith & Dean (2018).
!
! authors: Lettie Roach, NIWA/VUW
!
subroutine wave_dep_growth (nfsd, local_wave_spec, &
wavefreq, dwavefreq, &
new_size)
integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories
real (kind=dbl_kind), dimension(:), intent(in) :: &
local_wave_spec ! ocean surface wave spectrum as a function of frequency
! power spectral density of surface elevation, E(f) (units m^2 s)
! dimension set in ice_forcing
real(kind=dbl_kind), dimension(:), intent(in) :: &
wavefreq, & ! wave frequencies (s^-1)
dwavefreq ! wave frequency bin widths (s^-1)
integer (kind=int_kind), intent(out) :: &
new_size ! index of floe size category in which new floes will grow
! local variables
real (kind=dbl_kind), parameter :: &
tensile_param = 0.167_dbl_kind ! tensile mode parameter (kg m^-1 s^-2)
! value from Roach, Smith & Dean (2018)
real (kind=dbl_kind) :: &
w_amp, & ! wave amplitude (m)
f_peak, & ! peak frequency (s^-1)
r_max ! floe radius (m)
integer (kind=int_kind) :: k
w_amp = c2* SQRT(SUM(local_wave_spec*dwavefreq)) ! sig wave amplitude
f_peak = wavefreq(MAXLOC(local_wave_spec, DIM=1)) ! peak frequency
! tensile failure
if (w_amp > puny .and. f_peak > puny) then
r_max = p5*SQRT(tensile_param*gravit/(pi**5*rhoi*w_amp*2))/f_peak**2
else
r_max = bignum
end if
new_size = nfsd
do k = nfsd-1, 1, -1
if (r_max <= floe_rad_h(k)) new_size = k
end do
end subroutine wave_dep_growth
!=======================================================================
!
! Floes are perimitted to weld together in freezing conditions, according
! to their geometric probability of overlap if placed randomly on the
! domain. The rate per unit area c_weld is the total number
! of floes that weld with another, per square meter, per unit time, in the
! case of a fully covered ice surface (aice=1), equal to twice the reduction
! in total floe number. See Roach, Smith & Dean (2018).
!
!
! authors: Lettie Roach, NIWA/VUW
!
subroutine fsd_weld_thermo (ncat, nfsd, &
dt, frzmlt, &
aicen, trcrn, &
d_afsd_weld)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd ! number of floe size categories
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen ! ice concentration
real (kind=dbl_kind), intent(in) :: &
frzmlt ! freezing/melting potential (W/m^2)
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn ! ice tracers
real (kind=dbl_kind), dimension (:), intent(inout) :: &
d_afsd_weld ! change in fsd due to welding
! local variables
real (kind=dbl_kind), parameter :: &
aminweld = p1 ! minimum ice concentration likely to weld
real (kind=dbl_kind), parameter :: &
c_weld = 1.0e-8_dbl_kind
! constant of proportionality for welding
! total number of floes that weld with another, per square meter,
! per unit time, in the case of a fully covered ice surface
! units m^-2 s^-1, see documentation for details
integer (kind=int_kind) :: &
nt , & ! time step index
n , & ! thickness category index
k, kx, ky, i, j ! floe size category indices
real (kind=dbl_kind), dimension(nfsd,ncat) :: &
afsdn ! floe size distribution tracer
real (kind=dbl_kind), dimension(nfsd,ncat) :: &
d_afsdn_weld ! change in afsdn due to welding
real (kind=dbl_kind), dimension(nfsd) :: &
stability , & ! check for stability
nfsd_tmp , & ! number fsd
afsd_init , & ! initial values
afsd_tmp , & ! work array
gain, loss ! welding tendencies
real(kind=dbl_kind) :: &
prefac , & ! multiplies kernel
kern , & ! kernel
subdt , & ! subcycling time step for stability (s)
elapsed_t ! elapsed subcycling time
character(len=*), parameter :: subname='(fsd_weld_thermo)'
afsdn (:,:) = c0
afsd_init(:) = c0
stability = c0
prefac = p5
do n = 1, ncat
d_afsd_weld (:) = c0
d_afsdn_weld(:,n) = c0
afsdn(:,n) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
if (icepack_warnings_aborted(subname)) return
! If there is some ice in the lower (nfsd-1) categories
! and there is freezing potential
if ((frzmlt > puny) .and. & ! freezing potential
(aicen(n) > aminweld) .and. & ! low concentrations unlikely to weld
(SUM(afsdn(1:nfsd-1,n)) > puny)) then ! some ice in nfsd-1 categories
afsd_init(:) = afsdn(:,n) ! save initial values
afsd_tmp (:) = afsd_init(:) ! work array
! in case of minor numerical errors
WHERE(afsd_tmp < puny) afsd_tmp = c0
afsd_tmp = afsd_tmp/SUM(afsd_tmp)
! adaptive sub-timestep
elapsed_t = c0
DO WHILE (elapsed_t < dt)
! calculate sub timestep
nfsd_tmp = afsd_tmp/floe_area_c
WHERE (afsd_tmp > puny) &
stability = nfsd_tmp/(c_weld*afsd_tmp*aicen(n))
WHERE (stability < puny) stability = bignum
subdt = MINVAL(stability)
subdt = MIN(subdt,dt)
loss(:) = c0
gain(:) = c0
do i = 1, nfsd ! consider loss from this category
do j = 1, nfsd ! consider all interaction partners
k = iweld(i,j) ! product of i+j
if (k > i) then
kern = c_weld * floe_area_c(i) * aicen(n)
loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j)
if (i.eq.j) prefac = c1 ! otherwise 0.5
gain(k) = gain(k) + prefac*kern*afsd_tmp(i)*afsd_tmp(j)
end if
end do
end do
! does largest category lose?
! if (loss(nfsd) > puny) stop 'weld, largest cat losing'
! if (gain(1) > puny) stop 'weld, smallest cat gaining'
! update afsd
afsd_tmp(:) = afsd_tmp(:) + subdt*(gain(:) - loss(:))
! in case of minor numerical errors
WHERE(afsd_tmp < puny) afsd_tmp = c0
afsd_tmp = afsd_tmp/SUM(afsd_tmp)
! update time
elapsed_t = elapsed_t + subdt
! stop if all in largest floe size cat
if (afsd_tmp(nfsd) > (c1-puny)) exit
END DO ! time
call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
if (icepack_warnings_aborted(subname)) return
do k = 1, nfsd
afsdn(k,n) = afsd_tmp(k)
trcrn(nt_fsd+k-1,n) = afsdn(k,n)
! history/diagnostics
d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k)
enddo
endif ! try to weld
enddo ! ncat
! history/diagnostics
do k = 1, nfsd
d_afsd_weld(k) = c0
do n = 1, ncat
d_afsd_weld(k) = d_afsd_weld(k) + aicen(n)*d_afsdn_weld(k,n)
end do ! n
end do ! k
end subroutine fsd_weld_thermo
!=======================================================================
!
! Adaptive timestepping (process agnostic)
! See reference: Horvat & Tziperman (2017) JGR, Appendix A
!
! authors: 2018 Lettie Roach, NIWA/VUW
!
!
function get_subdt_fsd(nfsd, afsd_init, d_afsd) &
result(subdt)
integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories
real (kind=dbl_kind), dimension (nfsd), intent(in) :: &
afsd_init, d_afsd ! floe size distribution tracer
! output
real (kind=dbl_kind) :: &
subdt ! subcycle timestep (s)
! local variables
real (kind=dbl_kind), dimension (nfsd) :: &
check_dt ! to compute subcycle timestep (s)
integer (kind=int_kind) :: k
check_dt(:) = bignum
do k = 1, nfsd
if (d_afsd(k) > puny) check_dt(k) = (1-afsd_init(k))/d_afsd(k)
if (d_afsd(k) < -puny) check_dt(k) = afsd_init(k)/ABS(d_afsd(k))
end do
subdt = MINVAL(check_dt)
end function get_subdt_fsd
!=======================================================================
end module icepack_fsd
!=======================================================================