!=======================================================================
!
! Thermo calculations after call to coupler, related to ITD:
! ice thickness redistribution, lateral growth and melting.
!
! NOTE: The thermodynamic calculation is split in two for load balancing.
! First icepack_therm_vertical computes vertical growth rates and coupler
! fluxes. Then icepack_therm_itd does thermodynamic calculations not
! needed for coupling.
!
! authors William H. Lipscomb, LANL
! C. M. Bitz, UW
! Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
! 2004: Block structure added by William Lipscomb
! 2006: Streamlined for efficiency by Elizabeth Hunke
! 2014: Column package created by Elizabeth Hunke
!
module icepack_therm_itd
use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10
use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
use icepack_parameters, only: rhosi, conserv_check
use icepack_parameters, only: kitd, ktherm, heat_capacity
use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin
use icepack_tracers, only: ntrcr, nbtrcr
use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd
use icepack_tracers, only: nt_alvl, nt_vlvl
use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
use icepack_tracers, only: n_aero, n_iso
use icepack_tracers, only: bio_index
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd, get_subdt_fsd
use icepack_itd, only: reduce_area, cleanup_itd
use icepack_itd, only: aggregate_area, shift_ice
use icepack_itd, only: column_sum, column_conservation_check
use icepack_isotope, only: isoice_alpha, isotope_frac_method
use icepack_mushy_physics, only: liquidus_temperature_mush, enthalpy_mush
use icepack_therm_shared, only: hi_min
use icepack_zbgc, only: add_new_ice_bgc
use icepack_zbgc, only: lateral_melt_bgc
implicit none
private
public :: linear_itd, &
add_new_ice, &
lateral_melt, &
icepack_step_therm2
!=======================================================================
contains
!=======================================================================
!
! ITD scheme that shifts ice among categories
!
! See Lipscomb, W. H. Remapping the thickness distribution in sea
! ice models. 2001, J. Geophys. Res., Vol 106, 13989--14000.
!
! Using the thermodynamic "velocities", interpolate to find the
! velocities in thickness space at the category boundaries, and
! compute the new locations of the boundaries. Then for each
! category, compute the thickness distribution function, g(h),
! between hL and hR, the left and right boundaries of the category.
! Assume g(h) is a linear polynomial that satisfies two conditions:
!
! (1) The ice area implied by g(h) equals aicen(n).
! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
!
! Given g(h), at each boundary compute the ice area and volume lying
! between the original and new boundary locations. Transfer area
! and volume across each boundary in the appropriate direction, thus
! restoring the original boundaries.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
subroutine linear_itd (ncat, hin_max, &
nilyr, nslyr, &
ntrcr, trcr_depend, &
trcr_base, n_trcr_strata,&
nt_strata, &
aicen_init, vicen_init, &
aicen, trcrn, &
vicen, vsnon, &
aice, aice0, &
fpond )
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nilyr , & ! number of ice layers
nslyr , & ! number of snow layers
ntrcr ! number of tracers in use
real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max ! category boundaries (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), dimension(:), intent(in) :: &
aicen_init, & ! initial ice concentration (before vertical thermo)
vicen_init ! initial ice volume (m)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
aicen , & ! ice concentration
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) :: &
aice , & ! concentration of ice
aice0 , & ! concentration of open water
fpond ! fresh water flux to ponds (kg/m^2/s)
! local variables
integer (kind=int_kind) :: &
n, nd , & ! category indices
k ! ice layer index
real (kind=dbl_kind) :: &
slope , & ! rate of change of dhice with hice
dh0 , & ! change in ice thickness at h = 0
da0 , & ! area melting from category 1
damax , & ! max allowed reduction in category 1 area
etamin, etamax,& ! left and right limits of integration
x1 , & ! etamax - etamin
x2 , & ! (etamax^2 - etamin^2) / 2
x3 , & ! (etamax^3 - etamin^3) / 3
wk1, wk2 ! temporary variables
real (kind=dbl_kind), dimension(0:ncat) :: &
hbnew ! new boundary locations
real (kind=dbl_kind), dimension(ncat) :: &
g0 , & ! constant coefficient in g(h)
g1 , & ! linear coefficient in g(h)
hL , & ! left end of range over which g(h) > 0
hR ! right end of range over which g(h) > 0
real (kind=dbl_kind), dimension(ncat) :: &
hicen , & ! ice thickness for each cat (m)
hicen_init , & ! initial ice thickness for each cat (pre-thermo)
dhicen , & ! thickness change for remapping (m)
daice , & ! ice area transferred across boundary
dvice ! ice volume transferred across boundary
real (kind=dbl_kind), dimension (ncat) :: &
eicen, & ! energy of melting for each ice layer (J/m^2)
esnon, & ! energy of melting for each snow layer (J/m^2)
vbrin, & ! ice volume defined by brine height (m)
sicen ! Bulk salt in h ice (ppt*m)
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 categories
esno_init, esno_final, & ! snow energy summed over categories
sice_init, sice_final, & ! ice bulk salinity summed over categories
vbri_init, vbri_final ! briny ice volume summed over categories
! NOTE: Third index of donor, daice, dvice should be ncat-1,
! except that compilers would have trouble when ncat = 1
integer (kind=int_kind), dimension(ncat) :: &
donor ! donor category index
logical (kind=log_kind) :: &
remap_flag ! remap ITD if remap_flag is true
character (len=char_len) :: &
fieldid ! field identifier
logical (kind=log_kind), parameter :: &
print_diags = .false. ! if true, prints when remap_flag=F
character(len=*),parameter :: subname='(linear_itd)'
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
do n = 1, ncat
donor(n) = 0
daice(n) = c0
dvice(n) = c0
enddo
!-----------------------------------------------------------------
! Compute volume and energy sums that linear remapping should
! conserve.
!-----------------------------------------------------------------
if (conserv_check) then
do n = 1, ncat
eicen(n) = c0
esnon(n) = c0
vbrin(n) = c0
sicen(n) = c0
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+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
if (tr_brine) then
vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
endif
do k = 1, nilyr
sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
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 ! conserv_check
!-----------------------------------------------------------------
! Initialize remapping flag.
! Remapping is done wherever remap_flag = .true.
! In rare cases the category boundaries may shift too far for the
! remapping algorithm to work, and remap_flag is set to .false.
! In these cases the simpler 'rebin' subroutine will shift ice
! between categories if needed.
!-----------------------------------------------------------------
remap_flag = .true.
!-----------------------------------------------------------------
! Compute thickness change in each category.
!-----------------------------------------------------------------
do n = 1, ncat
if (aicen_init(n) > puny) then
hicen_init(n) = vicen_init(n) / aicen_init(n)
else
hicen_init(n) = c0
endif ! aicen_init > puny
if (aicen (n) > puny) then
hicen (n) = vicen(n) / aicen(n)
dhicen(n) = hicen(n) - hicen_init(n)
else
hicen (n) = c0
dhicen(n) = c0
endif ! aicen > puny
enddo ! n
!-----------------------------------------------------------------
! Compute new category boundaries, hbnew, based on changes in
! ice thickness from vertical thermodynamics.
!-----------------------------------------------------------------
hbnew(0) = hin_max(0)
do n = 1, ncat-1
if (hicen_init(n) > puny .and. &
hicen_init(n+1) > puny) then
! interpolate between adjacent category growth rates
slope = (dhicen(n+1) - dhicen(n)) / &
(hicen_init(n+1) - hicen_init(n))
hbnew(n) = hin_max(n) + dhicen(n) &
+ slope * (hin_max(n) - hicen_init(n))
elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0
hbnew(n) = hin_max(n) + dhicen(n)
elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0
hbnew(n) = hin_max(n) + dhicen(n+1)
else
hbnew(n) = hin_max(n)
endif
!-----------------------------------------------------------------
! Check that each boundary lies between adjacent values of hicen.
! If not, set remap_flag = .false.
! Write diagnosis outputs if remap_flag was changed to false
!-----------------------------------------------------------------
if (aicen(n) > puny .and. hicen(n) >= hbnew(n)) then
remap_flag = .false.
if (print_diags) then
write(warnstr,*) subname, 'ITD: hicen(n) > hbnew(n)'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'cat ',n
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hicen(n) =', hicen(n)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
call icepack_warnings_add(warnstr)
endif
elseif (aicen(n+1) > puny .and. hicen(n+1) <= hbnew(n)) then
remap_flag = .false.
if (print_diags) then
write(warnstr,*) subname, 'ITD: hicen(n+1) < hbnew(n)'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'cat ',n
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hicen(n+1) =', hicen(n+1)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
call icepack_warnings_add(warnstr)
endif
endif
!-----------------------------------------------------------------
! Check that hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
! If not, set remap_flag = .false.
! (In principle we could allow this, but it would make the code
! more complicated.)
! Write diagnosis outputs if remap_flag was changed to false
!-----------------------------------------------------------------
if (hbnew(n) > hin_max(n+1)) then
remap_flag = .false.
if (print_diags) then
write(warnstr,*) subname, 'ITD hbnew(n) > hin_max(n+1)'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'cat ',n
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hin_max(n+1) =', hin_max(n+1)
call icepack_warnings_add(warnstr)
endif
endif
if (hbnew(n) < hin_max(n-1)) then
remap_flag = .false.
if (print_diags) then
write(warnstr,*) subname, 'ITD: hbnew(n) < hin_max(n-1)'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'cat ',n
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'hin_max(n-1) =', hin_max(n-1)
call icepack_warnings_add(warnstr)
endif
endif
enddo ! boundaries, 1 to ncat-1
!-----------------------------------------------------------------
! Compute hbnew(ncat)
!-----------------------------------------------------------------
if (aicen(ncat) > puny) then
hbnew(ncat) = c3*hicen(ncat) - c2*hbnew(ncat-1)
else
hbnew(ncat) = hin_max(ncat)
endif
hbnew(ncat) = max(hbnew(ncat),hin_max(ncat-1))
!-----------------------------------------------------------------
! Identify cells where the ITD is to be remapped
!-----------------------------------------------------------------
if (remap_flag) then
!-----------------------------------------------------------------
! Compute g(h) for category 1 at start of time step
! (hicen = hicen_init)
!-----------------------------------------------------------------
call fit_line(aicen(1), hicen_init(1), &
hbnew(0), hin_max (1), &
g0 (1), g1 (1), &
hL (1), hR (1))
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Find area lost due to melting of thin (category 1) ice
!-----------------------------------------------------------------
if (aicen(1) > puny) then
dh0 = dhicen(1)
if (dh0 < c0) then ! remove area from category 1
dh0 = min(-dh0,hin_max(1)) ! dh0 --> |dh0|
!-----------------------------------------------------------------
! Integrate g(1) from 0 to dh0 to estimate area melted
!-----------------------------------------------------------------
! right integration limit (left limit = 0)
etamax = min(dh0,hR(1)) - hL(1)
if (etamax > c0) then
x1 = etamax
x2 = p5 * etamax*etamax
da0 = g1(1)*x2 + g0(1)*x1 ! ice area removed
! constrain new thickness <= hicen_init
damax = aicen(1) * (c1-hicen(1)/hicen_init(1)) ! damax > 0
da0 = min (da0, damax)
! remove area, conserving volume
hicen(1) = hicen(1) * aicen(1) / (aicen(1)-da0)
aicen(1) = aicen(1) - da0
if (tr_pond_topo) &
fpond = fpond - (da0 * trcrn(nt_apnd,1) &
* trcrn(nt_hpnd,1))
endif ! etamax > 0
else ! dh0 >= 0
hbnew(0) = min(dh0,hin_max(1)) ! shift hbnew(0) to right
endif
endif ! aicen(1) > puny
!-----------------------------------------------------------------
! Compute g(h) for each ice thickness category.
!-----------------------------------------------------------------
do n = 1, ncat
call fit_line(aicen(n), hicen(n), &
hbnew(n-1), hbnew(n), &
g0 (n), g1 (n), &
hL (n), hR (n))
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Compute area and volume to be shifted across each boundary.
!-----------------------------------------------------------------
donor(n) = 0
daice(n) = c0
dvice(n) = c0
enddo
do n = 1, ncat-1
if (hbnew(n) > hin_max(n)) then ! transfer from n to n+1
! left and right integration limits in eta space
etamin = max(hin_max(n), hL(n)) - hL(n)
etamax = min(hbnew(n), hR(n)) - hL(n)
donor(n) = n
else ! hbnew(n) <= hin_max(n); transfer from n+1 to n
! left and right integration limits in eta space
etamin = c0
etamax = min(hin_max(n), hR(n+1)) - hL(n+1)
donor(n) = n+1
endif ! hbnew(n) > hin_max(n)
if (etamax > etamin) then
x1 = etamax - etamin
wk1 = etamin*etamin
wk2 = etamax*etamax
x2 = p5 * (wk2 - wk1)
wk1 = wk1*etamin
wk2 = wk2*etamax
x3 = p333 * (wk2 - wk1)
nd = donor(n)
daice(n) = g1(nd)*x2 + g0(nd)*x1
dvice(n) = g1(nd)*x3 + g0(nd)*x2 + daice(n)*hL(nd)
endif ! etamax > etamin
! If daice or dvice is very small, shift no ice.
nd = donor(n)
if (daice(n) < aicen(nd)*puny) then
daice(n) = c0
dvice(n) = c0
donor(n) = 0
endif
if (dvice(n) < vicen(nd)*puny) then
daice(n) = c0
dvice(n) = c0
donor(n) = 0
endif
! If daice is close to aicen or dvice is close to vicen,
! shift entire category
if (daice(n) > aicen(nd)*(c1-puny)) then
daice(n) = aicen(nd)
dvice(n) = vicen(nd)
endif
if (dvice(n) > vicen(nd)*(c1-puny)) then
daice(n) = aicen(nd)
dvice(n) = vicen(nd)
endif
enddo ! boundaries, 1 to ncat-1
!-----------------------------------------------------------------
! Shift ice between categories as necessary
!-----------------------------------------------------------------
! maintain qsno negative definiteness
do n = 1, ncat
do k = nt_qsno, nt_qsno+nslyr-1
trcrn(k,n) = trcrn(k,n) + rhos*Lfresh
enddo
enddo
call shift_ice (ntrcr, ncat, &
trcr_depend, &
trcr_base, &
n_trcr_strata, &
nt_strata, &
aicen, trcrn, &
vicen, vsnon, &
hicen, donor, &
daice, dvice )
if (icepack_warnings_aborted(subname)) return
! maintain qsno negative definiteness
do n = 1, ncat
do k = nt_qsno, nt_qsno+nslyr-1
trcrn(k,n) = trcrn(k,n) - rhos*Lfresh
enddo
enddo
!-----------------------------------------------------------------
! Make sure hice(1) >= minimum ice thickness hi_min.
!-----------------------------------------------------------------
if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then
da0 = aicen(1) * (c1 - hicen(1)/hi_min)
aicen(1) = aicen(1) - da0
hicen(1) = hi_min
if (tr_pond_topo) &
fpond = fpond - (da0 * trcrn(nt_apnd,1) &
* trcrn(nt_hpnd,1))
endif
endif ! remap_flag
!-----------------------------------------------------------------
! Update fractional ice area in each grid cell.
!-----------------------------------------------------------------
call aggregate_area (ncat, aicen, aice, aice0)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Check volume and energy conservation.
!-----------------------------------------------------------------
if (conserv_check) then
do n = 1, ncat
eicen(n) = c0
esnon(n) = c0
vbrin(n) = c0
sicen(n) = c0
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+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
if (tr_brine) then
vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
endif
do k = 1, nilyr
sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
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
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//':sicen'
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 ! conservation check
end subroutine linear_itd
!=======================================================================
!
! Fit g(h) with a line, satisfying area and volume constraints.
! To reduce roundoff errors caused by large values of g0 and g1,
! we actually compute g(eta), where eta = h - hL, and hL is the
! left boundary.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
subroutine fit_line (aicen, hice, &
hbL, hbR, &
g0, g1, &
hL, hR)
real (kind=dbl_kind), intent(in) :: &
aicen ! concentration of ice
real (kind=dbl_kind), intent(in) :: &
hbL, hbR , & ! left and right category boundaries
hice ! ice thickness
real (kind=dbl_kind), intent(out):: &
g0, g1 , & ! coefficients in linear equation for g(eta)
hL , & ! min value of range over which g(h) > 0
hR ! max value of range over which g(h) > 0
! local variables
real (kind=dbl_kind) :: &
h13 , & ! hbL + 1/3 * (hbR - hbL)
h23 , & ! hbL + 2/3 * (hbR - hbL)
dhr , & ! 1 / (hR - hL)
wk1, wk2 ! temporary variables
character(len=*),parameter :: subname='(fit_line)'
!-----------------------------------------------------------------
! Compute g0, g1, hL, and hR for each category to be remapped.
!-----------------------------------------------------------------
if (aicen > puny .and. hbR - hbL > puny) then
! Initialize hL and hR
hL = hbL
hR = hbR
! Change hL or hR if hicen(n) falls outside central third of range
h13 = p333 * (c2*hL + hR)
h23 = p333 * (hL + c2*hR)
if (hice < h13) then
hR = c3*hice - c2*hL
elseif (hice > h23) then
hL = c3*hice - c2*hR
endif
! Compute coefficients of g(eta) = g0 + g1*eta
dhr = c1 / (hR - hL)
wk1 = c6 * aicen * dhr
wk2 = (hice - hL) * dhr
g0 = wk1 * (p666 - wk2)
g1 = c2*dhr * wk1 * (wk2 - p5)
else
g0 = c0
g1 = c0
hL = c0
hR = c0
endif ! aicen > puny
end subroutine fit_line
!=======================================================================
!
! Given some added new ice to the base of the existing ice, recalculate
! vertical tracer so that new grid cells are all the same size.
!
! author: A. K. Turner, LANL
!
subroutine update_vertical_tracers(nilyr, trc, h1, h2, trc0)
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real (kind=dbl_kind), dimension(:), intent(inout) :: &
trc ! vertical tracer
real (kind=dbl_kind), intent(in) :: &
h1, & ! old thickness
h2, & ! new thickness
trc0 ! tracer value of added ice on ice bottom
! local variables
real(kind=dbl_kind), dimension(nilyr) :: trc2 ! updated tracer temporary
! vertical indices for old and new grid
integer :: k1, k2
real (kind=dbl_kind) :: &
z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom
z2a, z2b, & ! upper, lower boundary of new cell
overlap , & ! overlap between old and new cell
rnilyr
character(len=*),parameter :: subname='(update_vertical_tracers)'
rnilyr = real(nilyr,dbl_kind)
! loop over new grid cells
do k2 = 1, nilyr
! initialize new tracer
trc2(k2) = c0
! calculate upper and lower boundary of new cell
z2a = ((k2 - 1) * h2) / rnilyr
z2b = (k2 * h2) / rnilyr
! loop over old grid cells
do k1 = 1, nilyr
! calculate upper and lower boundary of old cell
z1a = ((k1 - 1) * h1) / rnilyr
z1b = (k1 * h1) / rnilyr
! calculate overlap between old and new cell
overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
! aggregate old grid cell contribution to new cell
trc2(k2) = trc2(k2) + overlap * trc(k1)
enddo ! k1
! calculate upper and lower boundary of added new ice at bottom
z1a = h1
z1b = h2
! calculate overlap between added ice and new cell
overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
! aggregate added ice contribution to new cell
trc2(k2) = trc2(k2) + overlap * trc0
! renormalize new grid cell
trc2(k2) = (rnilyr * trc2(k2)) / h2
enddo ! k2
! update vertical tracer array with the adjusted tracer
trc = trc2
end subroutine update_vertical_tracers
!=======================================================================
!
! Given the fraction of ice melting laterally in each grid cell
! (computed in subroutine frzmlt_bottom_lateral), melt ice.
!
! author: C. M. Bitz, UW
! 2003: Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
! 2016 Lettie Roach, NIWA/VUW, added floe size dependence
!
subroutine lateral_melt (dt, ncat, &
nilyr, nslyr, &
n_aero, &
fpond, &
fresh, fsalt, &
fhocn, faero_ocn, &
fiso_ocn, &
rside, meltl, &
fside, sss, &
aicen, vicen, &
vsnon, trcrn, &
fzsal, flux_bio, &
nbtrcr, nblyr, &
nfsd, d_afsd_latm,&
floe_rad_c, floe_binwidth)
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd , & ! number of floe size categories
nilyr , & ! number of ice layers
nblyr , & ! number of bio layers
nslyr , & ! number of snow layers
n_aero , & ! number of aerosol tracers
nbtrcr ! number of bio tracers
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 ! tracer array
real (kind=dbl_kind), intent(in) :: &
rside , & ! fraction of ice that melts laterally
fside ! lateral heat flux (W/m^2)
real (kind=dbl_kind), intent(inout) :: &
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)
meltl , & ! lateral ice melt (m/step-->cm/day)
fzsal ! salt flux from zsalinity (kg/m2/s)
real (kind=dbl_kind), dimension (:), intent(in) :: &
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(out) :: &
d_afsd_latm ! change in fsd due to lateral melt (m)
real (kind=dbl_kind), dimension(nbtrcr), &
intent(inout) :: &
flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
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
integer (kind=int_kind) :: &
n , & ! thickness category index
k , & ! layer index
nsubt ! sub timesteps for FSD tendency
real (kind=dbl_kind) :: &
dfhocn , & ! change in fhocn
dfpond , & ! change in fpond
dfresh , & ! change in fresh
dfsalt , & ! change in fsalt
dvssl , & ! snow surface layer volume
dvint , & ! snow interior layer
cat1_arealoss, tmp !
logical (kind=log_kind) :: &
flag ! .true. if there could be lateral melting
real (kind=dbl_kind), dimension (ncat) :: &
aicen_init, & ! initial area fraction
vicen_init, & ! volume per unit area of ice (m)
G_radialn , & ! rate of lateral melt (m/s)
delta_an , & ! change in the ITD
qin , & ! enthalpy
rsiden ! delta_an/aicen
real (kind=dbl_kind), dimension (nfsd,ncat) :: &
afsdn , & ! floe size distribution tracer
afsdn_init ! initial value
real (kind=dbl_kind), dimension (nfsd) :: &
df_flx, & ! finite difference for FSD
afsd_tmp, d_afsd_tmp
real (kind=dbl_kind), dimension(nfsd+1) :: &
f_flx !
!echmod - for average qin
real (kind=dbl_kind), intent(in) :: &
sss
real (kind=dbl_kind) :: &
Ti, Si0, qi0, &
elapsed_t, & ! FSD subcycling
subdt ! FSD timestep (s)
character(len=*), parameter :: subname='(lateral_melt)'
flag = .false.
dfhocn = c0
dfpond = c0
dfresh = c0
dfsalt = c0
dvssl = c0
dvint = c0
cat1_arealoss = c0
tmp = c0
vicen_init = c0
G_radialn = c0
delta_an = c0
qin = c0
rsiden = c0
afsdn = c1
afsdn_init = c0
df_flx = c0
f_flx = c0
if (tr_fsd) then
call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
if (icepack_warnings_aborted(subname)) return
afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
aicen_init = aicen
afsdn_init = afsdn ! for diagnostics
d_afsd_latm(:) = c0
end if
if (tr_fsd .and. fside < c0) then
flag = .true.
! echmod - using category values would be preferable to the average value
! Compute average enthalpy of ice (taken from add_new_ice)
if (sss > c2 * dSin0_frazil) then
Si0 = sss - dSin0_frazil
else
Si0 = sss**2 / (c4*dSin0_frazil)
endif
Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
qi0 = enthalpy_mush(Ti, Si0)
do n = 1, ncat
if (ktherm == 2) then ! mushy
do k = 1, nilyr
!qin(n) = qin(n) &
! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind)
qin(n) = qi0
enddo
else
qin(n) = -rhoi*Lfresh
endif
if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative
if (G_radialn(n) < -puny) then
if (any(afsdn(:,n) < c0)) print*,&
'lateral_melt B afsd < 0',n
cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
* G_radialn(n) / floe_binwidth(1)
delta_an(n) = c0
do k = 1, nfsd
delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
* trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
end do
! add negative area loss from fsd
delta_an(n) = delta_an(n) - cat1_arealoss
if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
! following original code, not necessary for fsd
if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)
end if ! G_radialn
enddo ! ncat
else if (rside > c0) then ! original, non-fsd implementation
flag = .true.
rsiden(:) = rside ! initialize
endif
if (flag) then ! grid cells with lateral melting.
do n = 1, ncat
!-----------------------------------------------------------------
! Melt the ice and increment fluxes.
!-----------------------------------------------------------------
! fluxes to coupler
! dfresh > 0, dfsalt > 0, dfpond > 0
dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt
dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
fresh = fresh + dfresh
fsalt = fsalt + dfsalt
if (tr_pond_topo) then
dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
fpond = fpond - dfpond
endif
! history diagnostics
meltl = meltl + vicen(n)*rsiden(n)
! state variables
vicen_init(n) = vicen(n)
aicen(n) = aicen(n) * (c1 - rsiden(n))
vicen(n) = vicen(n) * (c1 - rsiden(n))
vsnon(n) = vsnon(n) * (c1 - rsiden(n))
! floe size distribution
if (tr_fsd) then
if (rsiden(n) > puny) then
if (aicen(n) > puny) then
! adaptive subtimestep
elapsed_t = c0
afsd_tmp(:) = afsdn_init(:,n)
d_afsd_tmp(:) = c0
nsubt = 0
DO WHILE (elapsed_t.lt.dt)
nsubt = nsubt + 1
if (nsubt.gt.100) &
print *, 'latm not converging'
! finite differences
df_flx(:) = c0
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
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*,'sum(df_flx)/=0'
! this term ensures area conservation
tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
! fsd tendency
do k = 1, nfsd
d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
* (c1/floe_rad_c(k) - tmp)
end do
! timestep required for this
subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
subdt = MIN(subdt, dt)
! update fsd and elapsed time
afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
elapsed_t = elapsed_t + subdt
END DO
afsdn(:,n) = afsd_tmp(:)
end if ! aicen
end if ! rside > 0, otherwise do nothing
end if ! tr_fsd
! fluxes
do k = 1, nilyr
! enthalpy tracers do not change (e/v constant)
! heat flux to coupler for ice melt (dfhocn < 0)
dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
* vicen(n)/real(nilyr,kind=dbl_kind)
fhocn = fhocn + dfhocn
enddo ! nilyr
do k = 1, nslyr
! heat flux to coupler for snow melt (dfhocn < 0)
dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
* vsnon(n)/real(nslyr,kind=dbl_kind)
fhocn = fhocn + dfhocn
enddo ! nslyr
if (tr_aero) then
do k = 1, n_aero
faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
*(trcrn(nt_aero +4*(k-1),n) &
+ trcrn(nt_aero+1+4*(k-1),n)) &
+ vicen(n) &
*(trcrn(nt_aero+2+4*(k-1),n) &
+ trcrn(nt_aero+3+4*(k-1),n))) &
* rsiden(n) / dt
enddo ! k
endif ! tr_aero
if (tr_iso) then
do k = 1, n_iso
fiso_ocn(k) = fiso_ocn(k) &
+ (vsnon(n)*trcrn(nt_isosno+k-1,n) &
+ vicen(n)*trcrn(nt_isoice+k-1,n)) &
* rside / dt
enddo ! k
endif ! tr_iso
!-----------------------------------------------------------------
! Biogeochemistry
!-----------------------------------------------------------------
if (z_tracers) then ! snow tracers
dvssl = min(p5*vsnon(n), hs_ssl*aicen(n)) !snow surface layer
dvint = vsnon(n)- dvssl !snow interior
do k = 1, nbtrcr
flux_bio(k) = flux_bio(k) &
+ (trcrn(bio_index(k)+nblyr+1,n)*dvssl &
+ trcrn(bio_index(k)+nblyr+2,n)*dvint) &
* rsiden(n) / dt
enddo
endif
enddo ! n
if (solve_zsal .or. z_tracers) &
call lateral_melt_bgc(dt, &
ncat, nblyr, &
rside, vicen_init, & !echmod: use rsiden
trcrn, fzsal, &
flux_bio, nbtrcr)
if (icepack_warnings_aborted(subname)) return
endif ! flag
if (tr_fsd) then
call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
if (icepack_warnings_aborted(subname)) return
! diagnostics
do k = 1, nfsd
d_afsd_latm(k) = c0
do n = 1, ncat
d_afsd_latm(k) = d_afsd_latm(k) &
+ afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
end do
end do
end if
end subroutine lateral_melt
!=======================================================================
!
! Given the volume of new ice grown in open water, compute its area
! and thickness and add it to the appropriate category or categories.
!
! NOTE: Usually all the new ice is added to category 1. An exception is
! made if there is no open water or if the new ice is too thick
! for category 1, in which case ice is distributed evenly over the
! entire cell. Subroutine rebin should be called in case the ice
! thickness lies outside category bounds after new ice formation.
!
! When ice must be added to categories above category 1, the mushy
! formulation (ktherm=2) adds it only to the bottom of the ice. When
! added to only category 1, all formulations combine the new ice and
! existing ice tracers as bulk quantities.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
! Adrian Turner, LANL
! Lettie Roach, NIWA/VUW
!
subroutine add_new_ice (ncat, nilyr, &
nfsd, nblyr, &
n_aero, dt, &
ntrcr, nltrcr, &
hin_max, ktherm, &
aicen, trcrn, &
vicen, vsnon1, &
aice0, aice, &
frzmlt, frazil, &
frz_onset, yday, &
update_ocn_f, &
fresh, fsalt, &
Tf, sss, &
salinz, phi_init, &
dSin0_frazil, &
bgrid, cgrid, igrid, &
nbtrcr, flux_bio, &
ocean_bio, fzsal, &
frazil_diag, &
fiso_ocn, &
HDO_ocn, H2_16O_ocn, &
H2_18O_ocn, &
wave_sig_ht, &
wave_spectrum, &
wavefreq, &
dwavefreq, &
d_afsd_latg, &
d_afsd_newi, &
floe_rad_c, floe_binwidth)
use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd , & ! number of floe size categories
nilyr , & ! number of ice layers
nblyr , & ! number of bio layers
ntrcr , & ! number of tracers
nltrcr, & ! number of zbgc tracers
n_aero, & ! number of aerosol tracers
ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
hin_max ! category boundaries (m)
real (kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
aice , & ! total concentration of ice
frzmlt, & ! freezing/melting potential (W/m^2)
Tf , & ! freezing temperature (C)
sss , & ! sea surface salinity (ppt)
vsnon1 ! category 1 snow volume per ice area (m)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn ! ice tracers
! 1: surface temperature
real (kind=dbl_kind), intent(inout) :: &
aice0 , & ! concentration of open water
frazil , & ! frazil ice growth (m/step-->cm/day)
frazil_diag,& ! frazil ice growth diagnostic (m/step-->cm/day)
fresh , & ! fresh water flux to ocean (kg/m^2/s)
fsalt ! salt flux to ocean (kg/m^2/s)
real (kind=dbl_kind), intent(inout), optional :: &
frz_onset ! day of year that freezing begins (congel or frazil)
real (kind=dbl_kind), intent(in), optional :: &
yday ! day of year
real (kind=dbl_kind), dimension(:), intent(in) :: &
salinz ! initial salinity profile
real (kind=dbl_kind), intent(in) :: &
phi_init , & ! initial frazil liquid fraction
dSin0_frazil ! initial frazil bulk salinity reduction from sss
logical (kind=log_kind), intent(in) :: &
update_ocn_f ! if true, update fresh water and salt fluxes
! BGC
real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
bgrid ! biology nondimensional vertical grid points
real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
igrid ! biology vertical interface points
real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
cgrid ! CICE vertical coordinate
integer (kind=int_kind), intent(in) :: &
nbtrcr ! number of biology tracers
real (kind=dbl_kind), dimension (:), intent(inout) :: &
flux_bio ! tracer flux to ocean from biology (mmol/m^2/s)
real (kind=dbl_kind), dimension (:), intent(in) :: &
ocean_bio ! ocean concentration of biological tracer
! zsalinity
real (kind=dbl_kind), intent(inout) :: &
fzsal ! salt flux to ocean from zsalinity (kg/m^2/s)
! water isotopes
real (kind=dbl_kind), dimension(:), intent(inout) :: &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)
real (kind=dbl_kind), intent(in) :: &
HDO_ocn , & ! ocean concentration of HDO (kg/kg)
H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
! floe size distribution
real (kind=dbl_kind), intent(in) :: &
wave_sig_ht ! significant height of waves globally (m)
real (kind=dbl_kind), dimension(:), intent(in) :: &
wave_spectrum ! ocean surface wave spectrum, E(f) (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) :: &
floe_rad_c , & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
real (kind=dbl_kind), dimension(ncat) :: & ! for now
! change in thickness distribution (area)
d_an_latg , & ! due to fsd lateral growth
d_an_newi ! new ice formation
real (kind=dbl_kind), dimension(:), intent(out) :: &
! change in thickness distribution (area)
d_afsd_latg , & ! due to fsd lateral growth
d_afsd_newi ! new ice formation
! local variables
integer (kind=int_kind) :: &
ncats , & ! max categories to which new ice is added, initially
n , & ! ice category index
k , & ! ice layer index
it ! aerosol tracer index
real (kind=dbl_kind) :: &
ai0new , & ! area of new ice added to cat 1
vi0new , & ! volume of new ice added to cat 1
vi0new_lat , & ! volume of new ice added laterally to fsd
hsurp , & ! thickness of new ice added to each cat
fnew , & ! heat flx to open water for new ice (W/m^2)
hi0new , & ! thickness of new ice
hi0max , & ! max allowed thickness of new ice
vsurp , & ! volume of new ice added to each cat
vtmp , & ! total volume of new and old ice
area1 , & ! starting fractional area of existing ice
alvl , & ! starting level ice area
dfresh , & ! change in fresh
dfsalt , & ! change in fsalt
vi0tmp , & ! frzmlt part of frazil
Ti , & ! frazil temperature
qi0new , & ! frazil ice enthalpy
Si0new , & ! frazil ice bulk salinity
vi0_init , & ! volume of new ice
vice1 , & ! starting volume of existing ice
vice_init, vice_final, & ! ice volume summed over categories
eice_init, eice_final ! ice energy summed over categories
real (kind=dbl_kind) :: frazil_conc
real (kind=dbl_kind), dimension (nilyr) :: &
Sprofile ! salinity profile used for new ice additions
character (len=char_len) :: &
fieldid ! field identifier
real (kind=dbl_kind), dimension (ncat) :: &
eicen, & ! energy of melting for each ice layer (J/m^2)
aicen_init, & ! fractional area of ice
vicen_init ! volume per unit area of ice (m)
! floe size distribution
real (kind=dbl_kind), dimension (nfsd,ncat) :: &
afsdn ! floe size distribution tracer (originally areal_mfstd_init)
! real (kind=dbl_kind), dimension (nfsd) :: &
! afsd , & ! fsd tracer for each thickness category
real (kind=dbl_kind), dimension (ncat) :: &
d_an_tot, & ! change in the ITD due to lateral growth and new ice
area2 ! area after lateral growth and before new ice formation
real (kind=dbl_kind), dimension (ncat) :: &
vin0new ! volume of new ice added to any thickness cat
real (kind=dbl_kind), dimension (nfsd) :: &
afsd_ni ! areal mFSTD after new ice added
real (kind=dbl_kind) :: &
tmp, &
latsurf_area, & ! fractional area of ice on sides of floes
lead_area , & ! fractional area of ice in lead region
G_radial , & ! lateral melt rate (m/s)
tot_latg , & ! total fsd lateral growth in open water
ai0mod ! ai0new - tot_latg
character(len=*),parameter :: subname='(add_new_ice)'
!-----------------------------------------------------------------
! initialize
!-----------------------------------------------------------------
hsurp = c0
hi0new = c0
ai0new = c0
afsdn(:,:) = c0
d_an_latg(:) = c0
d_an_tot(:) = c0
d_an_newi(:) = c0
vin0new(:) = c0
if (tr_fsd) then
d_afsd_latg(:) = c0 ! diagnostics
d_afsd_newi(:) = c0
end if
area2(:) = aicen(:)
lead_area = c0
latsurf_area = c0
G_radial = c0
tot_latg = c0
if (ncat > 1) then
hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary
else
hi0max = bignum ! big number
endif
if (tr_fsd) then
call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
if (icepack_warnings_aborted(subname)) return
endif
do n = 1, ncat
aicen_init(n) = aicen(n)
vicen_init(n) = vicen(n)
eicen(n) = c0
if (tr_fsd) then
do k = 1, nfsd
afsdn(k,n) = trcrn(nt_fsd+k-1,n)
enddo
endif
enddo
if (conserv_check) then
do n = 1, ncat
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
enddo
call column_sum (ncat, vicen, vice_init)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, eicen, eice_init)
if (icepack_warnings_aborted(subname)) return
endif ! conserv_check
!-----------------------------------------------------------------
! Compute average enthalpy of new ice.
! Sprofile is the salinity profile used when adding new ice to
! all categories, for ktherm/=2, and to category 1 for all ktherm.
!
! NOTE: POP assumes new ice is fresh!
!-----------------------------------------------------------------
if (ktherm == 2) then ! mushy
if (sss > c2 * dSin0_frazil) then
Si0new = sss - dSin0_frazil
else
Si0new = sss**2 / (c4*dSin0_frazil)
endif
do k = 1, nilyr
Sprofile(k) = Si0new
enddo
Ti = min(liquidus_temperature_mush(Si0new/phi_init), -p1)
qi0new = enthalpy_mush(Ti, Si0new)
else
do k = 1, nilyr
Sprofile(k) = salinz(k)
enddo
qi0new = -rhoi*Lfresh
endif ! ktherm
!-----------------------------------------------------------------
! Compute the volume, area, and thickness of new ice.
!-----------------------------------------------------------------
fnew = max (frzmlt, c0) ! fnew > 0 iff frzmlt > 0
vi0new = -fnew*dt / qi0new ! note sign convention, qi < 0
vi0_init = vi0new ! for bgc
! increment ice volume and energy
if (conserv_check) then
vice_init = vice_init + vi0new
eice_init = eice_init + vi0new*qi0new
endif
! history diagnostics
frazil = vi0new
if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss
if (present(frz_onset) .and. present(yday)) then
if (frazil > puny .and. frz_onset < puny) frz_onset = yday
endif
!-----------------------------------------------------------------
! Update fresh water and salt fluxes.
!
! NOTE: POP assumes fresh water and salt flux due to frzmlt > 0
! is NOT included in fluxes fresh and fsalt.
!-----------------------------------------------------------------
if (update_ocn_f) then
if (ktherm <= 1) then
dfresh = -rhoi*vi0new/dt
dfsalt = ice_ref_salinity*p001*dfresh
fresh = fresh + dfresh
fsalt = fsalt + dfsalt
! elseif (ktherm == 2) the fluxes are added elsewhere
endif
else ! update_ocn_f = false
if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
vi0tmp = fnew*dt / (rhoi*Lfresh)
dfresh = -rhoi*(vi0new - vi0tmp)/dt
dfsalt = ice_ref_salinity*p001*dfresh
fresh = fresh + dfresh
fsalt = fsalt + dfsalt
frazil_diag = frazil - vi0tmp
! elseif ktherm==1 do nothing
endif
endif
!-----------------------------------------------------------------
! Decide how to distribute the new ice.
!-----------------------------------------------------------------
if (vi0new > c0) then
if (tr_fsd) & ! lateral growth of existing ice
! calculate change in conc due to lateral growth
! update vi0new, without change to afsdn or aicen
call 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)
if (icepack_warnings_aborted(subname)) return
ai0mod = aice0
! separate frazil ice growth from lateral ice growth
if (tr_fsd) ai0mod = aice0-tot_latg
! new ice area and thickness
! hin_max(0) < new ice thickness < hin_max(1)
if (ai0mod > puny) then
hi0new = max(vi0new/ai0mod, hfrazilmin)
if (hi0new > hi0max .and. ai0mod+puny < c1) then
! distribute excess volume over all categories (below)
hi0new = hi0max
ai0new = ai0mod
vsurp = vi0new - ai0new*hi0new
hsurp = vsurp / aice
vi0new = ai0new*hi0new
else
! put ice in a single category, with hsurp = 0
ai0new = vi0new/hi0new
endif
else ! aice0 < puny
hsurp = vi0new/aice ! new thickness in each cat
vi0new = c0
endif ! aice0 > puny
! volume added to each category from lateral growth
do n = 1, ncat
if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
end do
! combine areal change from new ice growth and lateral growth
d_an_newi(1) = ai0new
d_an_tot(2:ncat) = d_an_latg(2:ncat)
d_an_tot(1) = d_an_latg(1) + d_an_newi(1)
if (tr_fsd) then
vin0new(1) = vin0new(1) + ai0new*hi0new ! not BFB
else
vin0new(1) = vi0new
endif
endif ! vi0new > puny
!-----------------------------------------------------------------
! Distribute excess ice volume among ice categories by increasing
! ice thickness, leaving ice area unchanged.
!
! NOTE: If new ice contains globally conserved tracers
! (e.g., isotopes from seawater), code must be added here.
!
! The mushy formulation (ktherm=2) puts the new ice only at the
! bottom of existing ice and adjusts the layers accordingly.
! The other formulations distribute the new ice throughout the
! existing ice column.
!-----------------------------------------------------------------
if (hsurp > c0) then ! add ice to all categories
do n = 1, ncat
vsurp = hsurp * aicen(n)
! update ice age due to freezing (new ice age = dt)
vtmp = vicen(n) + vsurp
if (tr_iage .and. vtmp > puny) &
trcrn(nt_iage,n) = &
(trcrn(nt_iage,n)*vicen(n) + dt*vsurp) / vtmp
if (tr_lvl .and. vicen(n) > puny) then
trcrn(nt_vlvl,n) = &
(trcrn(nt_vlvl,n)*vicen(n) + &
trcrn(nt_alvl,n)*vsurp) / vtmp
endif
if (tr_aero .and. vtmp > puny) then
do it = 1, n_aero
trcrn(nt_aero+2+4*(it-1),n) = &
trcrn(nt_aero+2+4*(it-1),n)*vicen(n) / vtmp
trcrn(nt_aero+3+4*(it-1),n) = &
trcrn(nt_aero+3+4*(it-1),n)*vicen(n) / vtmp
enddo
endif
frazil_conc = c0
if (tr_iso .and. vtmp > puny) then
do it=1,n_iso
if (it==1) &
frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
if (it==2) &
frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
if (it==3) &
frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
! dilution and uptake in the ice
trcrn(nt_isoice+it-1,n) &
= (trcrn(nt_isoice+it-1,n)*vicen(n) &
+ frazil_conc*rhoi*vsurp) &
/ vtmp
fiso_ocn(it) = fiso_ocn(it) &
- frazil_conc*rhoi*vsurp/dt
enddo
endif
! update category volumes
vicen(n) = vtmp
if (ktherm == 2) then
vsurp = hsurp * aicen(n) ! note - save this above?
vtmp = vicen(n) - vsurp ! vicen is the new volume
if (vicen(n) > c0) then
call update_vertical_tracers(nilyr, &
trcrn(nt_qice:nt_qice+nilyr-1,n), &
vtmp, vicen(n), qi0new)
if (icepack_warnings_aborted(subname)) return
call update_vertical_tracers(nilyr, &
trcrn(nt_sice:nt_sice+nilyr-1,n), &
vtmp, vicen(n), Si0new)
if (icepack_warnings_aborted(subname)) return
endif
else
do k = 1, nilyr
! factor of nilyr cancels out
vsurp = hsurp * aicen(n) ! note - save this above?
vtmp = vicen(n) - vsurp ! vicen is the new volume
if (vicen(n) > c0) then
! enthalpy
trcrn(nt_qice+k-1,n) = &
(trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n)
! salinity
if (.not. solve_zsal) &
trcrn(nt_sice+k-1,n) = &
(trcrn(nt_sice+k-1,n)*vtmp + Sprofile(k)*vsurp) / vicen(n)
endif
enddo ! k
endif ! ktherm
enddo ! n
endif ! hsurp > 0
!-----------------------------------------------------------------
! Combine new ice grown in open water with ice categories.
! Using the floe size distribution, ice is added laterally to all
! categories; otherwise it is added to category 1.
! Assume that vsnon and esnon are unchanged.
! The mushy formulation assumes salt from frazil is added uniformly
! to category 1, while the others use a salinity profile.
!-----------------------------------------------------------------
ncats = 1 ! add new ice to category 1 by default
if (tr_fsd) ncats = ncat ! add new ice laterally to all categories
do n = 1, ncats
if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n
area1 = aicen(n) ! save
vice1 = vicen(n) ! save
area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi
aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth
aice0 = aice0 - d_an_tot(n)
vicen(n) = vicen(n) + vin0new(n)
trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n)
trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0)
if (tr_FY) then
trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n)
trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1)
endif
if (tr_fsd) & ! evolve the floe size distribution
! both new frazil ice and lateral growth
call 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)
if (icepack_warnings_aborted(subname)) return
if (vicen(n) > puny) then
if (tr_iage) &
trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n)
if (tr_aero) then
do it = 1, n_aero
trcrn(nt_aero+2+4*(it-1),n) = &
trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n)
trcrn(nt_aero+3+4*(it-1),n) = &
trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n)
enddo
endif
frazil_conc = c0
if (tr_iso) then
do it=1,n_iso
if (it==1) &
frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
if (it==2) &
frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
if (it==3) &
frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
trcrn(nt_isoice+it-1,1) &
= (trcrn(nt_isoice+it-1,1)*vice1) &
+ frazil_conc*rhoi*vi0new/vicen(1)
fiso_ocn(it) = fiso_ocn(it) &
- frazil_conc*rhoi*vi0new/dt
enddo
endif
if (tr_lvl) then
alvl = trcrn(nt_alvl,n)
trcrn(nt_alvl,n) = &
(trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n)
trcrn(nt_vlvl,n) = &
(trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n)
endif
if (tr_pond_cesm .or. tr_pond_topo) then
trcrn(nt_apnd,n) = &
trcrn(nt_apnd,n)*area1/aicen(n)
elseif (tr_pond_lvl) then
if (trcrn(nt_alvl,n) > puny) then
trcrn(nt_apnd,n) = &
trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n))
endif
endif
endif
do k = 1, nilyr
if (vicen(n) > c0) then
! factor of nilyr cancels out
! enthalpy
trcrn(nt_qice+k-1,n) = &
(trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n)
! salinity
if (.NOT. solve_zsal)&
trcrn(nt_sice+k-1,n) = &
(trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n)
endif
enddo
endif ! vi0new > 0
enddo ! ncats
if (conserv_check) then
do n = 1, ncat
eicen(n) = c0
do k = 1, nilyr
eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
* vicen(n)/real(nilyr,kind=dbl_kind)
enddo
enddo
call column_sum (ncat, vicen, vice_final)
if (icepack_warnings_aborted(subname)) return
call column_sum (ncat, eicen, eice_final)
if (icepack_warnings_aborted(subname)) return
fieldid = subname//':vice'
call column_conservation_check (fieldid, &
vice_init, vice_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
endif ! conserv_check
!-----------------------------------------------------------------
! Biogeochemistry
!-----------------------------------------------------------------
if (tr_brine .or. nbtrcr > 0) &
call add_new_ice_bgc(dt, nblyr, &
ncat, nilyr, nltrcr, &
bgrid, cgrid, igrid, &
aicen_init, vicen_init, vi0_init, &
aicen, vicen, vsnon1, &
vi0new, ntrcr, trcrn, &
nbtrcr, sss, ocean_bio,&
flux_bio, hsurp)
if (icepack_warnings_aborted(subname)) return
end subroutine add_new_ice
!=======================================================================
!autodocument_start icepack_step_therm2
! Driver for thermodynamic changes not needed for coupling:
! transport in thickness space, lateral growth and melting.
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
nilyr, nslyr, &
hin_max, nblyr, &
aicen, &
vicen, vsnon, &
aicen_init, vicen_init, &
trcrn, &
aice0, aice, &
trcr_depend, &
trcr_base, n_trcr_strata, &
nt_strata, &
Tf, sss, &
salinz, &
rside, meltl, &
fside, &
frzmlt, frazil, &
frain, fpond, &
fresh, fsalt, &
fhocn, update_ocn_f, &
bgrid, cgrid, &
igrid, faero_ocn, &
first_ice, fzsal, &
flux_bio, ocean_bio, &
frazil_diag, &
frz_onset, yday, &
fiso_ocn, HDO_ocn, &
H2_16O_ocn, H2_18O_ocn, &
nfsd, wave_sig_ht, &
wave_spectrum, &
wavefreq, &
dwavefreq, &
d_afsd_latg, d_afsd_newi, &
d_afsd_latm, d_afsd_weld, &
floe_rad_c, floe_binwidth)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
nfsd , & ! number of floe size categories
nltrcr , & ! number of zbgc tracers
nblyr , & ! number of bio layers
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical (kind=log_kind), intent(in) :: &
update_ocn_f ! if true, update fresh water and salt fluxes
real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max ! category boundaries (m)
real (kind=dbl_kind), intent(in) :: &
dt , & ! time step
Tf , & ! freezing temperature (C)
sss , & ! sea surface salinity (ppt)
rside , & ! fraction of ice that melts laterally
fside , & ! lateral heat flux (W/m^2)
frzmlt , & ! freezing/melting potential (W/m^2)
wave_sig_ht ! significant height of waves in ice (m)
real (kind=dbl_kind), dimension(:), intent(in) :: &
wave_spectrum ! ocean surface wave spectrum E(f) (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) :: &
floe_rad_c , & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
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 (nblyr+2), intent(in) :: &
bgrid ! biology nondimensional vertical grid points
real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
igrid ! biology vertical interface points
real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
cgrid ! CICE vertical coordinate
real (kind=dbl_kind), dimension(:), intent(in) :: &
salinz , & ! initial salinity profile
ocean_bio ! ocean concentration of biological tracer
real (kind=dbl_kind), intent(inout) :: &
aice , & ! sea ice concentration
aice0 , & ! concentration of open water
frain , & ! rainfall rate (kg/m^2 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 , & ! salt flux to ocean from zsalinity (kg/m^2/s)
meltl , & ! lateral ice melt (m/step-->cm/day)
frazil , & ! frazil ice growth (m/step-->cm/day)
frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
aicen_init,& ! initial concentration of ice
vicen_init,& ! initial volume per unit area of ice (m)
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon , & ! volume per unit area of snow (m)
faero_ocn, & ! aerosol flux to ocean (kg/m^2/s)
flux_bio ! all bio fluxes to ocean
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
trcrn ! tracers
logical (kind=log_kind), dimension(:), intent(inout) :: &
first_ice ! true until ice forms
real (kind=dbl_kind), dimension(:), intent(out) :: &
! change in floe size distribution (area)
d_afsd_latg , & ! due to fsd lateral growth
d_afsd_newi , & ! new ice formation
d_afsd_latm , & ! lateral melt
d_afsd_weld ! welding
real (kind=dbl_kind), intent(inout), optional :: &
frz_onset ! day of year that freezing begins (congel or frazil)
real (kind=dbl_kind), intent(in), optional :: &
yday ! day of year
! water isotopes
real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)
real (kind=dbl_kind), optional, intent(in) :: &
HDO_ocn , & ! ocean concentration of HDO (kg/kg)
H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
!autodocument_end
! local variables
! water isotopes
real (kind=dbl_kind), dimension(:), allocatable :: &
l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
real (kind=dbl_kind) :: &
l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
l_H2_16O_ocn , & ! local ocean concentration of H2_16O (kg/kg)
l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg)
character(len=*),parameter :: subname='(icepack_step_therm2)'
!-----------------------------------------------------------------
! Check optional arguments and set local values
!-----------------------------------------------------------------
if (present(fiso_ocn)) then
allocate(l_fiso_ocn(size(fiso_ocn)))
l_fiso_ocn(:) = fiso_ocn(:)
else
allocate(l_fiso_ocn(1))
l_fiso_ocn(:) = c0
endif
l_HDO_ocn = c0
l_H2_16O_ocn = c0
l_H2_18O_ocn = c0
if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
!-----------------------------------------------------------------
! Let rain drain through to the ocean.
!-----------------------------------------------------------------
fresh = fresh + frain * aice
!-----------------------------------------------------------------
! Given thermodynamic growth rates, transport ice between
! thickness categories.
!-----------------------------------------------------------------
! call ice_timer_start(timer_catconv) ! category conversions
!-----------------------------------------------------------------
! Compute fractional ice area in each grid cell.
!-----------------------------------------------------------------
call aggregate_area (ncat, aicen, aice, aice0)
if (icepack_warnings_aborted(subname)) return
if (kitd == 1) then
!-----------------------------------------------------------------
! Identify grid cells with ice.
!-----------------------------------------------------------------
if (aice > puny) then
call linear_itd (ncat, hin_max, &
nilyr, nslyr, &
ntrcr, trcr_depend, &
trcr_base, &
n_trcr_strata, &
nt_strata, &
aicen_init, &
vicen_init, &
aicen, &
trcrn, &
vicen, &
vsnon, &
aice , &
aice0 , &
fpond )
if (icepack_warnings_aborted(subname)) return
endif ! aice > puny
endif ! kitd = 1
! call ice_timer_stop(timer_catconv) ! category conversions
!-----------------------------------------------------------------
! Add frazil ice growing in leads.
!-----------------------------------------------------------------
! identify ice-ocean cells
call add_new_ice (ncat, nilyr, &
nfsd, nblyr, &
n_aero, dt, &
ntrcr, nltrcr, &
hin_max, ktherm, &
aicen, trcrn, &
vicen, vsnon(1), &
aice0, aice, &
frzmlt, frazil, &
frz_onset, yday, &
update_ocn_f, &
fresh, fsalt, &
Tf, sss, &
salinz, phi_init, &
dSin0_frazil, bgrid, &
cgrid, igrid, &
nbtrcr, flux_bio, &
ocean_bio, fzsal, &
frazil_diag, l_fiso_ocn, &
l_HDO_ocn, l_H2_16O_ocn, &
l_H2_18O_ocn, &
wave_sig_ht, &
wave_spectrum, &
wavefreq, dwavefreq, &
d_afsd_latg, d_afsd_newi, &
floe_rad_c, floe_binwidth)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Melt ice laterally.
!-----------------------------------------------------------------
call lateral_melt (dt, ncat, &
nilyr, nslyr, &
n_aero, fpond, &
fresh, fsalt, &
fhocn, faero_ocn, &
l_fiso_ocn, &
rside, meltl, &
fside, sss, &
aicen, vicen, &
vsnon, trcrn, &
fzsal, flux_bio, &
nbtrcr, nblyr, &
nfsd, d_afsd_latm, &
floe_rad_c,floe_binwidth)
if (icepack_warnings_aborted(subname)) return
! Floe welding during freezing conditions
if (tr_fsd) &
call fsd_weld_thermo (ncat, nfsd, &
dt, frzmlt, &
aicen, trcrn, &
d_afsd_weld)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! For the special case of a single category, adjust the area and
! volume (assuming that half the volume change decreases the
! thickness, and the other half decreases the area).
!-----------------------------------------------------------------
!echmod: test this
if (ncat==1) &
call reduce_area (hin_max (0), &
aicen (1), vicen (1), &
aicen_init(1), vicen_init(1))
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! ITD cleanup: Rebin thickness categories if necessary, and remove
! categories with very small areas.
!-----------------------------------------------------------------
call cleanup_itd (dt, ntrcr, &
nilyr, nslyr, &
ncat, hin_max, &
aicen, trcrn(1:ntrcr,:), &
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)) then
fiso_ocn(:) = l_fiso_ocn(:)
endif
deallocate(l_fiso_ocn)
end subroutine icepack_step_therm2
!=======================================================================
end module icepack_therm_itd
!=======================================================================