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.

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

Merged
18 changes: 9 additions & 9 deletions columnphysics/icepack_therm_mushy.F90
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 @@ -72,15 +72,15 @@ subroutine temperature_changes_salinity(dt, &
lhcoef , & ! transfer coefficient for latent heat
Tbot , & ! ice bottom surfce temperature (deg C)
sss ! sea surface salinity (PSU)

davidclemenssewall marked this conversation as resolved.
Show resolved Hide resolved
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
apond , & ! melt pond area fraction of category
hpond ! melt pond depth (m)

real (kind=dbl_kind), dimension (:), intent(inout) :: &
Expand Down Expand Up @@ -3077,7 +3077,7 @@ subroutine flushing_velocity(zTin, phi, &
real(kind=dbl_kind), intent(in) :: &
hilyr , & ! ice layer thickness (m)
hpond , & ! melt pond thickness (m)
apond , & ! melt pond area (-)
apond , & ! melt pond area fraction of category (-)
hsn , & ! snow thickness (m)
hin , & ! ice thickness (m)
dt ! time step (s)
Expand Down Expand Up @@ -3138,8 +3138,8 @@ 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

hocn = (ice_mass + hpond * apond * rhofresh + hsn * rhos) / rhow
davidclemenssewall marked this conversation as resolved.
Show resolved Hide resolved
! calculate brine height above bottom of ice
hbrine = hin + hpond

Expand All @@ -3151,7 +3151,7 @@ subroutine flushing_velocity(zTin, phi, &

! maximum down flow to drain pond
w_down_max = (hpond * apond) / dt

davidclemenssewall marked this conversation as resolved.
Show resolved Hide resolved
! limit flow
w = min(w,w_down_max)

Expand All @@ -3178,7 +3178,7 @@ subroutine flush_pond(w, hpond, apond, dt)

real(kind=dbl_kind), intent(in) :: &
w , & ! vertical flushing Darcy flow rate (m s-1)
apond , & ! melt pond area (-)
apond , & ! melt pond area fraction of category (-)
dt ! time step (s)

real(kind=dbl_kind), intent(inout) :: &
Expand Down
26 changes: 23 additions & 3 deletions columnphysics/icepack_therm_vertical.F90
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ 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
apond , & ! melt pond area fraction of category
hpond ! melt pond depth (m)
! iage ! ice age (s)

Expand Down Expand Up @@ -2336,7 +2336,7 @@ subroutine icepack_step_therm1(dt, &
Tsfc , & ! ice/snow surface temperature, Tsfcn
alvl , & ! level ice area fraction
vlvl , & ! level ice volume fraction
apnd , & ! melt pond area fraction
apnd , & ! melt pond area fraction tracer
hpnd , & ! melt pond depth (m)
ipnd , & ! melt pond refrozen lid thickness (m)
iage , & ! volume-weighted ice age
Expand Down Expand Up @@ -2437,6 +2437,7 @@ subroutine icepack_step_therm1(dt, &
smliq ! tracer for mass of liquid in snow (kg/m^3)

real (kind=dbl_kind), dimension(ncat) :: &
apond , & ! melt pond area fraction of category
l_meltsliqn ! mass of snow melt local (kg/m^2)

real (kind=dbl_kind) :: &
Expand Down Expand Up @@ -2526,6 +2527,17 @@ subroutine icepack_step_therm1(dt, &
massicen(:,:) = c0
massliqn(:,:) = c0

!-----------------------------------------------------------------
! Initialize pond area fractions
!-----------------------------------------------------------------
do n= 1, ncat
if (tr_pond_lvl) then
apond(n) = apnd(n) * alvl(n)
else
apond(n) = apnd(n)
endif
enddo

!-----------------------------------------------------------------
! Initialize rate of snow loss to leads
!-----------------------------------------------------------------
Expand Down Expand Up @@ -2556,6 +2568,7 @@ subroutine icepack_step_therm1(dt, &
!-----------------------------------------------------------------

if (formdrag) then
!!!! This should take apond, not apnd, will fix in the next commit
call neutral_drag_coeffs (apnd , &
hpnd , ipnd , &
alvl , vlvl , &
Expand Down Expand Up @@ -2704,7 +2717,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), &
apond=apond (n), hpond=hpnd (n), &
flw=flw, potT=potT, &
Qa=Qa, rhoa=rhoa, &
fsnow=fsnow, fpond=fpond, &
Expand Down Expand Up @@ -2736,6 +2749,13 @@ subroutine icepack_step_therm1(dt, &
return
endif

! Translate changes in apond into apnd tracer
if (tr_pond_lvl) then
apnd(n) = apond(n) / alvl(n)
else
apnd(n) = apond(n)
endif

if (snwgrain) then
rsnwn (:,n) = rsnw (:)
smicen(:,n) = smice(:)
Expand Down