!=======================================================================
module icepack_therm_mushy
use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c8, c10
use icepack_parameters, only: p01, p05, p1, p2, p5, pi, bignum, puny
use icepack_parameters, only: viscosity_dyn, rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh, gravit
use icepack_parameters, only: hs_min
use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode
use icepack_parameters, only: aspect_rapid_mode, dSdt_slow_mode, phi_c_slow_mode
use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, enthalpy_snow
use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction
use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
use icepack_mushy_physics, only: temperature_snow, temperature_mush_liquid_fraction
use icepack_mushy_physics, only: liquidus_brine_salinity_mush, liquidus_temperature_mush
use icepack_mushy_physics, only: conductivity_mush_array, conductivity_snow_array
use icepack_tracers, only: tr_pond
use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf
use icepack_therm_shared, only: ferrmax
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
implicit none
private
public :: &
temperature_changes_salinity, &
permeability
real(kind=dbl_kind), parameter :: &
dTemp_errmax = 5.0e-4_dbl_kind ! max allowed change in temperature
! between iterations
!=======================================================================
contains
!=======================================================================
subroutine temperature_changes_salinity(dt, &
nilyr, nslyr, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fswsfc, fswint, &
Sswabs, Iswabs, &
hilyr, hslyr, &
apond, hpond, &
zqin, zTin, &
zqsn, zTsn, &
zSin, &
Tsf, Tbot, &
sss, &
fsensn, flatn, &
flwoutn, fsurfn, &
fcondtop, fcondbot, &
fadvheat, snoice)
! solve the enthalpy and bulk salinity of the ice for a single column
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
real (kind=dbl_kind), intent(in) :: &
rhoa , & ! air density (kg/m^3)
flw , & ! incoming longwave radiation (W/m^2)
potT , & ! air potential temperature (K)
Qa , & ! specific humidity (kg/kg)
shcoef , & ! transfer coefficient for sensible heat
lhcoef , & ! transfer coefficient for latent heat
Tbot , & ! ice bottom surfce temperature (deg C)
sss ! sea surface salinity (PSU)
real (kind=dbl_kind), intent(inout) :: &
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
fswint ! SW absorbed in ice interior below surface (W m-2)
real (kind=dbl_kind), intent(inout) :: &
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
apond , & ! melt pond area fraction
hpond ! melt pond depth (m)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
Sswabs , & ! SW radiation absorbed in snow layers (W m-2)
Iswabs ! SW radiation absorbed in ice layers (W m-2)
real (kind=dbl_kind), intent(inout):: &
fsurfn , & ! net flux to top surface, excluding fcondtopn
fcondtop , & ! downward cond flux at top surface (W m-2)
fsensn , & ! surface downward sensible heat (W m-2)
flatn , & ! surface downward latent heat (W m-2)
flwoutn ! upward LW at surface (W m-2)
real (kind=dbl_kind), intent(out):: &
fcondbot , & ! downward cond flux at bottom surface (W m-2)
fadvheat , & ! flow of heat to ocean due to advection (W m-2)
snoice ! snow ice formation
real (kind=dbl_kind), intent(inout):: &
Tsf ! ice/snow surface temperature (C)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
zqin , & ! ice layer enthalpy (J m-3)
zTin , & ! internal ice layer temperatures
zSin , & ! internal ice layer salinities
zqsn , & ! snow layer enthalpy (J m-3)
zTsn ! internal snow layer temperatures
! local variables
real(kind=dbl_kind), dimension(1:nilyr) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
zSin0 , & ! internal ice layer salinities (ppt) at start of timestep
phi , & ! liquid fraction
km , & ! ice conductivity (W m-1 K-1)
dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
real(kind=dbl_kind), dimension(1:nilyr+1) :: &
Sbr , & ! brine salinity (ppt)
qbr ! brine enthalpy (J m-3)
real(kind=dbl_kind), dimension(0:nilyr) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(1:nslyr) :: &
zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
ks ! snow conductivity (W m-1 K-1)
real(kind=dbl_kind) :: &
Tsf0 , & ! ice/snow surface temperature (C) at start of timestep
hin , & ! ice thickness (m)
hsn , & ! snow thickness (m)
hslyr_min , & ! minimum snow layer thickness (m)
w , & ! vertical flushing Darcy velocity (m/s)
qocn , & ! ocean brine enthalpy (J m-3)
qpond , & ! melt pond brine enthalpy (J m-3)
Spond ! melt pond salinity (ppt)
real(kind=dbl_kind) :: Tmlt, Iswabs_tmp, dt_rhoi_hlyr, ci, dswabs
integer(kind=int_kind) :: &
k ! ice/snow layer index
logical(kind=log_kind) :: &
lsnow ! snow presence: T: has snow, F: no snow
character(len=*),parameter :: subname='(temperature_changes_salinity)'
fadvheat = c0
snoice = c0
Tsf0 = Tsf
zqsn0 = zqsn
zqin0 = zqin
zSin0 = zSin
Spond = c0
qpond = enthalpy_brine(c0)
hslyr_min = hs_min / real(nslyr, dbl_kind)
lsnow = (hslyr > hslyr_min)
hin = hilyr * real(nilyr,dbl_kind)
qocn = enthalpy_brine(Tbot)
if (lsnow) then
hsn = hslyr * real(nslyr,dbl_kind)
else
hsn = c0
endif
do k = 1, nilyr
phi(k) = icepack_mushy_liquid_fraction(icepack_mushy_temperature_mush(zqin(k),zSin(k)),zSin(k))
enddo ! k
! calculate vertical bulk darcy flow
call flushing_velocity(zTin, &
phi, nilyr, &
hin, hsn, &
hilyr, &
hpond, apond, &
dt, w)
if (icepack_warnings_aborted(subname)) return
! calculate quantities related to drainage
call explicit_flow_velocities(nilyr, zSin, &
zTin, Tsf, &
Tbot, q, &
dSdt, Sbr, &
qbr, dt, &
sss, qocn, &
hilyr, hin)
if (icepack_warnings_aborted(subname)) return
! calculate the conductivities
call conductivity_mush_array(nilyr, zqin0, zSin0, km)
if (icepack_warnings_aborted(subname)) return
!-----------------------------------------------------------------
! Check for excessive absorbed solar radiation that may result in
! temperature overshoots. Convergence is particularly difficult
! if the starting temperature is already very close to the melting
! temperature and extra energy is added. In that case, or if the
! amount of energy absorbed is greater than the amount needed to
! melt through a given fraction of a layer, we put the extra
! energy into the surface.
! NOTE: This option is not available if the atmosphere model
! has already computed fsurf. (Unless we adjust fsurf here)
!-----------------------------------------------------------------
!mclaren: Should there be an if calc_Tsfc statement here then??
if (sw_redist) then
dt_rhoi_hlyr = dt / (rhoi*hilyr)
do k = 1, nilyr
Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
Tmlt = liquidus_temperature_mush(zSin(k))
if (zTin(k) <= Tmlt - sw_dtemp) then
ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2)
Iswabs_tmp = min(Iswabs(k), &
sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr)
endif
if (Iswabs_tmp < puny) Iswabs_tmp = c0
dswabs = min(Iswabs(k) - Iswabs_tmp, fswint)
fswsfc = fswsfc + dswabs
fswint = fswint - dswabs
Iswabs(k) = Iswabs_tmp
enddo
endif
if (lsnow) then
! case with snow
! calculate the snow conductivities
call conductivity_snow_array(ks)
if (icepack_warnings_aborted(subname)) return
! run the two stage solver
call two_stage_solver_snow(nilyr, nslyr, &
Tsf, Tsf0, &
zqsn, zqsn0, &
zqin, zqin0, &
zSin, zSin0, &
zTsn, &
zTin, &
phi, Tbot, &
km, ks, &
q, dSdt, &
w, dt, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
Iswabs, Sswabs, &
qpond, qocn, &
Spond, sss, &
hilyr, hslyr, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn )
if (icepack_warnings_aborted(subname)) then
write(warnstr,*) subname, "temperature_changes_salinity: Picard solver non-convergence (snow)"
call icepack_warnings_add(warnstr)
return
endif
! given the updated enthalpy and bulk salinity calculate other quantities
do k = 1, nslyr
zTsn(k) = temperature_snow(zqsn(k))
enddo ! k
do k = 1, nilyr
zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
qbr(k) = enthalpy_brine(zTin(k))
enddo ! k
else
! case without snow
! run the two stage solver
call two_stage_solver_nosnow(nilyr, nslyr, &
Tsf, Tsf0, &
zqsn, &
zqin, zqin0, &
zSin, zSin0, &
zTsn, &
zTin, &
phi, Tbot, &
km, ks, &
q, dSdt, &
w, dt, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
Iswabs, Sswabs, &
qpond, qocn, &
Spond, sss, &
hilyr, hslyr, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn )
if (icepack_warnings_aborted(subname)) then
write(warnstr,*) subname, "temperature_changes_salinity: Picard solver non-convergence (no snow)"
call icepack_warnings_add(warnstr)
return
endif
! given the updated enthalpy and bulk salinity calculate other quantities
do k = 1, nilyr
zTin(k) = temperature_mush_liquid_fraction(zqin(k), phi(k))
Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
qbr(k) = enthalpy_brine(zTin(k))
enddo ! k
endif
! drain ponds from flushing
call flush_pond(w, hpond, apond, dt)
if (icepack_warnings_aborted(subname)) return
! flood snow ice
call flood_ice(hsn, hin, &
nslyr, nilyr, &
hslyr, hilyr, &
zqsn, zqin, &
phi, dt, &
zSin, Sbr, &
sss, qocn, &
snoice, fadvheat)
if (icepack_warnings_aborted(subname)) return
end subroutine temperature_changes_salinity
!=======================================================================
subroutine two_stage_solver_snow(nilyr, nslyr, &
Tsf, Tsf0, &
zqsn, zqsn0, &
zqin, zqin0, &
zSin, zSin0, &
zTsn, &
zTin, &
phi, Tbot, &
km, ks, &
q, dSdt, &
w, dt, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
Iswabs, Sswabs, &
qpond, qocn, &
Spond, sss, &
hilyr, hslyr, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn )
! solve the vertical temperature and salt change for case with snow
! 1) determine what type of surface condition existed previously - cold or melting
! 2) solve the system assuming this condition persists
! 3) check the consistency of the surface condition of the solution
! 4) If the surface condition is inconsistent resolve for the other surface condition
! 5) If neither solution is consistent the resolve the inconsistency
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), intent(inout) :: &
Tsf ! snow surface temperature (C)
real(kind=dbl_kind), intent(out) :: &
fcondtop , & ! downward cond flux at top surface (W m-2)
fcondbot , & ! downward cond flux at bottom surface (W m-2)
flwoutn , & ! upward LW at surface (W m-2)
fsensn , & ! surface downward sensible heat (W m-2)
flatn , & ! surface downward latent heat (W m-2)
fsurfn , & ! net flux to top surface, excluding fcondtop
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), intent(in) :: &
Tsf0 ! snow surface temperature (C) at beginning of timestep
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqsn , & ! snow layer enthalpy (J m-3)
zTsn ! snow layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqsn0 , & ! snow layer enthalpy (J m-3) at beginning of timestep
ks , & ! snow conductivity (W m-1 K-1)
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqin , & ! ice layer enthalpy (J m-3)
zSin , & ! ice layer bulk salinity (ppt)
zTin , & ! ice layer temperature (C)
phi ! ice layer liquid fraction
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beginning of timestep
zSin0 , & ! ice layer bulk salinity (ppt) at beginning of timestep
km , & ! ice conductivity (W m-1 K-1)
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
Tbot , & ! ice bottom surfce temperature (deg C)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
fswint , & ! SW absorbed in ice interior below surface (W m-2)
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
rhoa , & ! air density (kg/m^3)
flw , & ! incoming longwave radiation (W/m^2)
potT , & ! air potential temperature (K)
Qa , & ! specific humidity (kg/kg)
shcoef , & ! transfer coefficient for sensible heat
lhcoef , & ! transfer coefficient for latent heat
w , & ! vertical flushing Darcy velocity (m/s)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
Spond , & ! melt pond salinity (ppt)
sss ! sea surface salinity (PSU)
real(kind=dbl_kind) :: &
fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
Tsf1 ! first stage ice surface temperature (C)
character(len=*),parameter :: subname='(two_stage_solver_snow)'
! determine if surface is initially cold or melting
if (Tsf < c0) then
! initially cold
! solve the system for cold and snow
call picard_solver(nilyr, nslyr, &
.true., .true., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent - surface should still be cold
if (Tsf < dTemp_errmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent - surface is warmer than melting
! resolve assuming surface is melting
Tsf1 = Tsf
! reset the solution to initial values
Tsf = c0
zqsn = zqsn0
zqin = zqin0
zSin = zSin0
! solve the system for melting and snow
call picard_solver(nilyr, nslyr, &
.true., .false., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent
! surface conductive heat flux should be less than
! incoming surface heat flux
if (fcondtop - fsurfn < ferrmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent
call two_stage_inconsistency(1, Tsf1, c0, fcondtop, fsurfn)
if (icepack_warnings_aborted(subname)) return
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" two_stage_solver_snow: two_stage_inconsistency: cold" )
return
endif ! surface flux consistency
endif ! surface temperature consistency
else
! initially melting
Tsf = c0
! solve the system for melting and snow
call picard_solver(nilyr, nslyr, &
.true., .false., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent
! surface conductive heat flux should be less than
! incoming surface heat flux
if (fcondtop - fsurfn < ferrmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent - resolve assuming other surface condition
! assume surface is cold
fcondtop1 = fcondtop
fsurfn1 = fsurfn
! reset the solution to initial values
Tsf = Tsf0
zqsn = zqsn0
zqin = zqin0
zSin = zSin0
! solve the system for cold and snow
call picard_solver(nilyr, nslyr, &
.true., .true., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent - surface should be cold
if (Tsf < dTemp_errmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent
! failed to find a solution so need to refine solutions until consistency found
call two_stage_inconsistency(2, Tsf, c0, fcondtop1, fsurfn1)
if (icepack_warnings_aborted(subname)) return
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" two_stage_solver_snow: two_stage_inconsistency: melting" )
return
endif ! surface temperature consistency
endif ! surface flux consistency
endif
end subroutine two_stage_solver_snow
!=======================================================================
subroutine two_stage_solver_nosnow(nilyr, nslyr, &
Tsf, Tsf0, &
zqsn, &
zqin, zqin0, &
zSin, zSin0, &
zTsn, &
zTin, &
phi, Tbot, &
km, ks, &
q, dSdt, &
w, dt, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
Iswabs, Sswabs, &
qpond, qocn, &
Spond, sss, &
hilyr, hslyr, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn )
! solve the vertical temperature and salt change for case with no snow
! 1) determine what type of surface condition existed previously - cold or melting
! 2) solve the system assuming this condition persists
! 3) check the consistency of the surface condition of the solution
! 4) If the surface condition is inconsistent resolve for the other surface condition
! 5) If neither solution is consistent the resolve the inconsistency
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), intent(inout) :: &
Tsf ! ice surface temperature (C)
real(kind=dbl_kind), intent(out) :: &
fcondtop , & ! downward cond flux at top surface (W m-2)
fcondbot , & ! downward cond flux at bottom surface (W m-2)
flwoutn , & ! upward LW at surface (W m-2)
fsensn , & ! surface downward sensible heat (W m-2)
flatn , & ! surface downward latent heat (W m-2)
fsurfn , & ! net flux to top surface, excluding fcondtop
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), intent(in) :: &
Tsf0 ! ice surface temperature (C) at beginning of timestep
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqsn , & ! snow layer enthalpy (J m-3)
zTsn ! snow layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(in) :: &
ks , & ! snow conductivity (W m-1 K-1)
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqin , & ! ice layer enthalpy (J m-3)
zSin , & ! ice layer bulk salinity (ppt)
zTin , & ! ice layer temperature (C)
phi ! ice layer liquid fraction
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beginning of timestep
zSin0 , & ! ice layer bulk salinity (ppt) at beginning of timestep
km , & ! ice conductivity (W m-1 K-1)
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
fswint , & ! SW absorbed in ice interior below surface (W m-2)
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
rhoa , & ! air density (kg/m^3)
flw , & ! incoming longwave radiation (W/m^2)
potT , & ! air potential temperature (K)
Qa , & ! specific humidity (kg/kg)
shcoef , & ! transfer coefficient for sensible heat
lhcoef , & ! transfer coefficient for latent heat
w , & ! vertical flushing Darcy velocity (m/s)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
Spond , & ! melt pond salinity (ppt)
sss ! sea surface salinity (PSU)
real(kind=dbl_kind) :: &
Tmlt , & ! upper ice layer melting temperature (C)
fcondtop1 , & ! first stage downward cond flux at top surface (W m-2)
fsurfn1 , & ! first stage net flux to top surface, excluding fcondtop
Tsf1 ! first stage ice surface temperature (C)
character(len=*),parameter :: subname='(two_stage_solver_nosnow)'
! initial surface melting temperature
Tmlt = liquidus_temperature_mush(zSin0(1))
! determine if surface is initially cold or melting
if (Tsf < Tmlt) then
! initially cold
! solve the system for cold and no snow
call picard_solver(nilyr, nslyr, &
.false., .true., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent - surface should still be cold
if (Tsf < Tmlt + dTemp_errmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent - surface is warmer than melting
! resolve assuming surface is melting
Tsf1 = Tsf
! reset the solution to initial values
Tsf = liquidus_temperature_mush(zSin0(1))
zqin = zqin0
zSin = zSin0
! solve the system for melt and no snow
call picard_solver(nilyr, nslyr, &
.false., .false., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent
! surface conductive heat flux should be less than
! incoming surface heat flux
if (fcondtop - fsurfn < ferrmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent
call two_stage_inconsistency(3, Tsf1, Tmlt, fcondtop, fsurfn)
if (icepack_warnings_aborted(subname)) return
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" two_stage_solver_nosnow: two_stage_inconsistency: cold" )
return
endif
endif
else
! initially melting
! solve the system for melt and no snow
Tsf = Tmlt
call picard_solver(nilyr, nslyr, &
.false., .false., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent
! surface conductive heat flux should be less than
! incoming surface heat flux
if (fcondtop - fsurfn < ferrmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent - resolve assuming other surface condition
! assume surface is cold
fcondtop1 = fcondtop
fsurfn1 = fsurfn
! reset the solution to initial values
Tsf = Tsf0
zqin = zqin0
zSin = zSin0
! solve the system for cold and no snow
call picard_solver(nilyr, nslyr, &
.false., .true., &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
if (icepack_warnings_aborted(subname)) return
! halt if solver failed
if (icepack_warnings_aborted(subname)) return
! check if solution is consistent - surface should be cold
if (Tsf < Tmlt + dTemp_errmax) then
! solution is consistent - have solution so finish
return
else
! solution is inconsistent
call two_stage_inconsistency(4, Tsf, Tmlt, fcondtop1, fsurfn1)
if (icepack_warnings_aborted(subname)) return
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//" two_stage_solver_nosnow: two_stage_inconsistency: melting" )
return
endif
endif
endif
end subroutine two_stage_solver_nosnow
!=======================================================================
subroutine two_stage_inconsistency(type, Tsf, Tmlt, fcondtop, fsurfn)
integer (kind=int_kind), intent(in) :: &
type
real(kind=dbl_kind), intent(in) :: &
Tsf, &
Tmlt, &
fcondtop, &
fsurfn
character(len=*),parameter :: subname='(two_stage_inconsistency)'
write(warnstr,*) subname, "icepack_therm_mushy: two stage inconsistency"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "type:", type
call icepack_warnings_add(warnstr)
if (type == 1) then
write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
call icepack_warnings_add(warnstr)
else if (type == 2) then
write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
call icepack_warnings_add(warnstr)
else if (type == 3) then
write(warnstr,*) subname, "First stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "Second stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
call icepack_warnings_add(warnstr)
else if (type == 4) then
write(warnstr,*) subname, "First stage : fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", fcondtop, fsurfn, ferrmax, fcondtop - fsurfn - ferrmax
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "Second stage : Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, " :", Tsf, Tmlt, dTemp_errmax, Tsf - Tmlt - dTemp_errmax
call icepack_warnings_add(warnstr)
endif
end subroutine two_stage_inconsistency
!=======================================================================
! Picard/TDMA based solver
!=======================================================================
subroutine prep_picard(nilyr, nslyr, &
lsnow, zqsn, &
zqin, zSin, &
hilyr, hslyr, &
km, ks, &
zTin, zTsn, &
Sbr, phi, &
dxp, kcstar, &
einit)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin , & ! ice layer enthalpy (J m-3)
zSin , & ! ice layer bulk salinity (ppt)
km , & ! ice conductivity (W m-1 K-1)
zqsn , & ! snow layer enthalpy (J m-3)
ks ! snow conductivity (W m-1 K-1)
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hslyr ! snow layer thickness (m)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zTin , & ! ice layer temperature (C)
Sbr , & ! ice layer brine salinity (ppt)
phi , & ! ice layer liquid fraction
zTsn , & ! snow layer temperature (C)
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
real(kind=dbl_kind), intent(out) :: &
einit ! initial total energy (J)
integer(kind=int_kind) :: k
character(len=*),parameter :: subname='(prep_picard)'
! calculate initial ice temperatures
do k = 1, nilyr
zTin(k) = icepack_mushy_temperature_mush(zqin(k), zSin(k))
Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
enddo ! k
if (lsnow) then
do k = 1, nslyr
zTsn(k) = temperature_snow(zqsn(k))
enddo ! k
endif ! lsnow
! interface distances
call calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
if (icepack_warnings_aborted(subname)) return
! interface conductivities
call calc_intercell_conductivity(lsnow, nilyr, nslyr, &
km, ks, hilyr, hslyr, kcstar)
if (icepack_warnings_aborted(subname)) return
! total energy content
call total_energy_content(lsnow, &
nilyr, nslyr, &
zqin, zqsn, &
hilyr, hslyr, &
einit)
if (icepack_warnings_aborted(subname)) return
end subroutine prep_picard
!=======================================================================
subroutine picard_solver(nilyr, nslyr, &
lsnow, lcold, &
Tsf, zqsn, &
zqin, zSin, &
zTin, zTsn, &
phi, dt, &
hilyr, hslyr, &
km, ks, &
Iswabs, Sswabs, &
Tbot, &
fswint, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
fcondtop, fcondbot, &
fadvheat, &
flwoutn, fsensn, &
flatn, fsurfn, &
qpond, qocn, &
Spond, sss, &
q, dSdt, &
w )
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical, intent(in) :: &
lsnow , & ! snow presence: T: has snow, F: no snow
lcold ! surface cold: T: surface is cold, F: surface is melting
real(kind=dbl_kind), intent(inout) :: &
Tsf ! snow surface temperature (C)
real(kind=dbl_kind), intent(out) :: &
fcondtop , & ! downward cond flux at top surface (W m-2)
fcondbot , & ! downward cond flux at bottom surface (W m-2)
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqin , & ! ice layer enthalpy (J m-3)
zSin , & ! ice layer bulk salinity (ppt)
zTin , & ! ice layer temperature (C)
phi , & ! ice layer liquid fraction
zqsn , & ! snow layer enthalpy (J m-3)
zTsn ! snow layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(in) :: &
km , & ! ice conductivity (W m-1 K-1)
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
ks , & ! snow conductivity (W m-1 K-1)
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), intent(out) :: &
flwoutn , & ! upward LW at surface (W m-2)
fsensn , & ! surface downward sensible heat (W m-2)
flatn , & ! surface downward latent heat (W m-2)
fsurfn ! net flux to top surface, excluding fcondtop
real(kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
fswint , & ! SW absorbed in ice interior below surface (W m-2)
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
rhoa , & ! air density (kg/m^3)
flw , & ! incoming longwave radiation (W/m^2)
potT , & ! air potential temperature (K)
Qa , & ! specific humidity (kg/kg)
shcoef , & ! transfer coefficient for sensible heat
lhcoef , & ! transfer coefficient for latent heat
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
Spond , & ! melt pond salinity (ppt)
sss , & ! sea surface salinity (ppt)
w ! vertical flushing Darcy velocity (m/s)
real(kind=dbl_kind), dimension(nilyr) :: &
Sbr , & ! ice layer brine salinity (ppt)
qbr , & ! ice layer brine enthalpy (J m-3)
zTin0 , & ! ice layer temperature (C) at start of timestep
zqin0 , & ! ice layer enthalpy (J m-3) at start of timestep
zSin0 , & ! ice layer bulk salinity (ppt) at start of timestep
zTin_prev ! ice layer temperature at previous iteration
real(kind=dbl_kind), dimension(nslyr) :: &
zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
zTsn0 , & ! snow layer temperature (C) at start of timestep
zTsn_prev ! snow layer temperature at previous iteration
real(kind=dbl_kind), dimension(nslyr+nilyr+1) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
real(kind=dbl_kind) :: &
Tsf0 , & ! snow surface temperature (C) at start of timestep
dfsurfn_dTsf , & ! derivative of net flux to top surface, excluding fcondtopn
dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature
dfsensn_dTsf , & ! derivative of sensible heat flux wrt surface temperature
dflatn_dTsf , & ! derivative of latent heat flux wrt surface temperature
Tsf_prev , & ! snow surface temperature at previous iteration
einit , & ! initial total energy (J)
fadvheat_nit ! heat to ocean due to advection (W m-2) during iteration
logical :: &
lconverged ! has Picard solver converged?
integer :: &
nit ! Picard iteration count
integer, parameter :: &
nit_max = 100 ! maximum number of Picard iterations
character(len=*),parameter :: subname='(picard_solver)'
lconverged = .false.
! prepare quantities for picard iteration
call prep_picard(nilyr, nslyr, &
lsnow, zqsn, &
zqin, zSin, &
hilyr, hslyr, &
km, ks, &
zTin, zTsn, &
Sbr, phi, &
dxp, kcstar, &
einit)
if (icepack_warnings_aborted(subname)) return
Tsf0 = Tsf
zqin0 = zqin
zqsn0 = zqsn
zTin0 = zTin
zTsn0 = zTsn
zSin0 = zSin
! set prev variables
Tsf_prev = Tsf
zTsn_prev = zTsn
zTin_prev = zTin
! picard iteration
picard: do nit = 1, nit_max
! surface heat flux
call surface_heat_flux(Tsf, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
flwoutn, fsensn, &
flatn, fsurfn)
if (icepack_warnings_aborted(subname)) return
! derivative of heat flux with respect to surface temperature
call dsurface_heat_flux_dTsf(Tsf, rhoa, &
shcoef, lhcoef, &
dfsurfn_dTsf, dflwoutn_dTsf, &
dfsensn_dTsf, dflatn_dTsf)
if (icepack_warnings_aborted(subname)) return
! tridiagonal solve of new temperatures
call solve_heat_conduction(lsnow, lcold, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
phi, dt, &
qpond, qocn, &
q, w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
fsurfn, dfsurfn_dTsf, &
zTin, zTsn)
if (icepack_warnings_aborted(subname)) return
! update brine enthalpy
call picard_updates_enthalpy(nilyr, zTin, qbr)
if (icepack_warnings_aborted(subname)) return
! drainage fluxes
call picard_drainage_fluxes(fadvheat_nit, q, &
qbr, qocn, &
nilyr)
if (icepack_warnings_aborted(subname)) return
! flushing fluxes
call picard_flushing_fluxes(nilyr, &
fadvheat_nit, w, &
qbr, &
qpond)
if (icepack_warnings_aborted(subname)) return
! perform convergence check
call check_picard_convergence(nilyr, nslyr, &
lsnow, &
lconverged, &
Tsf, Tsf_prev, &
zTin, zTin_prev,&
zTsn, zTsn_prev,&
phi, Tbot, &
zqin, zqsn, &
km, ks, &
hilyr, hslyr, &
fswint, &
einit, dt, &
fcondtop, fcondbot, &
fadvheat_nit)
if (icepack_warnings_aborted(subname)) return
if (lconverged) exit
Tsf_prev = Tsf
zTsn_prev = zTsn
zTin_prev = zTin
enddo picard
fadvheat = fadvheat_nit
! update the picard iterants
call picard_updates(nilyr, zTin, &
Sbr, qbr)
if (icepack_warnings_aborted(subname)) return
! solve for the salinity
call solve_salinity(zSin, Sbr, &
Spond, sss, &
q, dSdt, &
w, hilyr, &
dt, nilyr)
if (icepack_warnings_aborted(subname)) return
! final surface heat flux
call surface_heat_flux(Tsf, fswsfc, &
rhoa, flw, &
potT, Qa, &
shcoef, lhcoef, &
flwoutn, fsensn, &
flatn, fsurfn)
if (icepack_warnings_aborted(subname)) return
! if not converged
if (.not. lconverged) then
call picard_nonconvergence(nilyr, nslyr,&
Tsf0, Tsf, &
zTsn0, zTsn, &
zTin0, zTin, &
zSin0, zSin, &
zqsn0, zqsn, &
zqin0, zqin, &
phi)
if (icepack_warnings_aborted(subname)) return
call icepack_warnings_add(subname//" picard_solver: Picard solver non-convergence" )
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
if (icepack_warnings_aborted(subname)) return
endif
end subroutine picard_solver
!=======================================================================
subroutine picard_nonconvergence(nilyr, nslyr,&
Tsf0, Tsf, &
zTsn0, zTsn, &
zTin0, zTin, &
zSin0, zSin, &
zqsn0, zqsn, &
zqin0, zqin, phi)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), intent(in) :: &
Tsf0 , & ! snow surface temperature (C) at beginning of timestep
Tsf ! snow surface temperature (C)
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTsn0 , & ! snow layer temperature (C) at beginning of timestep
zTsn , & ! snow layer temperature (C)
zqsn0 , &
zqsn , &
zTin0 , & ! ice layer temperature (C)
zTin , & ! ice layer temperature (C)
zSin0 , & ! ice layer bulk salinity (ppt)
zSin , & ! ice layer bulk salinity (ppt)
zqin0 , &
zqin , &
phi ! ice layer liquid fraction
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(picard_nonconvergence)'
write(warnstr,*) subname, "-------------------------------------"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, "picard convergence failed!"
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 0, Tsf0, Tsf
call icepack_warnings_add(warnstr)
do k = 1, nslyr
write(warnstr,*) subname, k, zTsn0(k), zTsn(k), zqsn(k), zqsn0(k)
call icepack_warnings_add(warnstr)
enddo ! k
do k = 1, nilyr
write(warnstr,*) subname, k, zTin0(k), zTin(k), zSin0(k), zSin(k), phi(k), zqin(k), zqin0(k)
call icepack_warnings_add(warnstr)
enddo ! k
write(warnstr,*) subname, "-------------------------------------"
call icepack_warnings_add(warnstr)
end subroutine picard_nonconvergence
!=======================================================================
subroutine check_picard_convergence(nilyr, nslyr, &
lsnow, &
lconverged, &
Tsf, Tsf_prev, &
zTin, zTin_prev,&
zTsn, zTsn_prev,&
phi, Tbot, &
zqin, zqsn, &
km, ks, &
hilyr, hslyr, &
fswint, &
einit, dt, &
fcondtop, fcondbot, &
fadvheat)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical, intent(inout) :: &
lconverged ! has Picard solver converged?
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
Tsf , & ! snow surface temperature (C)
Tsf_prev , & ! snow surface temperature at previous iteration
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
fswint , & ! SW absorbed in ice interior below surface (W m-2)
einit , & ! initial total energy (J)
Tbot , & ! ice bottom surfce temperature (deg C)
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTin , & ! ice layer temperature (C)
zTin_prev, & ! ice layer temperature at previous iteration
phi , & ! ice layer liquid fraction
km ! ice conductivity (W m-1 K-1)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqin ! ice layer enthalpy (J m-3)
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTsn , & ! snow layer temperature (C)
zTsn_prev, & ! snow layer temperature at previous iteration
ks ! snow conductivity (W m-1 K-1)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqsn ! snow layer enthalpy (J m-3)
real(kind=dbl_kind), intent(out) :: &
fcondtop , & ! downward cond flux at top surface (W m-2)
fcondbot ! downward cond flux at bottom surface (W m-2)
real(kind=dbl_kind) :: &
ferr , & ! energy flux error
efinal , & ! initial total energy (J) at iteration
dzTsn , & ! change in snow temperature (C) between iterations
dzTin , & ! change in ice temperature (C) between iterations
dTsf ! change in surface temperature (C) between iterations
character(len=*),parameter :: subname='(check_picard_convergence)'
call picard_final(lsnow, &
nilyr, nslyr, &
zqin, zqsn, &
zTin, zTsn, &
phi)
if (icepack_warnings_aborted(subname)) return
call total_energy_content(lsnow, &
nilyr, nslyr, &
zqin, zqsn, &
hilyr, hslyr, &
efinal)
if (icepack_warnings_aborted(subname)) return
call maximum_variables_changes(lsnow, &
Tsf, Tsf_prev, dTsf, &
zTsn, zTsn_prev, dzTsn, &
zTin, zTin_prev, dzTin)
if (icepack_warnings_aborted(subname)) return
fcondbot = c2 * km(nilyr) * ((zTin(nilyr) - Tbot) / hilyr)
if (lsnow) then
fcondtop = c2 * ks(1) * ((Tsf - zTsn(1)) / hslyr)
else
fcondtop = c2 * km(1) * ((Tsf - zTin(1)) / hilyr)
endif
ferr = (efinal - einit) / dt - (fcondtop - fcondbot + fswint - fadvheat)
lconverged = (dTsf < dTemp_errmax .and. &
dzTsn < dTemp_errmax .and. &
dzTin < dTemp_errmax .and. &
abs(ferr) < 0.9_dbl_kind*ferrmax)
end subroutine check_picard_convergence
!=======================================================================
subroutine picard_drainage_fluxes(fadvheat, q, &
qbr, qocn, &
nilyr)
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), intent(out) :: &
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
qbr ! ice layer brine enthalpy (J m-3)
real(kind=dbl_kind), intent(in) :: &
qocn ! ocean brine enthalpy (J m-3)
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(picard_drainage_fluxes)'
fadvheat = c0
! calculate fluxes from base upwards
do k = 1, nilyr-1
fadvheat = fadvheat - q(k) * (qbr(k+1) - qbr(k))
enddo ! k
k = nilyr
fadvheat = fadvheat - q(k) * (qocn - qbr(k))
end subroutine picard_drainage_fluxes
!=======================================================================
subroutine picard_flushing_fluxes(nilyr, &
fadvheat, w, &
qbr, &
qpond)
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), intent(inout) :: &
fadvheat ! flow of heat to ocean due to advection (W m-2)
real(kind=dbl_kind), dimension(:), intent(in) :: &
qbr ! ice layer brine enthalpy (J m-3)
real(kind=dbl_kind), intent(in) :: &
w , & ! vertical flushing Darcy velocity (m/s)
qpond ! melt pond brine enthalpy (J m-3)
character(len=*),parameter :: subname='(picard_flushing_fluxes)'
fadvheat = fadvheat + w * (qbr(nilyr) - qpond)
end subroutine picard_flushing_fluxes
!=======================================================================
subroutine maximum_variables_changes(lsnow, &
Tsf, Tsf_prev, dTsf, &
zTsn, zTsn_prev, dzTsn, &
zTin, zTin_prev, dzTin)
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! snow surface temperature (C)
Tsf_prev ! snow surface temperature at previous iteration
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTsn , & ! snow layer temperature (C)
zTsn_prev, & ! snow layer temperature at previous iteration
zTin , & ! ice layer temperature (C)
zTin_prev ! ice layer temperature at previous iteration
real(kind=dbl_kind), intent(out) :: &
dTsf , & ! change in surface temperature (C) between iterations
dzTsn , & ! change in snow temperature (C) between iterations
dzTin ! change in surface temperature (C) between iterations
character(len=*),parameter :: subname='(maximum_variables_changes)'
dTsf = abs(Tsf - Tsf_prev)
if (lsnow) then
dzTsn = maxval(abs(zTsn - zTsn_prev))
else ! lsnow
dzTsn = c0
endif ! lsnow
dzTin = maxval(abs(zTin - zTin_prev))
end subroutine maximum_variables_changes
!=======================================================================
subroutine total_energy_content(lsnow, &
nilyr, nslyr, &
zqin, zqsn, &
hilyr, hslyr, &
energy)
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin , & ! ice layer enthalpy (J m-3)
zqsn ! snow layer enthalpy (J m-3)
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hslyr ! snow layer thickness (m)
real(kind=dbl_kind), intent(out) :: &
energy ! total energy of ice and snow
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(total_energy_content)'
energy = c0
if (lsnow) then
do k = 1, nslyr
energy = energy + hslyr * zqsn(k)
enddo ! k
endif ! lsnow
do k = 1, nilyr
energy = energy + hilyr * zqin(k)
enddo ! k
end subroutine total_energy_content
!=======================================================================
subroutine picard_updates(nilyr, zTin, &
Sbr, qbr)
! update brine salinity and liquid fraction based on new temperatures
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTin ! ice layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
Sbr , & ! ice layer brine salinity (ppt)
qbr ! ice layer brine enthalpy (J m-3)
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(picard_updates)'
do k = 1, nilyr
Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
qbr(k) = enthalpy_brine(zTin(k))
enddo ! k
end subroutine picard_updates
!=======================================================================
subroutine picard_updates_enthalpy(nilyr, zTin, qbr)
! update brine salinity and liquid fraction based on new temperatures
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTin ! ice layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
qbr ! ice layer brine enthalpy (J m-3)
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(picard_updates_enthalpy)'
do k = 1, nilyr
qbr(k) = enthalpy_brine(zTin(k))
enddo ! k
end subroutine picard_updates_enthalpy
!=======================================================================
subroutine picard_final(lsnow, &
nilyr, nslyr, &
zqin, zqsn, &
zTin, zTsn, &
phi)
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nslyr ! number of snow layers
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqin, & ! ice layer enthalpy (J m-3)
zqsn ! snow layer enthalpy (J m-3)
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTin, & ! ice layer temperature (C)
phi , & ! ice layer liquid fraction
zTsn ! snow layer temperature (C)
integer :: &
k ! vertical layer index
character(len=*),parameter :: subname='(picard_final)'
do k = 1, nilyr
zqin(k) = enthalpy_mush_liquid_fraction(zTin(k), phi(k))
enddo ! k
if (lsnow) then
do k = 1, nslyr
zqsn(k) = enthalpy_snow(zTsn(k))
enddo ! k
endif ! lsnow
end subroutine picard_final
!=======================================================================
subroutine calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nslyr ! number of snow layers
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hslyr ! snow layer thickness (m)
real(kind=dbl_kind), dimension(:), intent(out) :: &
dxp ! distances between grid points (m)
integer :: &
l ! vertical index
character(len=*),parameter :: subname='(calc_intercell_thickness)'
if (lsnow) then
dxp(1) = hslyr / c2
do l = 2, nslyr
dxp(l) = hslyr
enddo ! l
dxp(nslyr+1) = (hilyr + hslyr) / c2
do l = nslyr+2, nilyr+nslyr
dxp(l) = hilyr
enddo ! l
dxp(nilyr+nslyr+1) = hilyr / c2
else ! lsnow
dxp(1) = hilyr / c2
do l = 2, nilyr
dxp(l) = hilyr
enddo ! l
dxp(nilyr+1) = hilyr / c2
do l = nilyr+2, nilyr+nslyr+1
dxp(l) = c0
enddo ! l
endif ! lsnow
end subroutine calc_intercell_thickness
!=======================================================================
subroutine calc_intercell_conductivity(lsnow, &
nilyr, nslyr, &
km, ks, &
hilyr, hslyr, &
kcstar)
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nslyr ! number of snow layers
logical, intent(in) :: &
lsnow ! snow presence: T: has snow, F: no snow
real(kind=dbl_kind), dimension(:), intent(in) :: &
km , & ! ice conductivity (W m-1 K-1)
ks ! snow conductivity (W m-1 K-1)
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hslyr ! snow layer thickness (m)
real(kind=dbl_kind), dimension(:), intent(out) :: &
kcstar ! interface conductivities (W m-1 K-1)
real(kind=dbl_kind) :: &
fe ! distance fraction at interface
integer :: &
k, & ! vertical layer index
l ! vertical index
character(len=*),parameter :: subname='(calc_intercell_conductivity)'
if (lsnow) then
kcstar(1) = ks(1)
do l = 2, nslyr
k = l
kcstar(l) = (c2 * ks(k) * ks(k-1)) / (ks(k) + ks(k-1))
enddo ! l
fe = hilyr / (hilyr + hslyr)
kcstar(nslyr+1) = c1 / ((c1 - fe) / ks(nslyr) + fe / km(1))
do k = 2, nilyr
l = k + nslyr
kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
enddo ! k
kcstar(nilyr+nslyr+1) = km(nilyr)
else ! lsnow
kcstar(1) = km(1)
do k = 2, nilyr
l = k
kcstar(l) = (c2 * km(k) * km(k-1)) / (km(k) + km(k-1))
enddo ! k
kcstar(nilyr+1) = km(nilyr)
do l = nilyr+2, nilyr+nslyr+1
kcstar(l) = c0
enddo ! l
endif ! lsnow
end subroutine calc_intercell_conductivity
!=======================================================================
subroutine solve_heat_conduction(lsnow, lcold, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
phi, dt, &
qpond, qocn, &
q, w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
fsurfn, dfsurfn_dTsf, &
zTin, zTsn)
logical, intent(in) :: &
lsnow , & ! snow presence: T: has snow, F: no snow
lcold ! surface cold: T: surface is cold, F: surface is melting
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
phi , & ! ice layer liquid fraction
zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), intent(inout) :: &
Tsf ! snow surface temperature (C)
real(kind=dbl_kind), intent(in) :: &
dt , & ! timestep (s)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
w , & ! vertical flushing Darcy velocity (m/s)
fsurfn , & ! net flux to top surface, excluding fcondtop
dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zTin ! ice layer temperature (C)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zTsn ! snow layer temperature (C)
real(kind=dbl_kind), dimension(nilyr+nslyr+1) :: &
Ap , & ! diagonal of tridiagonal matrix
As , & ! lower off-diagonal of tridiagonal matrix
An , & ! upper off-diagonal of tridiagonal matrix
b , & ! right hand side of matrix solve
T ! ice and snow temperatures
integer :: &
nyn ! matrix size
character(len=*),parameter :: subname='(solve_heat_conduction)'
! set up matrix and right hand side - snow
if (lsnow) then
if (lcold) then
call matrix_elements_snow_cold(Ap, As, An, b, nyn, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
fsurfn, dfsurfn_dTsf, &
dt)
if (icepack_warnings_aborted(subname)) return
else ! lcold
call matrix_elements_snow_melt(Ap, As, An, b, nyn, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
dt)
if (icepack_warnings_aborted(subname)) return
endif ! lcold
else ! lsnow
if (lcold) then
call matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
nilyr, &
Tsf, Tbot, &
zqin0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, &
dxp, kcstar, &
Iswabs, &
fsurfn, dfsurfn_dTsf, &
dt)
if (icepack_warnings_aborted(subname)) return
else ! lcold
call matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
nilyr, &
Tsf, Tbot, &
zqin0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, &
dxp, kcstar, &
Iswabs, &
dt)
if (icepack_warnings_aborted(subname)) return
endif ! lcold
endif ! lsnow
! tridiag to get new temperatures
call tdma_solve_sparse(nilyr, nslyr, &
An(1:nyn), Ap(1:nyn), As(1:nyn), b(1:nyn), T(1:nyn), nyn)
if (icepack_warnings_aborted(subname)) return
call update_temperatures(lsnow, lcold, &
nilyr, nslyr, &
T, Tsf, &
zTin, zTsn)
if (icepack_warnings_aborted(subname)) return
end subroutine solve_heat_conduction
!=======================================================================
subroutine update_temperatures(lsnow, lcold, &
nilyr, nslyr, &
T, Tsf, &
zTin, zTsn)
logical, intent(in) :: &
lsnow , & ! snow presence: T: has snow, F: no snow
lcold ! surface cold: T: surface is cold, F: surface is melting
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
T ! matrix solution vector
real(kind=dbl_kind), intent(inout) :: &
Tsf ! snow surface temperature (C)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zTin , & ! ice layer temperature (C)
zTsn ! snow layer temperature (C)
integer :: &
l , & ! vertical index
k ! vertical layer index
character(len=*),parameter :: subname='(update_temperatures)'
if (lsnow) then
if (lcold) then
Tsf = T(1)
do k = 1, nslyr
l = k + 1
zTsn(k) = T(l)
enddo ! k
do k = 1, nilyr
l = k + nslyr + 1
zTin(k) = T(l)
enddo ! k
else ! lcold
do k = 1, nslyr
l = k
zTsn(k) = T(l)
enddo ! k
do k = 1, nilyr
l = k + nslyr
zTin(k) = T(l)
enddo ! k
endif ! lcold
else ! lsnow
if (lcold) then
Tsf = T(1)
do k = 1, nilyr
l = k + 1
zTin(k) = T(l)
enddo ! k
else ! lcold
do k = 1, nilyr
l = k
zTin(k) = T(l)
enddo ! k
endif ! lcold
endif ! lsnow
end subroutine update_temperatures
!=======================================================================
subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn, &
nilyr, &
Tsf, Tbot, &
zqin0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, &
dxp, kcstar, &
Iswabs, &
dt)
real(kind=dbl_kind), dimension(:), intent(out) :: &
Ap , & ! diagonal of tridiagonal matrix
As , & ! lower off-diagonal of tridiagonal matrix
An , & ! upper off-diagonal of tridiagonal matrix
b ! right hand side of matrix solve
integer, intent(out) :: &
nyn ! matrix size
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
phi ! ice layer liquid fraction
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! snow surface temperature (C)
dt , & ! timestep (s)
hilyr , & ! ice layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
w ! downwards vertical flushing Darcy velocity (m/s)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
integer :: &
k , & ! vertical layer index
l ! vertical index
character(len=*),parameter :: subname='(matrix_elements_nosnow_melt)'
! surface layer
k = 1
l = k
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(k+1) / dxp(k+1) - &
q(k) * cp_ocn * rhow
An(l) = c0
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
(kcstar(k) / dxp(k)) * Tsf + &
w * qpond
! interior ice layers
do k = 2, nilyr-1
l = k
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(k+1) / dxp(k+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(k) / dxp(k) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
enddo ! k
! bottom layer
k = nilyr
l = k
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = c0
An(l) = -kcstar(k) / dxp(k) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
(kcstar(k+1) * Tbot) / dxp(k+1) + &
q(k) * qocn
nyn = nilyr
end subroutine matrix_elements_nosnow_melt
!=======================================================================
subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn, &
nilyr, &
Tsf, Tbot, &
zqin0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, &
dxp, kcstar, &
Iswabs, &
fsurfn, dfsurfn_dTsf, &
dt)
real(kind=dbl_kind), dimension(:), intent(out) :: &
Ap , & ! diagonal of tridiagonal matrix
As , & ! lower off-diagonal of tridiagonal matrix
An , & ! upper off-diagonal of tridiagonal matrix
b ! right hand side of matrix solve
integer, intent(out) :: &
nyn ! matrix size
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
phi ! ice layer liquid fraction
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! snow surface temperature (C)
dt , & ! timestep (s)
hilyr , & ! ice layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
w , & ! downwards vertical flushing Darcy velocity (m/s)
fsurfn , & ! net flux to top surface, excluding fcondtop
dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
integer :: &
k , & ! vertical layer index
l ! vertical index
character(len=*),parameter :: subname='(matrix_elements_nosnow_cold)'
! surface temperature
l = 1
Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
As(l) = kcstar(1) / dxp(1)
An(l) = c0
b (l) = dfsurfn_dTsf * Tsf - fsurfn
! surface layer
k = 1
l = k + 1
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(k+1) / dxp(k+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(k) / dxp(k)
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
w * qpond
! interior ice layers
do k = 2, nilyr-1
l = k + 1
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(k+1) / dxp(k+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(k) / dxp(k) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
enddo ! k
! bottom layer
k = nilyr
l = k + 1
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(k+1) / dxp(k+1) + &
kcstar(k) / dxp(k) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = c0
An(l) = -kcstar(k) / dxp(k) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
(kcstar(k+1) * Tbot) / dxp(k+1) + &
q(k) * qocn
nyn = nilyr + 1
end subroutine matrix_elements_nosnow_cold
!=======================================================================
subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
dt)
real(kind=dbl_kind), dimension(:), intent(out) :: &
Ap , & ! diagonal of tridiagonal matrix
As , & ! lower off-diagonal of tridiagonal matrix
An , & ! upper off-diagonal of tridiagonal matrix
b ! right hand side of matrix solve
integer, intent(out) :: &
nyn ! matrix size
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
phi , & ! ice layer liquid fraction
zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! snow surface temperature (C)
dt , & ! timestep (s)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
w ! downwards vertical flushing Darcy velocity (m/s)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
integer :: &
k , & ! vertical layer index
l ! vertical index
character(len=*),parameter :: subname='(matrix_elements_snow_melt)'
! surface layer
k = 1
l = k
Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
kcstar(l+1) / dxp(l+1) + &
kcstar(l) / dxp(l)
As(l) = -kcstar(l+1) / dxp(l+1)
An(l) = c0
b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k) + &
(kcstar(l) * Tsf) / dxp(l)
! interior snow layers
do k = 2, nslyr
l = k
Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
kcstar(l+1) / dxp(l+1) + &
kcstar(l) / dxp(l)
As(l) = -kcstar(l+1) / dxp(l+1)
An(l) = -kcstar(l) / dxp(l)
b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
enddo ! k
! top ice layer
k = 1
l = nslyr + k
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(l+1) / dxp(l+1) + &
kcstar(l) / dxp(l) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(l+1) / dxp(l+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(l) / dxp(l)
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
w * qpond
! interior ice layers
do k = 2, nilyr-1
l = nslyr + k
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(l+1) / dxp(l+1) + &
kcstar(l) / dxp(l) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(l+1) / dxp(l+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(l) / dxp(l) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
enddo ! k
! bottom layer
k = nilyr
l = nilyr + nslyr
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(l+1) / dxp(l+1) + &
kcstar(l) / dxp(l) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = c0
An(l) = -kcstar(l) / dxp(l) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
(kcstar(l+1) * Tbot) / dxp(l+1) + &
q(k) * qocn
nyn = nilyr + nslyr
end subroutine matrix_elements_snow_melt
!=======================================================================
subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn, &
nilyr, nslyr, &
Tsf, Tbot, &
zqin0, zqsn0, &
qpond, qocn, &
phi, q, &
w, &
hilyr, hslyr, &
dxp, kcstar, &
Iswabs, Sswabs, &
fsurfn, dfsurfn_dTsf, &
dt)
real(kind=dbl_kind), dimension(:), intent(out) :: &
Ap , & ! diagonal of tridiagonal matrix
As , & ! lower off-diagonal of tridiagonal matrix
An , & ! upper off-diagonal of tridiagonal matrix
b ! right hand side of matrix solve
integer, intent(out) :: &
nyn ! matrix size
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqin0 , & ! ice layer enthalpy (J m-3) at beggining of timestep
Iswabs , & ! SW radiation absorbed in ice layers (W m-2)
phi , & ! ice layer liquid fraction
zqsn0 , & ! snow layer enthalpy (J m-3) at start of timestep
Sswabs ! SW radiation absorbed in snow layers (W m-2)
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! snow surface temperature (C)
dt , & ! timestep (s)
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
Tbot , & ! ice bottom surfce temperature (deg C)
qpond , & ! melt pond brine enthalpy (J m-3)
qocn , & ! ocean brine enthalpy (J m-3)
w , & ! downwards vertical flushing Darcy velocity (m/s)
fsurfn , & ! net flux to top surface, excluding fcondtop
dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(in) :: &
dxp , & ! distances between grid points (m)
kcstar ! interface conductivities (W m-1 K-1)
integer :: &
k , & ! vertical layer index
l , & ! matrix index
m ! vertical index
character(len=*),parameter :: subname='(matrix_elements_snow_cold)'
! surface temperature
l = 1
Ap(l) = dfsurfn_dTsf - kcstar(1) / dxp(1)
As(l) = kcstar(1) / dxp(1)
An(l) = c0
b (l) = dfsurfn_dTsf * Tsf - fsurfn
! surface layer
k = 1
l = k + 1
m = 1
Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
kcstar(m+1) / dxp(m+1) + &
kcstar(m) / dxp(m)
As(l) = -kcstar(m+1) / dxp(m+1)
An(l) = -kcstar(m) / dxp(m)
b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
! interior snow layers
do k = 2, nslyr
l = k + 1
m = k
Ap(l) = ((rhos * cp_ice) / dt) * hslyr + &
kcstar(m+1) / dxp(m+1) + &
kcstar(m) / dxp(m)
As(l) = -kcstar(m+1) / dxp(m+1)
An(l) = -kcstar(m) / dxp(m)
b (l) = ((rhos * Lfresh + zqsn0(k)) / dt) * hslyr + Sswabs(k)
enddo ! k
! top ice layer
k = 1
l = nslyr + k + 1
m = k + nslyr
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(m+1) / dxp(m+1) + &
kcstar(m) / dxp(m) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(m+1) / dxp(m+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(m) / dxp(m)
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
w * qpond
! interior ice layers
do k = 2, nilyr-1
l = nslyr + k + 1
m = k + nslyr
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(m+1) / dxp(m+1) + &
kcstar(m) / dxp(m) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = -kcstar(m+1) / dxp(m+1) - &
q(k) * cp_ocn * rhow
An(l) = -kcstar(m) / dxp(m) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k)
enddo ! k
! bottom layer
k = nilyr
l = nilyr + nslyr + 1
m = k + nslyr
Ap(l) = ((phi(k) * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice) / dt) * hilyr + &
kcstar(m+1) / dxp(m+1) + &
kcstar(m) / dxp(m) + &
q(k) * cp_ocn * rhow + &
w * cp_ocn * rhow
As(l) = c0
An(l) = -kcstar(m) / dxp(m) - &
w * cp_ocn * rhow
b (l) = (((c1 - phi(k)) * rhoi * Lfresh + zqin0(k)) / dt) * hilyr + Iswabs(k) + &
(kcstar(m+1) * Tbot) / dxp(m+1) + &
q(k) * qocn
nyn = nilyr + nslyr + 1
end subroutine matrix_elements_snow_cold
!=======================================================================
subroutine solve_salinity(zSin, Sbr, &
Spond, sss, &
q, dSdt, &
w, hilyr, &
dt, nilyr)
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zSin ! ice layer bulk salinity (ppt)
real(kind=dbl_kind), dimension(:), intent(in) :: &
Sbr , & ! ice layer brine salinity (ppt)
dSdt ! gravity drainage desalination rate for slow mode (ppt s-1)
real(kind=dbl_kind), dimension(0:nilyr), intent(in) :: &
q ! upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), intent(in) :: &
Spond , & ! melt pond salinity (ppt)
sss , & ! sea surface salinity (ppt)
w , & ! vertical flushing Darcy velocity (m/s)
hilyr , & ! ice layer thickness (m)
dt ! timestep (s)
integer :: &
k ! vertical layer index
real(kind=dbl_kind), parameter :: &
S_min = p01
real(kind=dbl_kind), dimension(nilyr) :: &
zSin0
character(len=*),parameter :: subname='(solve_salinity)'
zSin0 = zSin
k = 1
zSin(k) = zSin(k) + max(S_min - zSin(k), &
((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
dSdt(k) + &
(w * (Spond - Sbr(k))) / hilyr) * dt)
do k = 2, nilyr-1
zSin(k) = zSin(k) + max(S_min - zSin(k), &
((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + &
dSdt(k) + &
(w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
enddo ! k
k = nilyr
zSin(k) = zSin(k) + max(S_min - zSin(k), &
((q(k) * (sss - Sbr(k))) / hilyr + &
dSdt(k) + &
(w * (Sbr(k-1) - Sbr(k))) / hilyr) * dt)
if (minval(zSin) < c0) then
write(warnstr,*) subname, (q(k) * (Sbr(k+1) - Sbr(k))) / hilyr, &
dSdt(k) , &
(w * (Spond - Sbr(k))) / hilyr
call icepack_warnings_add(warnstr)
do k = 1, nilyr
write(warnstr,*) subname, k, zSin(k), zSin0(k)
call icepack_warnings_add(warnstr)
enddo
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
end subroutine solve_salinity
!=======================================================================
subroutine tdma_solve_sparse(nilyr, nslyr, a, b, c, d, x, n)
! perform a tri-diagonal solve with TDMA using a sparse tridiagoinal matrix
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nslyr ! number of snow layers
integer(kind=int_kind), intent(in) :: &
n ! matrix size
real(kind=dbl_kind), dimension(:), intent(in) :: &
a , & ! matrix lower off-diagonal
b , & ! matrix diagonal
c , & ! matrix upper off-diagonal
d ! right hand side vector
real(kind=dbl_kind), dimension(:), intent(out) :: &
x ! solution vector
real(kind=dbl_kind), dimension(nilyr+nslyr+1) :: &
cp , & ! modified upper off-diagonal vector
dp ! modified right hand side vector
integer(kind=int_kind) :: &
i ! vector index
character(len=*),parameter :: subname='(tdma_solve_sparse)'
! forward sweep
cp(1) = c(1) / b(1)
do i = 2, n-1
cp(i) = c(i) / (b(i) - cp(i-1)*a(i))
enddo
dp(1) = d(1) / b(1)
do i = 2, n
dp(i) = (d(i) - dp(i-1)*a(i)) / (b(i) - cp(i-1)*a(i))
enddo
! back substitution
x(n) = dp(n)
do i = n-1,1,-1
x(i) = dp(i) - cp(i)*x(i+1)
enddo
end subroutine tdma_solve_sparse
!=======================================================================
! Effect of salinity
!=======================================================================
function permeability(phi) result(perm)
! given the liquid fraction calculate the permeability
! See Golden et al. 2007
real(kind=dbl_kind), intent(in) :: &
phi ! liquid fraction
real(kind=dbl_kind) :: &
perm ! permeability (m2)
real(kind=dbl_kind), parameter :: &
phic = p05 ! critical liquid fraction for impermeability
character(len=*),parameter :: subname='(permeability)'
perm = 3.0e-8_dbl_kind * max(phi - phic, c0)**3
end function permeability
!=======================================================================
subroutine explicit_flow_velocities(nilyr, zSin, &
zTin, Tsf, &
Tbot, q, &
dSdt, Sbr, &
qbr, dt, &
sss, qocn, &
hilyr, hin)
! calculate the rapid gravity drainage mode Darcy velocity and the
! slow mode drainage rate
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zSin, & ! ice layer bulk salinity (ppt)
zTin ! ice layer temperature (C)
real(kind=dbl_kind), intent(in) :: &
Tsf , & ! ice/snow surface temperature (C)
Tbot , & ! ice bottom temperature (C)
dt , & ! time step (s)
sss , & ! sea surface salinty (ppt)
qocn , & ! ocean enthalpy (J m-3)
hilyr , & ! ice layer thickness (m)
hin ! ice thickness (m)
real(kind=dbl_kind), dimension(0:nilyr), intent(out) :: &
q ! rapid mode upward interface vertical Darcy flow (m s-1)
real(kind=dbl_kind), dimension(:), intent(out) :: &
dSdt ! slow mode drainage rate (ppt s-1)
real(kind=dbl_kind), dimension(:), intent(out) :: &
Sbr , & ! ice layer brine salinity (ppt)
qbr ! ice layer brine enthalpy (J m-3)
real(kind=dbl_kind), parameter :: &
kappal = 8.824e-8_dbl_kind, & ! heat diffusivity of liquid
fracmax = p2 , & ! limiting advective layer fraction
zSin_min = p1 , & ! minimum bulk salinity (ppt)
safety_factor = c10 ! to prevent negative salinities
real(kind=dbl_kind), dimension(1:nilyr) :: &
phi ! ice layer liquid fraction
real(kind=dbl_kind), dimension(0:nilyr) :: &
rho ! ice layer brine density (kg m-3)
real(kind=dbl_kind) :: &
rho_ocn , & ! ocean density (kg m-3)
perm_min , & ! minimum permeability from layer to ocean (m2)
perm_harm , & ! harmonic mean of permeability from layer to ocean (m2)
rho_sum , & ! sum of the brine densities from layer to ocean (kg m-3)
rho_pipe , & ! density of the brine in the channel (kg m-3)
z , & ! distance to layer from top surface (m)
perm , & ! ice layer permeability (m2)
drho , & ! brine density difference between layer and ocean (kg m-3)
Ra , & ! local mush Rayleigh number
rn , & ! real value of number of layers considered
L , & ! thickness of the layers considered (m)
dx , & ! horizontal size of convective flow (m)
dx2 , & ! square of the horizontal size of convective flow (m2)
Am , & ! A parameter for mush
Bm , & ! B parameter for mush
Ap , & ! A parameter for channel
Bp , & ! B parameter for channel
qlimit , & ! limit to vertical Darcy flow for numerical stability
dS_guess , & ! expected bulk salinity without limits
alpha , & ! desalination limiting factor
ra_constants ! for Rayleigh number
integer(kind=int_kind) :: &
k ! ice layer index
character(len=*),parameter :: subname='(explicit_flow_velocities)'
! initial downward sweep - determine derived physical quantities
do k = 1, nilyr
Sbr(k) = liquidus_brine_salinity_mush(zTin(k))
phi(k) = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
qbr(k) = enthalpy_brine(zTin(k))
rho(k) = icepack_mushy_density_brine(Sbr(k))
enddo ! k
rho(0) = rho(1)
! ocean conditions
Sbr(nilyr+1) = sss
qbr(nilyr+1) = qocn
rho_ocn = icepack_mushy_density_brine(sss)
! initialize accumulated quantities
perm_min = bignum
perm_harm = c0
rho_sum = c0
! limit to q for numerical stability
qlimit = (fracmax * hilyr) / dt
! no flow through ice top surface
q(0) = c0
ra_constants = gravit / (viscosity_dyn * kappal)
! first iterate over layers going up
do k = nilyr, 1, -1
! vertical position from ice top surface
z = ((real(k, dbl_kind) - p5) / real(nilyr, dbl_kind)) * hin
! permeabilities
perm = permeability(phi(k))
perm_min = min(perm_min,perm)
perm_harm = perm_harm + (c1 / max(perm,1.0e-30_dbl_kind))
! densities
rho_sum = rho_sum + rho(k)
!rho_pipe = rho(k)
rho_pipe = p5 * (rho(k) + rho(k-1))
drho = max(rho(k) - rho_ocn, c0)
! mush Rayleigh number
Ra = drho * (hin-z) * perm_min * ra_constants
! height of mush layer to layer k
rn = real(nilyr-k+1,dbl_kind)
L = rn * hilyr
! horizontal size of convection
dx = L * c2 * aspect_rapid_mode
dx2 = dx**2
! determine vertical Darcy flow
Am = (dx2 * rn) / (viscosity_dyn * perm_harm)
Bm = (-gravit * rho_sum) / rn
Ap = (pi * a_rapid_mode**4) / (c8 * viscosity_dyn)
Bp = -rho_pipe * gravit
q(k) = max((Am / dx2) * ((-Ap*Bp - Am*Bm) / (Am + Ap) + Bm), 1.0e-30_dbl_kind)
! modify by Rayleigh number and advection limit
q(k) = min(q(k) * (max(Ra - Rac_rapid_mode, c0) / (Ra+puny)), qlimit)
! late stage drainage
dSdt(k) = dSdt_slow_mode * (max((zSin(k) - phi_c_slow_mode*Sbr(k)), c0) &
* max((Tbot - Tsf), c0)) / (hin + 0.001_dbl_kind)
dSdt(k) = max(dSdt(k), (-zSin(k) * 0.5_dbl_kind) / dt)
! restrict flows to prevent too much salt loss
dS_guess = (((q(k) * (Sbr(k+1) - Sbr(k))) / hilyr + dSdt(k)) * dt) * safety_factor
if (abs(dS_guess) < puny) then
alpha = c1
else
alpha = (zSin_min - zSin(k)) / dS_guess
endif
if (alpha < c0 .or. alpha > c1) alpha = c1
q(k) = q(k) * alpha
dSdt(k) = dSdt(k) * alpha
enddo ! k
end subroutine explicit_flow_velocities
!=======================================================================
! Flushing
!=======================================================================
subroutine flushing_velocity(zTin, &
phi, nilyr, &
hin, hsn, &
hilyr, &
hpond, apond, &
dt, w)
! calculate the vertical flushing Darcy velocity (positive downward)
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real(kind=dbl_kind), dimension(:), intent(in) :: &
zTin , & ! ice layer temperature (C)
phi ! ice layer liquid fraction
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hpond , & ! melt pond thickness (m)
apond , & ! melt pond area (-)
hsn , & ! snow thickness (m)
hin , & ! ice thickness (m)
dt ! time step (s)
real(kind=dbl_kind), intent(out) :: &
w ! vertical flushing Darcy flow rate (m s-1)
real(kind=dbl_kind), parameter :: &
advection_limit = 0.005_dbl_kind ! limit to fraction of brine in
! any layer that can be advected
real(kind=dbl_kind) :: &
perm , & ! ice layer permeability (m2)
ice_mass , & ! mass of ice (kg m-2)
perm_harm , & ! harmonic mean of ice permeability (m2)
hocn , & ! ocean surface height above ice base (m)
hbrine , & ! brine surface height above ice base (m)
w_down_max , & ! maximum downward flushing Darcy flow rate (m s-1)
phi_min , & ! minimum porosity in the mush
wlimit , & ! limit to w to avoid advecting all brine in layer
dhhead ! hydraulic head (m)
integer(kind=int_kind) :: &
k ! ice layer index
character(len=*),parameter :: subname='(flushing_velocity)'
! initialize
w = c0
! only flush if ponds are active
if (tr_pond) then
ice_mass = c0
perm_harm = c0
phi_min = c1
do k = 1, nilyr
! liquid fraction
!phi = icepack_mushy_liquid_fraction(zTin(k), zSin(k))
phi_min = min(phi_min,phi(k))
! permeability
perm = permeability(phi(k))
! ice mass
ice_mass = ice_mass + phi(k) * icepack_mushy_density_brine(liquidus_brine_salinity_mush(zTin(k))) + &
(c1 - phi(k)) * rhoi
! permeability harmonic mean
perm_harm = perm_harm + c1 / (perm + 1e-30_dbl_kind)
enddo ! k
ice_mass = ice_mass * hilyr
perm_harm = real(nilyr,dbl_kind) / perm_harm
! calculate ocean surface height above bottom of ice
hocn = (ice_mass + hpond * apond * rhow + hsn * rhos) / rhow
! calculate brine height above bottom of ice
hbrine = hin + hpond
! pressure head
dhhead = max(hbrine - hocn,c0)
! darcy flow through ice
w = (perm_harm * rhow * gravit * (dhhead / hin)) / viscosity_dyn
! maximum down flow to drain pond
w_down_max = (hpond * apond) / dt
! limit flow
w = min(w,w_down_max)
! limit amount of brine that can be advected out of any particular layer
wlimit = (advection_limit * phi_min * hilyr) / dt
if (abs(w) > puny) then
w = w * max(min(abs(wlimit/w),c1),c0)
else
w = c0
endif
w = max(w, c0)
endif
end subroutine flushing_velocity
!=======================================================================
subroutine flush_pond(w, hpond, apond, dt)
! given a flushing velocity drain the meltponds
real(kind=dbl_kind), intent(in) :: &
w , & ! vertical flushing Darcy flow rate (m s-1)
apond , & ! melt pond area (-)
dt ! time step (s)
real(kind=dbl_kind), intent(inout) :: &
hpond ! melt pond thickness (m)
real(kind=dbl_kind), parameter :: &
lambda_pond = c1 / (10.0_dbl_kind * 24.0_dbl_kind * 3600.0_dbl_kind), &
hpond0 = 0.01_dbl_kind
character(len=*),parameter :: subname='(flush_pond)'
if (tr_pond) then
if (apond > c0 .and. hpond > c0) then
! flush pond through mush
hpond = hpond - w * dt / apond
hpond = max(hpond, c0)
! exponential decay of pond
hpond = hpond - lambda_pond * dt * (hpond + hpond0)
hpond = max(hpond, c0)
endif
endif
end subroutine flush_pond
!=======================================================================
subroutine flood_ice(hsn, hin, &
nslyr, nilyr, &
hslyr, hilyr, &
zqsn, zqin, &
phi, dt, &
zSin, Sbr, &
sss, qocn, &
snoice, fadvheat)
! given upwards flushing brine flow calculate amount of snow ice and
! convert snow to ice with appropriate properties
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real(kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
hsn , & ! snow thickness (m)
hin , & ! ice thickness (m)
sss , & ! sea surface salinity (ppt)
qocn ! ocean brine enthalpy (J m-2)
real(kind=dbl_kind), dimension(:), intent(inout) :: &
zqsn , & ! snow layer enthalpy (J m-2)
zqin , & ! ice layer enthalpy (J m-2)
zSin , & ! ice layer bulk salinity (ppt)
phi ! ice liquid fraction
real(kind=dbl_kind), dimension(:), intent(in) :: &
Sbr ! ice layer brine salinity (ppt)
real(kind=dbl_kind), intent(inout) :: &
hslyr , & ! snow layer thickness (m)
hilyr ! ice layer thickness (m)
real(kind=dbl_kind), intent(out) :: &
snoice ! snow ice formation
real(kind=dbl_kind), intent(inout) :: &
fadvheat ! advection heat flux to ocean
real(kind=dbl_kind) :: &
hin2 , & ! new ice thickness (m)
hsn2 , & ! new snow thickness (m)
hilyr2 , & ! new ice layer thickness (m)
hslyr2 , & ! new snow layer thickness (m)
dh , & ! thickness of snowice formation (m)
phi_snowice , & ! liquid fraction of new snow ice
rho_snowice , & ! density of snowice (kg m-3)
zSin_snowice , & ! bulk salinity of new snowice (ppt)
zqin_snowice , & ! ice enthalpy of new snowice (J m-2)
zqsn_snowice , & ! snow enthalpy of snow thats becoming snowice (J m-2)
freeboard_density , & ! negative of ice surface freeboard times the ocean density (kg m-2)
ice_mass , & ! mass of the ice (kg m-2)
rho_ocn , & ! density of the ocean (kg m-3)
ice_density , & ! density of ice layer (kg m-3)
hadded , & ! thickness rate of water used from ocean (m/s)
wadded , & ! mass rate of water used from ocean (kg/m^2/s)
eadded ! energy rate of water used from ocean (W/m^2)
! real(kind=dbl_kind) :: &
! sadded ! salt rate of water used from ocean (kg/m^2/s)
integer :: &
k ! vertical index
character(len=*),parameter :: subname='(flood_ice)'
snoice = c0
! check we have snow
if (hsn > puny) then
rho_ocn = icepack_mushy_density_brine(sss)
! ice mass
ice_mass = c0
do k = 1, nilyr
ice_density = min(phi(k) * icepack_mushy_density_brine(Sbr(k)) + (c1 - phi(k)) * rhoi,rho_ocn)
ice_mass = ice_mass + ice_density
enddo ! k
ice_mass = ice_mass * hilyr
! negative freeboard times ocean density
freeboard_density = max(ice_mass + hsn * rhos - hin * rho_ocn, c0)
! check if have flooded ice
if (freeboard_density > c0) then
! sea ice fraction of newly formed snow ice
phi_snowice = (c1 - rhos / rhoi)
! density of newly formed snowice
rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
! calculate thickness of new ice added
dh = freeboard_density / (rho_ocn - rho_snowice + rhos)
dh = max(min(dh,hsn),c0)
! enthalpy of snow that becomes snowice
call enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
if (icepack_warnings_aborted(subname)) return
! change thicknesses
hin2 = hin + dh
hsn2 = hsn - dh
hilyr2 = hin2 / real(nilyr,dbl_kind)
hslyr2 = hsn2 / real(nslyr,dbl_kind)
! properties of new snow ice
zSin_snowice = phi_snowice * sss
zqin_snowice = phi_snowice * qocn + zqsn_snowice
! change snow properties
call update_vertical_tracers_snow(nslyr, zqsn, hslyr, hslyr2)
if (icepack_warnings_aborted(subname)) return
! change ice properties
call update_vertical_tracers_ice(nilyr, zqin, hilyr, hilyr2, &
hin, hin2, zqin_snowice)
if (icepack_warnings_aborted(subname)) return
call update_vertical_tracers_ice(nilyr, zSin, hilyr, hilyr2, &
hin, hin2, zSin_snowice)
if (icepack_warnings_aborted(subname)) return
call update_vertical_tracers_ice(nilyr, phi, hilyr, hilyr2, &
hin, hin2, phi_snowice)
if (icepack_warnings_aborted(subname)) return
! change thicknesses
hilyr = hilyr2
hslyr = hslyr2
snoice = dh
hadded = (dh * phi_snowice) / dt
wadded = hadded * rhoi
eadded = hadded * qocn
! sadded = wadded * ice_ref_salinity * p001
! conservation
fadvheat = fadvheat - eadded
endif
endif
end subroutine flood_ice
!=======================================================================
subroutine enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
! determine enthalpy of the snow being converted to snow ice
integer (kind=int_kind), intent(in) :: &
nslyr ! number of snow layers
real(kind=dbl_kind), intent(in) :: &
dh , & ! thickness of new snowice formation (m)
hsn ! initial snow thickness
real(kind=dbl_kind), dimension(:), intent(in) :: &
zqsn ! snow layer enthalpy (J m-2)
real(kind=dbl_kind), intent(out) :: &
zqsn_snowice ! enthalpy of snow becoming snowice (J m-2)
real(kind=dbl_kind) :: &
rnlyr ! real value of number of snow layers turning to snowice
integer(kind=int_kind) :: &
nlyr , & ! no of snow layers completely converted to snowice
k ! snow layer index
character(len=*),parameter :: subname='(enthalpy_snow_snowice)'
zqsn_snowice = c0
! snow depth and snow layers affected by snowice formation
if (hsn > puny) then
rnlyr = (dh / hsn) * nslyr
nlyr = min(floor(rnlyr),nslyr-1) ! nlyr=0 if nslyr=1
! loop over full snow layers affected
! not executed if nlyr=0
do k = nslyr, nslyr-nlyr+1, -1
zqsn_snowice = zqsn_snowice + zqsn(k) / rnlyr
enddo ! k
! partially converted snow layer
zqsn_snowice = zqsn_snowice + &
((rnlyr - real(nlyr,dbl_kind)) / rnlyr) * zqsn(nslyr-nlyr)
endif
end subroutine enthalpy_snow_snowice
!=======================================================================
subroutine update_vertical_tracers_snow(nslyr, trc, hlyr1, hlyr2)
! given some snow ice formation regrid snow layers
integer (kind=int_kind), intent(in) :: &
nslyr ! number of snow layers
real(kind=dbl_kind), dimension(:), intent(inout) :: &
trc ! vertical tracer
real(kind=dbl_kind), intent(in) :: &
hlyr1 , & ! old cell thickness
hlyr2 ! new cell thickness
real(kind=dbl_kind), dimension(1:nslyr) :: &
trc2 ! temporary array for updated tracer
! vertical indexes for old and new grid
integer(kind=int_kind) :: &
k1 , & ! vertical index for old grid
k2 ! vertical index for new grid
real(kind=dbl_kind) :: &
z1a , & ! lower boundary of old cell
z1b , & ! upper boundary of old cell
z2a , & ! lower boundary of new cell
z2b , & ! upper boundary of new cell
overlap ! overlap between old and new cell
character(len=*),parameter :: subname='(update_vertical_tracers_snow)'
! loop over new grid cells
do k2 = 1, nslyr
! initialize new tracer
trc2(k2) = c0
! calculate upper and lower boundary of new cell
z2a = (k2 - 1) * hlyr2
z2b = k2 * hlyr2
! loop over old grid cells
do k1 = 1, nslyr
! calculate upper and lower boundary of old cell
z1a = (k1 - 1) * hlyr1
z1b = k1 * hlyr1
! 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
! renormalize new grid cell
trc2(k2) = trc2(k2) / hlyr2
enddo ! k2
! update vertical tracer array with the adjusted tracer
trc = trc2
end subroutine update_vertical_tracers_snow
!=======================================================================
subroutine update_vertical_tracers_ice(nilyr, trc, hlyr1, hlyr2, &
h1, h2, trc0)
! given some snow ice formation regrid ice layers
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) :: &
hlyr1 , & ! old cell thickness
hlyr2 , & ! new cell thickness
h1 , & ! old total thickness
h2 , & ! new total thickness
trc0 ! tracer value of added snow ice on ice top
real(kind=dbl_kind), dimension(1:nilyr) :: &
trc2 ! temporary array for updated tracer
integer(kind=int_kind) :: &
k1 , & ! vertical indexes for old grid
k2 ! vertical indexes for new grid
real(kind=dbl_kind) :: &
z1a , & ! lower boundary of old cell
z1b , & ! upper boundary of old cell
z2a , & ! lower boundary of new cell
z2b , & ! upper boundary of new cell
overlap ! overlap between old and new cell
character(len=*),parameter :: subname='(update_vertical_tracers_ice)'
! 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) * hlyr2
z2b = k2 * hlyr2
! calculate upper and lower boundary of added snow ice at top
z1a = c0
z1b = h2 - h1
! 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
! loop over old grid cells
do k1 = 1, nilyr
! calculate upper and lower boundary of old cell
z1a = (k1 - 1) * hlyr1 + h2 - h1
z1b = k1 * hlyr1 + h2 - h1
! 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
! renormalize new grid cell
trc2(k2) = trc2(k2) / hlyr2
enddo ! k2
! update vertical tracer array with the adjusted tracer
trc = trc2
end subroutine update_vertical_tracers_ice
!=======================================================================
end module icepack_therm_mushy
!=======================================================================