From baaa3f4f9c3242412b12499029b25fec6bdd3714 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Mon, 25 Sep 2023 13:39:18 -0500 Subject: [PATCH 1/2] bring some mushy functions into icepack interface --- columnphysics/icepack_mushy_physics.F90 | 8 ++++---- columnphysics/icepack_therm_itd.F90 | 4 ++-- columnphysics/icepack_therm_shared.F90 | 4 ++-- columnphysics/icepack_therm_vertical.F90 | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/columnphysics/icepack_mushy_physics.F90 b/columnphysics/icepack_mushy_physics.F90 index 466c5d727..3f572e7cd 100644 --- a/columnphysics/icepack_mushy_physics.F90 +++ b/columnphysics/icepack_mushy_physics.F90 @@ -16,7 +16,7 @@ module icepack_mushy_physics conductivity_snow_array, & enthalpy_snow, & enthalpy_brine, & - enthalpy_mush, & + icepack_enthalpy_mush, & enthalpy_mush_liquid_fraction, & enthalpy_of_melting, & temperature_snow, & @@ -262,7 +262,7 @@ end function liquidus_temperature_mush !======================================================================= ! Enthalpy of mush from mush temperature and bulk salinity - function enthalpy_mush(zTin, zSin) result(zqin) + function icepack_enthalpy_mush(zTin, zSin) result(zqin) real(kind=dbl_kind), intent(in) :: & zTin, & ! ice layer temperature (C) @@ -274,14 +274,14 @@ function enthalpy_mush(zTin, zSin) result(zqin) real(kind=dbl_kind) :: & phi ! ice liquid fraction - character(len=*),parameter :: subname='(enthalpy_mush)' + character(len=*),parameter :: subname='(icepack_enthalpy_mush)' phi = icepack_mushy_liquid_fraction(zTin, zSin) zqin = phi * (cp_ocn * rhow - cp_ice * rhoi) * zTin + & rhoi * cp_ice * zTin - (c1 - phi) * rhoi * Lfresh - end function enthalpy_mush + end function icepack_enthalpy_mush !======================================================================= ! Enthalpy of mush from mush temperature and liquid fraction diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index e43e0c08e..b6048dbce 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -49,7 +49,7 @@ module icepack_therm_itd use icepack_itd, only: aggregate_area, shift_ice use icepack_itd, only: column_sum, column_conservation_check use icepack_isotope, only: isoice_alpha, isotope_frac_method - use icepack_mushy_physics, only: liquidus_temperature_mush, enthalpy_mush + use icepack_mushy_physics, only: liquidus_temperature_mush, icepack_enthalpy_mush use icepack_zbgc, only: add_new_ice_bgc use icepack_zbgc, only: lateral_melt_bgc @@ -1579,7 +1579,7 @@ subroutine add_new_ice (ncat, nilyr, & Sprofile(k) = Si0new enddo Ti = min(liquidus_temperature_mush(Si0new/phi_init), Tliquidus_max) - qi0new = enthalpy_mush(Ti, Si0new) + qi0new = icepack_enthalpy_mush(Ti, Si0new) else do k = 1, nilyr Sprofile(k) = salinz(k) diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 8a53a7665..47c29ce86 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -17,7 +17,7 @@ module icepack_therm_shared use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted - use icepack_mushy_physics, only: enthalpy_mush + use icepack_mushy_physics, only: icepack_enthalpy_mush use icepack_mushy_physics, only: temperature_snow use icepack_mushy_physics, only: enthalpy_snow use icepack_mushy_physics, only: icepack_mushy_temperature_mush @@ -317,7 +317,7 @@ subroutine icepack_init_trcr(Tair, Tf, & Ti = Tsfc + slope*(real(k,kind=dbl_kind)-p5) & /real(nilyr,kind=dbl_kind) if (ktherm == 2) then - qin(k) = enthalpy_mush(Ti, Sprofile(k)) + qin(k) = icepack_enthalpy_mush(Ti, Sprofile(k)) else qin(k) = -(rhoi * (cp_ice*(Tprofile(k)-Ti) & + Lfresh*(c1-Tprofile(k)/Ti) - cp_ocn*Tprofile(k))) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 59e18f89f..2d0c20cea 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -45,7 +45,7 @@ module icepack_therm_vertical use icepack_mushy_physics, only: icepack_mushy_temperature_mush use icepack_mushy_physics, only: liquidus_temperature_mush - use icepack_mushy_physics, only: enthalpy_mush, enthalpy_of_melting + use icepack_mushy_physics, only: icepack_enthalpy_mush, enthalpy_of_melting use icepack_aerosol, only: update_aerosol use icepack_isotope, only: update_isotope @@ -1257,7 +1257,7 @@ subroutine thickness_changes (nilyr, nslyr, & if (ktherm == 2) then - qbotm = enthalpy_mush(Tbot, sss) + qbotm = icepack_enthalpy_mush(Tbot, sss) qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy) qbot0 = qbotm - qbotp From 1f0b307fa3d171f1ccf1af7ad4708771ae5e5eaf Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Mon, 25 Sep 2023 17:51:25 -0500 Subject: [PATCH 2/2] elevate snow/mush enthalpy and salinity profile to interface --- columnphysics/icepack_intfc.F90 | 4 +- columnphysics/icepack_mushy_physics.F90 | 8 ++-- columnphysics/icepack_therm_mushy.F90 | 4 +- columnphysics/icepack_therm_shared.F90 | 54 ++++++++++++++----------- 4 files changed, 39 insertions(+), 31 deletions(-) diff --git a/columnphysics/icepack_intfc.F90 b/columnphysics/icepack_intfc.F90 index d76be4ff8..0d8be0b31 100644 --- a/columnphysics/icepack_intfc.F90 +++ b/columnphysics/icepack_intfc.F90 @@ -112,10 +112,12 @@ module icepack_intfc use icepack_therm_shared , only: icepack_snow_temperature use icepack_therm_shared , only: icepack_liquidus_temperature use icepack_therm_shared , only: icepack_sea_freezing_temperature - use icepack_therm_shared , only: icepack_enthalpy_snow use icepack_therm_shared , only: icepack_init_thermo + use icepack_therm_shared , only: icepack_salinity_profile use icepack_therm_shared , only: icepack_init_trcr + use icepack_mushy_physics , only: icepack_enthalpy_snow + use icepack_mushy_physics , only: icepack_enthalpy_mush use icepack_mushy_physics , only: icepack_mushy_density_brine use icepack_mushy_physics , only: icepack_mushy_liquid_fraction use icepack_mushy_physics , only: icepack_mushy_temperature_mush diff --git a/columnphysics/icepack_mushy_physics.F90 b/columnphysics/icepack_mushy_physics.F90 index 3f572e7cd..03700d56c 100644 --- a/columnphysics/icepack_mushy_physics.F90 +++ b/columnphysics/icepack_mushy_physics.F90 @@ -14,7 +14,7 @@ module icepack_mushy_physics public :: & conductivity_mush_array, & conductivity_snow_array, & - enthalpy_snow, & + icepack_enthalpy_snow, & enthalpy_brine, & icepack_enthalpy_mush, & enthalpy_mush_liquid_fraction, & @@ -143,7 +143,7 @@ end subroutine conductivity_snow_array !======================================================================= ! Enthalpy of snow from snow temperature - function enthalpy_snow(zTsn) result(zqsn) + function icepack_enthalpy_snow(zTsn) result(zqsn) real(kind=dbl_kind), intent(in) :: & zTsn ! snow layer temperature (C) @@ -151,11 +151,11 @@ function enthalpy_snow(zTsn) result(zqsn) real(kind=dbl_kind) :: & zqsn ! snow layer enthalpy (J m-3) - character(len=*),parameter :: subname='(enthalpy_snow)' + character(len=*),parameter :: subname='(icepack_enthalpy_snow)' zqsn = -rhos * (-cp_ice * zTsn + Lfresh) - end function enthalpy_snow + end function icepack_enthalpy_snow !======================================================================= ! Temperature of snow from the snow enthalpy diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index 999c92af7..e182e7db2 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -10,7 +10,7 @@ module icepack_therm_mushy use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode 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_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, enthalpy_snow + 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 use icepack_mushy_physics, only: temperature_snow, temperature_mush_liquid_fraction @@ -1802,7 +1802,7 @@ subroutine picard_final(lsnow, & if (lsnow) then do k = 1, nslyr - zqsn(k) = enthalpy_snow(zTsn(k)) + zqsn(k) = icepack_enthalpy_snow(zTsn(k)) enddo ! k endif ! lsnow diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 47c29ce86..b3a752f1b 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -19,7 +19,6 @@ module icepack_therm_shared use icepack_mushy_physics, only: icepack_enthalpy_mush use icepack_mushy_physics, only: temperature_snow - use icepack_mushy_physics, only: enthalpy_snow use icepack_mushy_physics, only: icepack_mushy_temperature_mush use icepack_mushy_physics, only: liquidus_temperature_mush @@ -30,12 +29,12 @@ module icepack_therm_shared surface_heat_flux, & dsurface_heat_flux_dTsf, & icepack_init_thermo, & + icepack_salinity_profile, & icepack_init_trcr, & icepack_ice_temperature, & icepack_snow_temperature, & icepack_liquidus_temperature, & icepack_sea_freezing_temperature, & - icepack_enthalpy_snow, & adjust_enthalpy real (kind=dbl_kind), parameter, public :: & @@ -226,10 +225,6 @@ subroutine icepack_init_thermo(nilyr, sprofile) !autodocument_end - real (kind=dbl_kind), parameter :: & - nsal = 0.407_dbl_kind, & - msal = 0.573_dbl_kind - integer (kind=int_kind) :: k ! ice layer index real (kind=dbl_kind) :: zn ! normalized ice thickness @@ -254,7 +249,8 @@ subroutine icepack_init_thermo(nilyr, sprofile) do k = 1, nilyr zn = (real(k,kind=dbl_kind)-p5) / & real(nilyr,kind=dbl_kind) - sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn)))) +! sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn)))) + sprofile(k)=icepack_salinity_profile(zn) sprofile(k) = max(sprofile(k), min_salin) enddo ! k sprofile(nilyr+1) = saltmax @@ -267,6 +263,33 @@ subroutine icepack_init_thermo(nilyr, sprofile) end subroutine icepack_init_thermo +!======================================================================= +!autodocument_start icepack_salinity_profile +! Initial salinity profile +! +! authors: C. M. Bitz, UW +! William H. Lipscomb, LANL + + function icepack_salinity_profile(zn) result(salinity) + + real(kind=dbl_kind), intent(in) :: & + zn ! depth + + real(kind=dbl_kind) :: & + salinity ! initial salinity profile + +!autodocument end + + real (kind=dbl_kind), parameter :: & + nsal = 0.407_dbl_kind, & + msal = 0.573_dbl_kind + + character(len=*),parameter :: subname='(icepack_init_thermo)' + + salinity = (saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn)))) + + end function icepack_salinity_profile + !======================================================================= !autodocument_start icepack_init_trcr ! @@ -443,23 +466,6 @@ function icepack_snow_temperature(qin) result(Tsn) end function icepack_snow_temperature -!======================================================================= -!autodocument_start icepack_enthalpy_snow -! compute snow enthalpy - - function icepack_enthalpy_snow(zTsn) result(qsn) - - real(kind=dbl_kind), intent(in) :: zTsn - real(kind=dbl_kind) :: qsn - -!autodocument_end - - character(len=*),parameter :: subname='(icepack_enthalpy_snow)' - - qsn = enthalpy_snow(zTsn) - - end function icepack_enthalpy_snow - !======================================================================= ! ! Conserving energy, compute the new enthalpy of equal-thickness ice