Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

hocn calculation and flushing velocity issue #504

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open

hocn calculation and flushing velocity issue #504

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
33aa692
fixed hocn and flushing velocity issue, passes basic restart test
davidclemenssewall Nov 8, 2024
0b09d92
Merge branch 'main' into hocn_mushy_thermo_fix
davidclemenssewall Nov 11, 2024
5f53385
fixed indentation and relocated alvl in calling arguments
davidclemenssewall Nov 12, 2024
29f1d2d
Merge branch 'main' into hocn_mushy_thermo_fix
davidclemenssewall Dec 6, 2024
f0e1cc2
moved apond = apnd*alvl calculation to icepack_step_therm1
davidclemenssewall Dec 6, 2024
2392676
changed apnd for apond for neutral_drag_coeffs
davidclemenssewall Dec 6, 2024
f35a56d
removed extraneous tr_pond_lvl
davidclemenssewall Dec 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 35 additions & 19 deletions columnphysics/icepack_therm_mushy.F90
View file
Open in desktop
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
eclare108213 marked this conversation as resolved.
Show resolved Hide resolved
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix indentation

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ subroutine thermo_vertical (dt, aicen, &
vicen, vsnon, &
Tsf, zSin, &
zqin, zqsn, &
apond, hpond, &
apnd, hpond, &
flw, potT, &
Qa, rhoa, &
fsnow, fpond, &
Expand All @@ -107,7 +107,7 @@ subroutine thermo_vertical (dt, aicen, &
congel, snoice, &
mlt_onset, frz_onset, &
yday, dsnow, &
prescribed_ice)
prescribed_ice, alvl)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggest moving alvl up in the calling arguments. Makes no difference technically, but seems like it belongs with apnd and hpond, for instance.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.


real (kind=dbl_kind), intent(in) :: &
dt , & ! time step
Expand All @@ -122,8 +122,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 @@ -314,7 +315,7 @@ subroutine thermo_vertical (dt, aicen, &
fswsfc, fswint, &
Sswabs, Iswabs, &
hilyr, hslyr, &
apond, hpond, &
apnd, hpond, &
zqin, zTin, &
zqsn, zTsn, &
zSin, &
Expand All @@ -324,7 +325,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 @@ -451,7 +453,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 @@ -2692,7 +2694,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 @@ -2716,7 +2718,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
Loading