diff --git a/components/mpas-seaice/src/Registry.xml b/components/mpas-seaice/src/Registry.xml
index 2b90267cc41d..a8767cdc77bc 100644
--- a/components/mpas-seaice/src/Registry.xml
+++ b/components/mpas-seaice/src/Registry.xml
@@ -4345,6 +4345,11 @@
dimensions="nCells Time"
units="1"
/>
+
0 .and. k <= natmiter)
ustar_prev = ustar
diff --git a/components/mpas-seaice/src/column/ice_colpkg.F90 b/components/mpas-seaice/src/column/ice_colpkg.F90
index 4eb4aa875aa4..e1ba360d7646 100644
--- a/components/mpas-seaice/src/column/ice_colpkg.F90
+++ b/components/mpas-seaice/src/column/ice_colpkg.F90
@@ -3640,7 +3640,7 @@ subroutine colpkg_atm_boundary(sfctype, &
Uref)
use ice_atmo, only: atmo_boundary_const, atmo_boundary_layer
- use ice_constants_colpkg, only: c0
+ use ice_constants_colpkg, only: c0, c1
character (len=3), intent(in) :: &
sfctype ! ice or ocean
@@ -3689,10 +3689,13 @@ subroutine colpkg_atm_boundary(sfctype, &
worku = uvel
endif
! should this be for vvel,workv?
- if (present(uvel)) then
- worku = uvel
+ if (present(vvel)) then
+ workv = vvel
endif
+ ! NJ keeps icepack/colpkg BFB when atmbndy = 'constant'
+ Cdn_atm_ratio_n = c1
+
if (trim(atmbndy) == 'constant') then
call atmo_boundary_const (sfctype, calc_strair, &
uatm, vatm, &
diff --git a/components/mpas-seaice/src/column/ice_shortwave.F90 b/components/mpas-seaice/src/column/ice_shortwave.F90
index a44f06f2612a..9e74987a03e8 100644
--- a/components/mpas-seaice/src/column/ice_shortwave.F90
+++ b/components/mpas-seaice/src/column/ice_shortwave.F90
@@ -60,6 +60,9 @@ module ice_shortwave
hpmin = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
hp0 = 0.200_dbl_kind ! pond depth below which transition to bare ice
+ real (kind=dbl_kind), parameter :: &
+ exp_argmax = c10 ! maximum argument of exponential
+
real (kind=dbl_kind) :: &
exp_min ! minimum exponential value
@@ -926,16 +929,11 @@ subroutine run_dEdd(dt, tr_aero, &
logical (kind=log_kind) :: &
linitonly ! local initonly value
- real (kind=dbl_kind), parameter :: &
- argmax = c10 ! maximum argument of exponential
-
linitonly = .false.
if (present(initonly)) then
linitonly = initonly
endif
- exp_min = exp(-argmax)
-
! cosine of the zenith angle
call compute_coszen (tlat, tlon, &
calendar_type, days_per_year, &
@@ -2397,13 +2395,13 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, &
! get grain size index:
! works for 25 < snw_rds < 1625 um:
- if (tmp_gs < 125) then
- tmp1 = tmp_gs/50
+ if (tmp_gs < 125.0_dbl_kind) then
+ tmp1 = tmp_gs/50.0_dbl_kind
k_bcini(k) = nint(tmp1)
- elseif (tmp_gs < 175) then
+ elseif (tmp_gs < 175.0_dbl_kind) then
k_bcini(k) = 2
else
- tmp1 = (tmp_gs/250)+2
+ tmp1 = (tmp_gs/250.0_dbl_kind)+c2
k_bcini(k) = nint(tmp1)
endif
else ! use the largest snow grain size for ice
@@ -3347,7 +3345,10 @@ subroutine solution_dEdd &
tdr , & ! tdir for gaussian integration
smr , & ! accumulator for rdif gaussian integration
smt ! accumulator for tdif gaussian integration
-
+
+ real (kind=dbl_kind) :: &
+ exp_min ! minimum exponential value
+
! Delta-Eddington solution expressions
alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
@@ -3385,7 +3386,7 @@ subroutine solution_dEdd &
! value below the fresnel level, i.e. the cosine solar zenith
! angle below the fresnel level for the refracted solar beam:
mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
-
+
! compute level of fresnel refraction
! if ponded sea ice, fresnel level is the top of the pond.
kfrsnl = 0
@@ -3434,7 +3435,9 @@ subroutine solution_dEdd &
! non-refracted beam instead
if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
- extins = max(exp_min, exp(-lm*ts))
+ !extins = max(exp_min, exp(-lm*ts))
+ exp_min = min(exp_argmax,lm*ts)
+ extins = exp(-exp_min)
ne = n(ue,extins)
! first calculation of rdif, tdif using Delta-Eddington formulas
@@ -3443,7 +3446,9 @@ subroutine solution_dEdd &
tdif_a(k) = c4*ue/ne
! evaluate rdir,tdir for direct beam
- trnlay(k) = max(exp_min, exp(-ts/mu0n))
+ !trnlay(k) = max(exp_min, exp(-ts/mu0n))
+ exp_min = min(exp_argmax,ts/mu0n)
+ trnlay(k) = exp(-exp_min)
alp = alpha(ws,mu0n,gs,lm)
gam = agamm(ws,mu0n,gs,lm)
apg = alp + gam
@@ -3464,7 +3469,9 @@ subroutine solution_dEdd &
mu = gauspt(ng)
gwt = gauswt(ng)
swt = swt + mu*gwt
- trn = max(exp_min, exp(-ts/mu))
+ !trn = max(exp_min, exp(-ts/mu))
+ exp_min = min(exp_argmax,ts/mu)
+ trn = exp(-exp_min)
alp = alpha(ws,mu,gs,lm)
gam = agamm(ws,mu,gs,lm)
apg = alp + gam
@@ -4648,13 +4655,13 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
! get grain size index:
! works for 25 < snw_rds < 1625 um:
- if (tmp_gs < 125) then
- tmp1 = tmp_gs/50
+ if (tmp_gs < 125._dbl_kind) then
+ tmp1 = tmp_gs/50._dbl_kind
k_bcini(k) = nint(tmp1)
- elseif (tmp_gs < 175) then
+ elseif (tmp_gs < 175._dbl_kind) then
k_bcini(k) = 2
else
- tmp1 = (tmp_gs/250)+2
+ tmp1 = (tmp_gs/250._dbl_kind) + c2
k_bcini(k) = nint(tmp1)
endif
else ! use the largest snow grain size for ice
diff --git a/components/mpas-seaice/src/icepack b/components/mpas-seaice/src/icepack
index a9fc4fafa834..6ca316f0306c 160000
--- a/components/mpas-seaice/src/icepack
+++ b/components/mpas-seaice/src/icepack
@@ -1 +1 @@
-Subproject commit a9fc4fafa83460a8b8582e5a016d97259f0bada8
+Subproject commit 6ca316f0306c33cb68d4b8b76342e91e73da5a86
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_column.F b/components/mpas-seaice/src/shared/mpas_seaice_column.F
index 6debe1e7c7db..dae49de50f33 100644
--- a/components/mpas-seaice/src/shared/mpas_seaice_column.F
+++ b/components/mpas-seaice/src/shared/mpas_seaice_column.F
@@ -6256,8 +6256,12 @@ subroutine seaice_init_column_constants()
seaiceIceSnowEmissivity, &
seaiceFreshWaterFreezingPoint, &
seaiceAirSpecificHeat, &
+ seaiceWaterVaporSpecificHeat, &
+ seaiceSeaWaterSpecificHeat, &
+ seaiceLatentHeatVaporization, &
seaiceLatentHeatSublimation, &
seaiceLatentHeatMelting, &
+ seaiceReferenceSalinity, &
seaiceOceanAlbedo, &
seaiceVonKarmanConstant, &
seaiceIceSurfaceRoughness, &
@@ -6281,8 +6285,12 @@ subroutine seaice_init_column_constants()
emissivity, &
Tffresh, &
cp_air, &
+ cp_wv, &
+ cp_ocn, &
+ Lvap, &
Lsub, &
Lfresh, &
+ ice_ref_salinity, &
Pstar, &
Cstar, &
dragio, &
@@ -6304,8 +6312,12 @@ subroutine seaice_init_column_constants()
seaiceIceSnowEmissivity = emissivity
seaiceFreshWaterFreezingPoint = Tffresh
seaiceAirSpecificHeat = cp_air
+ seaiceWaterVaporSpecificHeat = cp_wv
+ seaiceSeaWaterSpecificHeat = cp_ocn
+ seaiceLatentHeatVaporization = Lvap
seaiceLatentHeatSublimation = Lsub
seaiceLatentHeatMelting = Lfresh
+ seaiceReferenceSalinity = ice_ref_salinity
seaiceOceanAlbedo = albocn
seaiceVonKarmanConstant = vonkar
seaiceIceSurfaceRoughness = iceruf
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
index 09aae35ad60a..2adf76da3aac 100644
--- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
+++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
@@ -1415,6 +1415,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
oceanStressCellV, &
freezingMeltingPotential, &
lateralIceMeltFraction, &
+ lateralIceMeltRate, &
lateralHeatFlux, &
snowfallRate, &
rainfallRate, &
@@ -1660,6 +1661,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
call MPAS_pool_get_array(drag, "dragFloeSeparation", dragFloeSeparation)
call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltFraction", lateralIceMeltFraction)
+ call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltRate", lateralIceMeltRate)
call MPAS_pool_get_array(melt_growth_rates, "lateralHeatFlux", lateralHeatFlux)
call MPAS_pool_get_array(melt_growth_rates, "surfaceIceMelt", surfaceIceMelt)
call MPAS_pool_get_array(melt_growth_rates, "surfaceIceMeltCategory", surfaceIceMeltCategory)
@@ -1861,6 +1863,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
Tsnice=snowIceInterfaceTemperature(iCell), &
frzmlt=freezingMeltingPotential(iCell), &
rside=lateralIceMeltFraction(iCell), &
+ wlat=lateralIceMeltRate(iCell), &
fside=lateralHeatFlux(iCell), &
fsnow=snowfallRate(iCell), &
frain=rainfallRate(iCell), &
@@ -2025,6 +2028,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
call mpas_log_write("oceanHeatFluxIceBottom: $r", messageType=MPAS_LOG_ERR, realArgs=(/oceanHeatFluxIceBottom(iCell)/))
call mpas_log_write("freezingMeltingPotential: $r", messageType=MPAS_LOG_ERR, realArgs=(/freezingMeltingPotential(iCell)/))
call mpas_log_write("lateralIceMeltFraction: $r", messageType=MPAS_LOG_ERR, realArgs=(/lateralIceMeltFraction(iCell)/))
+ call mpas_log_write("lateralIceMeltRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/lateralIceMeltRate(iCell)/))
call mpas_log_write("snowfallRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/snowfallRate(iCell)/))
call mpas_log_write("rainfallRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/rainfallRate(iCell)/))
call mpas_log_write("pondFreshWaterFlux: $r", messageType=MPAS_LOG_ERR, realArgs=(/pondFreshWaterFlux(iCell)/))
@@ -2189,6 +2193,7 @@ subroutine column_itd_thermodynamics(domain, clock)
seaFreezingTemperature, &
seaSurfaceSalinity, &
lateralIceMeltFraction, &
+ lateralIceMeltRate, &
lateralIceMelt, &
freezingMeltingPotential, &
frazilFormation, &
@@ -2323,6 +2328,7 @@ subroutine column_itd_thermodynamics(domain, clock)
call MPAS_pool_get_array(ocean_fluxes, "oceanHeatFlux", oceanHeatFlux)
call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltFraction", lateralIceMeltFraction)
+ call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltRate", lateralIceMeltRate)
call MPAS_pool_get_array(melt_growth_rates, "lateralIceMelt", lateralIceMelt)
call MPAS_pool_get_array(melt_growth_rates, "frazilFormation", frazilFormation)
call MPAS_pool_get_array(melt_growth_rates, "frazilGrowthDiagnostic", frazilGrowthDiagnostic)
@@ -2422,6 +2428,7 @@ subroutine column_itd_thermodynamics(domain, clock)
seaSurfaceSalinity(iCell), &
initialSalinityProfile(:,iCell), &
lateralIceMeltFraction(iCell), &
+!in icepack, not colpkg lateralIceMeltRate(iCell), &
lateralIceMelt(iCell), &
freezingMeltingPotential(iCell), &
frazilFormation(iCell), &