Skip to content

Commit

Permalink
Corrects thin ice/snow treatment of enthalpy and other tracers (CICE-…
Browse files Browse the repository at this point in the history
…Consortium#454)

* Corrects thin ice/snow treatment of enthalpy and other tracers

Adopted from E3SM-Project/E3SM#5630

"This fix redistributes enthalpy and other tracers evenly in snow and ice when their
respective thicknesses are < 1e-15 . Previously, these tracers were zero'd non-conservatively.
Also corrects a bug in the zeroing of snow thicknesses, and removes snow in thickness_changes
if the ice vanishes."

* Clean up modifications
  • Loading branch information
apcraig authored Aug 30, 2023
1 parent f5e093f commit 23b6c12
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 29 deletions.
76 changes: 48 additions & 28 deletions columnphysics/icepack_therm_shared.F90
Original file line number Diff line number Diff line change
Expand Up @@ -497,7 +497,8 @@ subroutine adjust_enthalpy (nlyr, &
hovlp ! overlap between old and new layers (m)

real (kind=dbl_kind) :: &
rhlyr ! 1./hlyr
rhlyr, & ! 1./hlyr
qtot ! total h*q in the column

real (kind=dbl_kind), dimension (nlyr) :: &
hq ! h * q for a layer
Expand All @@ -509,36 +510,55 @@ subroutine adjust_enthalpy (nlyr, &
!-----------------------------------------------------------------

rhlyr = c0
if (hn > puny) rhlyr = c1 / hlyr

!-----------------------------------------------------------------
! Compute h*q for new layers (k2) given overlap with old layers (k1)
!-----------------------------------------------------------------

do k2 = 1, nlyr
hq(k2) = c0
enddo ! k
k1 = 1
k2 = 1
do while (k1 <= nlyr .and. k2 <= nlyr)
hovlp = min (z1(k1+1), z2(k2+1)) &
- max (z1(k1), z2(k2))
hovlp = max (hovlp, c0)
hq(k2) = hq(k2) + hovlp*qn(k1)
if (z1(k1+1) > z2(k2+1)) then
k2 = k2 + 1
if (hn > puny) then
rhlyr = c1 / hlyr

!-----------------------------------------------------------------
! Compute h*q for new layers (k2) given overlap with old layers (k1)
!-----------------------------------------------------------------

do k2 = 1, nlyr
hq(k2) = c0
enddo ! k
k1 = 1
k2 = 1
do while (k1 <= nlyr .and. k2 <= nlyr)
hovlp = min (z1(k1+1), z2(k2+1)) &
- max (z1(k1), z2(k2))
hovlp = max (hovlp, c0)
hq(k2) = hq(k2) + hovlp*qn(k1)
if (z1(k1+1) > z2(k2+1)) then
k2 = k2 + 1
else
k1 = k1 + 1
endif
enddo ! while

!-----------------------------------------------------------------
! Compute new enthalpies.
!-----------------------------------------------------------------

do k = 1, nlyr
qn(k) = hq(k) * rhlyr
enddo ! k

else

qtot = c0
do k = 1, nlyr
qtot = qtot + qn(k) * (z1(k+1)-z1(k))
enddo
if (hn > c0) then
do k = 1, nlyr
qn(k) = qtot/hn
enddo
else
k1 = k1 + 1
do k = 1, nlyr
qn(k) = c0
enddo
endif
enddo ! while

!-----------------------------------------------------------------
! Compute new enthalpies.
!-----------------------------------------------------------------

do k = 1, nlyr
qn(k) = hq(k) * rhlyr
enddo ! k
endif

end subroutine adjust_enthalpy

Expand Down
5 changes: 4 additions & 1 deletion columnphysics/icepack_therm_vertical.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1706,7 +1706,7 @@ subroutine thickness_changes (nilyr, nslyr, &
!-----------------------------------------------------------------

if (ktherm == 2) then
if (hsn <= puny) then
if (hsn <= puny .or. hin <= c0) then
do k = 1, nslyr
fhocnn = fhocnn &
+ zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
Expand All @@ -1715,8 +1715,11 @@ subroutine thickness_changes (nilyr, nslyr, &
meltsliq = meltsliq + massice(k) ! add to meltponds
smice(k) = rhos
smliq(k) = c0
rsnw(k) = rsnw_fall
endif
enddo
melts = melts + hsn
hsn = c0
hslyr = c0
endif
endif
Expand Down

0 comments on commit 23b6c12

Please sign in to comment.