Skip to content

Commit

Permalink
fixed hocn and flushing velocity issue, passes basic restart test
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidclemenssewall
davidclemenssewall committed Nov 8, 2024
1 parent 6da5668 commit 33aa692
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 28 deletions.
54 changes: 35 additions & 19 deletions columnphysics/icepack_therm_mushy.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,12 @@ 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: viscosity_dyn, rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh, gravit, rhofresh
use icepack_parameters, only: hs_min, snwgrain
use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode, tscale_pnd_drain
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_tracers, only: nilyr, nslyr, tr_pond
use icepack_tracers, only: nilyr, nslyr, tr_pond, tr_pond_lvl
use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, icepack_enthalpy_snow
use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction
use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
Expand Down Expand Up @@ -46,7 +46,7 @@ subroutine temperature_changes_salinity(dt, &
fswsfc, fswint, &
Sswabs, Iswabs, &
hilyr, hslyr, &
apond, hpond, &
apnd, hpond, &
zqin, zTin, &
zqsn, zTsn, &
zSin, &
Expand All @@ -56,7 +56,8 @@ subroutine temperature_changes_salinity(dt, &
flwoutn, fsurfn, &
fcondtop, fcondbot, &
fadvheat, snoice, &
smice, smliq)
smice, smliq, &
alvl)

! solve the enthalpy and bulk salinity of the ice for a single column

Expand All @@ -71,7 +72,8 @@ subroutine temperature_changes_salinity(dt, &
shcoef , & ! transfer coefficient for sensible heat
lhcoef , & ! transfer coefficient for latent heat
Tbot , & ! ice bottom surfce temperature (deg C)
sss ! sea surface salinity (PSU)
sss , & ! sea surface salinity (PSU)
alvl ! melt pond area fraction

real (kind=dbl_kind), intent(inout) :: &
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
Expand All @@ -80,7 +82,7 @@ subroutine temperature_changes_salinity(dt, &
real (kind=dbl_kind), intent(inout) :: &
hilyr , & ! ice layer thickness (m)
hslyr , & ! snow layer thickness (m)
apond , & ! melt pond area fraction
apnd , & ! melt pond area fraction tracer
hpond ! melt pond depth (m)

real (kind=dbl_kind), dimension (:), intent(inout) :: &
Expand Down Expand Up @@ -182,8 +184,8 @@ subroutine temperature_changes_salinity(dt, &
! calculate vertical bulk darcy flow
call flushing_velocity(zTin, phi, &
hin, hsn, &
hilyr, &
hpond, apond, &
hilyr, alvl, &
hpond, apnd, &
dt, w)
if (icepack_warnings_aborted(subname)) return

Expand Down Expand Up @@ -328,7 +330,7 @@ subroutine temperature_changes_salinity(dt, &
endif

! drain ponds from flushing
call flush_pond(w, hpond, apond, dt)
call flush_pond(w, hpond, apnd, dt, alvl)
if (icepack_warnings_aborted(subname)) return

! flood snow ice
Expand Down Expand Up @@ -3064,8 +3066,8 @@ end subroutine explicit_flow_velocities

subroutine flushing_velocity(zTin, phi, &
hin, hsn, &
hilyr, &
hpond, apond, &
hilyr, alvl, &
hpond, apnd, &
dt, w)

! calculate the vertical flushing Darcy velocity (positive downward)
Expand All @@ -3076,8 +3078,9 @@ subroutine flushing_velocity(zTin, phi, &

real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
alvl , & ! level ice area tracer
hpond , & ! melt pond thickness (m)
apond , & ! melt pond area (-)
apnd , & ! melt pond area (-)
hsn , & ! snow thickness (m)
hin , & ! ice thickness (m)
dt ! time step (s)
Expand Down Expand Up @@ -3138,7 +3141,11 @@ subroutine flushing_velocity(zTin, phi, &
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
if (tr_pond_lvl) then
hocn = (ice_mass + hpond * apnd * rhofresh * alvl + hsn * rhos) / rhow
else
hocn = (ice_mass + hpond * apnd * rhofresh + hsn * rhos) / rhow
endif

! calculate brine height above bottom of ice
hbrine = hin + hpond
Expand All @@ -3150,7 +3157,11 @@ subroutine flushing_velocity(zTin, phi, &
w = (perm_harm * rhow * gravit * (dhhead / hin)) / viscosity_dyn

! maximum down flow to drain pond
w_down_max = (hpond * apond) / dt
if (tr_pond_lvl) then
w_down_max = (hpond * apnd * alvl) / dt
else
w_down_max = (hpond * apnd) / dt
endif

! limit flow
w = min(w,w_down_max)
Expand All @@ -3172,14 +3183,15 @@ end subroutine flushing_velocity

!=======================================================================

subroutine flush_pond(w, hpond, apond, dt)
subroutine flush_pond(w, hpond, apnd, dt, alvl)

! 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)
apnd , & ! melt pond area (-)
dt , & ! time step (s)
alvl ! level ice area tracer

real(kind=dbl_kind), intent(inout) :: &
hpond ! melt pond thickness (m)
Expand All @@ -3193,10 +3205,14 @@ subroutine flush_pond(w, hpond, apond, dt)
character(len=*),parameter :: subname='(flush_pond)'

if (tr_pond) then
if (apond > c0 .and. hpond > c0) then
if (apnd > c0 .and. hpond > c0) then

! flush pond through mush
hpond = hpond - w * dt / apond
if (tr_pond_lvl) then
hpond = hpond - w * dt / (apnd * alvl)
else
hpond = hpond - w * dt / apnd
endif

hpond = max(hpond, c0)

Expand Down
20 changes: 11 additions & 9 deletions columnphysics/icepack_therm_vertical.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ subroutine thermo_vertical (dt, aicen, &
vicen, vsnon, &
Tsf, zSin, &
zqin, zqsn, &
apond, hpond, &
apnd, hpond, &
flw, potT, &
Qa, rhoa, &
fsnow, fpond, &
Expand All @@ -104,7 +104,7 @@ subroutine thermo_vertical (dt, aicen, &
congel, snoice, &
mlt_onset, frz_onset, &
yday, dsnow, &
prescribed_ice)
prescribed_ice, alvl)

real (kind=dbl_kind), intent(in) :: &
dt , & ! time step
Expand All @@ -119,8 +119,9 @@ subroutine thermo_vertical (dt, aicen, &
! tracers
real (kind=dbl_kind), intent(inout) :: &
Tsf , & ! ice/snow top surface temp, same as Tsfcn (deg C)
apond , & ! melt pond area fraction
hpond ! melt pond depth (m)
apnd , & ! melt pond area fraction
hpond , & ! melt pond depth (m)
alvl ! level ice area fraction
! iage ! ice age (s)

logical (kind=log_kind), intent(in), optional :: &
Expand Down Expand Up @@ -311,7 +312,7 @@ subroutine thermo_vertical (dt, aicen, &
fswsfc, fswint, &
Sswabs, Iswabs, &
hilyr, hslyr, &
apond, hpond, &
apnd, hpond, &
zqin, zTin, &
zqsn, zTsn, &
zSin, &
Expand All @@ -321,7 +322,8 @@ subroutine thermo_vertical (dt, aicen, &
flwoutn, fsurfn, &
fcondtopn, fcondbotn, &
fadvocn, snoice, &
smice, smliq)
smice, smliq, &
alvl)
if (icepack_warnings_aborted(subname)) return

else ! ktherm
Expand Down Expand Up @@ -448,7 +450,7 @@ subroutine thermo_vertical (dt, aicen, &
fhocnn = fhocnn + fadvocn ! for ktherm=2

if (hin == c0) then
if (tr_pond_topo) fpond = fpond - aicen * apond * hpond
if (tr_pond_topo) fpond = fpond - aicen * apnd * hpond
endif

!-----------------------------------------------------------------
Expand Down Expand Up @@ -2620,7 +2622,7 @@ subroutine icepack_step_therm1(dt, &
vicen=vicen (n), vsnon=vsnon (n), &
Tsf=Tsfc (n), zSin=zSin (:,n), &
zqin=zqin (:,n), zqsn=zqsn (:,n), &
apond=apnd (n), hpond=hpnd (n), &
apnd=apnd (n), hpond=hpnd (n), &
flw=flw, potT=potT, &
Qa=Qa, rhoa=rhoa, &
fsnow=fsnow, fpond=fpond, &
Expand All @@ -2644,7 +2646,7 @@ subroutine icepack_step_therm1(dt, &
congel=congeln (n), snoice=snoicen (n), &
mlt_onset=mlt_onset, frz_onset=frz_onset, &
yday=yday, dsnow=dsnown (n), &
prescribed_ice=prescribed_ice)
prescribed_ice=prescribed_ice, alvl=alvl (n))

if (icepack_warnings_aborted(subname)) then
write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n
Expand Down

0 comments on commit 33aa692

Please sign in to comment.