Skip to content

fix: correct physics bugs in 3-equation ice cavity melt model #875

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

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

fix: correct physics bugs in 3-equation ice cavity melt model #875

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2b65e36
fix: correct physics bugs in 3-equation ice cavity melt model
JanStreffing Apr 2, 2026
1be3dfa
Merge branch 'main' into fix/cavity-3eq-issues
JanStreffing Apr 2, 2026
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
12 changes: 6 additions & 6 deletions src/cavity_param.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,8 +225,8 @@ subroutine cavity_heat_water_fluxes_3eq(ice, dynamics, tracers, partit, mesh)
real(kind=WP),parameter :: tob= -20. !temperatur at the ice surface
real(kind=WP),parameter :: rhoi= 920. !mean ice density
real(kind=WP),parameter :: cpw = 4180.0 !Barnier et al. (1995)
real(kind=WP),parameter :: lhf = 3.33e+5 !latent heat of fusion
real(kind=WP),parameter :: tdif= 1.54e-6 !thermal conductivity of ice shelf !RG4190 / RG44027
real(kind=WP),parameter :: lhf = 334000. !latent heat of fusion [J/kg]
real(kind=WP),parameter :: tdif= 1.54e-6 !thermal diffusivity of ice shelf [m2/s] !RG4190 / RG44027
real(kind=WP),parameter :: atk = 273.15 !0 deg C in Kelvin
real(kind=WP),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers"

Expand Down Expand Up @@ -262,7 +262,7 @@ subroutine cavity_heat_water_fluxes_3eq(ice, dynamics, tracers, partit, mesh)
!_______________________________________________________________________
! Calculate the in-situ temperature tin
!call potit(s(i,j,N,lrhs)+35.0,t(i,j,N,lrhs),-zice(i,j),rp,tin)
call potit(sal,temp,abs(zice),rp,tin)
call potit(sal,temp,abs(zice)*1.025_WP,rp,tin) ! convert depth [m] to pressure [dbar]

!_______________________________________________________________________
! Calculate or prescribe the turbulent heat and salt transfer coeff. GAT and GAS
Expand All @@ -275,7 +275,7 @@ subroutine cavity_heat_water_fluxes_3eq(ice, dynamics, tracers, partit, mesh)
vt1 = sqrt(UVnode(1,nzmin,node)*UVnode(1,nzmin,node)+UVnode(2,nzmin,node)*UVnode(2,nzmin,node))
vt1 = max(vt1,0.001_WP)
!vt1 = max(vt1,0.005) ! CW
re = 10._WP/un !vt1*re (=velocity times length scale over kinematic viscosity) is the Reynolds number
re = hnode(nzmin,node)/un !vt1*re is the Reynolds number; use actual top-layer thickness [m]

gats1= sak1*vt1
gats2= 2.12_WP*log(gats1*re)-9._WP
Expand All @@ -301,7 +301,7 @@ subroutine cavity_heat_water_fluxes_3eq(ice, dynamics, tracers, partit, mesh)
ep1 = cpw*gat
ep2 = cpi*gas
ep3 = lhf*gas
ep31 = -rhor*cpi*tdif/zice !RG4190 / RG44027 ! CW
ep31 = -rhor*rhor*cpi*tdif/zice !RG4190 / RG44027 ! CW; rhor^2 because D_ice = D_draft*(rhow/rhoi) = D_draft/rhor
ep4 = b+c*zice
ep5 = gas/rhor

Expand Down Expand Up @@ -336,7 +336,7 @@ subroutine cavity_heat_water_fluxes_3eq(ice, dynamics, tracers, partit, mesh)
ex4 = ex2/ex1
ex5 = ex3/ex1

sr1 = 0.25_WP*ex4*ex4-ex5
sr1 = max(0.25_WP*ex4*ex4-ex5, 0._WP) ! guard against negative discriminant
!sr2 = -0.5*ex4
sr2 = ex6*ex4 ! modified for RG4190 / RG44027 ! CW
sf1 = sr2+sqrt(sr1)
Expand Down
Loading